dmri_connections_regions_calc

  CALCULATE DMRI CONNECTIONS BETWEEN REGIONS 
 [Usage]
  [CC,DD,xyz0,Keys0]=dmri_connections_regions_calc(dmrifile, brainfile, areafile);

 [Input]
  dmrifile  : connectivity matrix obtained from DMRI processor
  brainfile : brainfile of VBMEG
  areafile  : areafile of coarse region parcellation such as AAL and
              brodmann. AAL is recommended.
 [Output]
  CC        : region connectivity matrix (number of connections average over region vertices)
  DD        : region distance matrix in unit meter  (distance average over region vertices). 

 2015/09/08 O.Yamashita

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function [CC,DD,xyz0,Keys0]=dmri_connections_regions_calc(dmrifile, brainfile, areafile, parcelfile)
0002 %  CALCULATE DMRI CONNECTIONS BETWEEN REGIONS
0003 % [Usage]
0004 %  [CC,DD,xyz0,Keys0]=dmri_connections_regions_calc(dmrifile, brainfile, areafile);
0005 %
0006 % [Input]
0007 %  dmrifile  : connectivity matrix obtained from DMRI processor
0008 %  brainfile : brainfile of VBMEG
0009 %  areafile  : areafile of coarse region parcellation such as AAL and
0010 %              brodmann. AAL is recommended.
0011 % [Output]
0012 %  CC        : region connectivity matrix (number of connections average over region vertices)
0013 %  DD        : region distance matrix in unit meter  (distance average over region vertices).
0014 %
0015 % 2015/09/08 O.Yamashita
0016 %
0017 % Copyright (C) 2011, ATR All Rights Reserved.
0018 % License : New BSD License(see VBMEG_LICENSE.txt)
0019 
0020 if nargin < 4
0021     parcelfile = [];
0022 end
0023 
0024 
0025 %---------------------
0026 % brain model
0027 %---------------------
0028 
0029 B=load(brainfile);
0030 
0031 V = B.V .* 1000; % [m] => [mm]
0032 F0 = B.F.F3;
0033 
0034 %------------------------------
0035 %  connectivity
0036 %------------------------------
0037 
0038 % IX : parcel indicies to region indices
0039 [tmp,IX,xyz,Keys] = inner_parcels_to_regions(areafile,dmrifile,parcelfile,V);
0040 
0041 % remove empty regions
0042 [IX0,xyz0,Keys0] = inner_remove_null_regions(IX,xyz,Keys);  
0043 
0044 % connectivity matrix, distance matrix
0045 [C,D] = inner_get_connectivity_matrix(dmrifile,parcelfile);
0046 
0047 % region connectivity
0048 [CC,DD] = inner_get_region_connectivity(C,D,IX0);
0049     
0050 
0051 %------------------------------------------------------
0052 % inner functions
0053 %
0054 
0055 function [C,IX,xyz,Keys] = inner_parcels_to_regions(areafile,dmrifile,parcelfile,V)
0056 
0057 if isempty(parcelfile)  % new format
0058     AC = load(dmrifile);
0059     C = AC.connectivity_matrix;
0060     ix = AC.vbmeg_area_ix;
0061 else
0062     AC = load(dmrifile);
0063     P = load(parcelfile);
0064     ix = P.ix_area_mod;
0065     C = full(AC.Ind);
0066 
0067 end
0068 
0069 % Corse representaion of connectivity matrix (require area file)
0070 A = load(areafile);
0071 Narea = length(A.Area);
0072 
0073 for aa = 1 : Narea
0074     tmpArea = A.Area{aa};
0075     xyz(aa,:) = mean(V(tmpArea.Iextract,:),1);  % center of region
0076     [tmp,jx] = intersect(ix, tmpArea.Iextract);
0077     IX{aa} = jx;
0078     Keys{aa} = tmpArea.key;
0079 end
0080 
0081 function [IX0,xyz0,Keys0] = inner_remove_null_regions(IX,xyz,Keys)
0082 % erase empty area
0083 nn = 1;
0084 for aa = 1 : size(xyz,1)
0085     if ~isempty(IX{aa})
0086         xyz0(nn,:) = xyz(aa,:);
0087         IX0{nn} = IX{aa};
0088         Keys0{nn} = Keys{aa};
0089         nn = nn + 1;
0090     end
0091 end
0092 
0093 function [CC,DD] = inner_get_region_connectivity(C,D,IX0)
0094 
0095 Narea0 = length(IX0);
0096 CC = zeros(Narea0, Narea0);
0097 DD = zeros(Narea0, Narea0);
0098 for aa1 = 1 : Narea0
0099     for aa2 = 1 : Narea0
0100         CC(aa1,aa2) = mean(mean(C(IX0{aa1},IX0{aa2})));
0101         DD(aa1,aa2) = mean(mean(D(IX0{aa1},IX0{aa2})));
0102         
0103     end
0104 end
0105 
0106 
0107 
0108 function [C,D] = inner_get_connectivity_matrix(dmrifile,parcelfile)
0109 
0110 if isempty(parcelfile)  % new format
0111     AC = load(dmrifile);
0112     C = AC.connectivity_matrix;
0113     
0114     D = zeros(size(AC.connectivity_matrix));
0115     tmp = AC.delay_ms_matrix ;
0116     ix = find(tmp > 0);
0117     D(ix) = (tmp(ix) - AC.dmri_parm.tau)/1000*AC.dmri_parm.v;   % distance matrix in s*(m/s) = m
0118     
0119 else
0120     AC = load(dmrifile);
0121     P = load(parcelfile);
0122     ix = P.ix_area_mod;
0123     C = full(AC.Ind);
0124    
0125     D = zeros(size(AC.connectivity_matrix));
0126     tmp = AC.Delta;
0127     ix = find(tmp > 0);
0128     D(ix) = (tmp(ix) - AC.dmri_parm.tau)/1000*AC.dmri_parm.v;   % distance matrix in s/(m/s) = m
0129     
0130     
0131 end