dmri_cortex_parcel_from_area

 Get indices of parcellated brain
 [Usage]
    dmri_cortex_parcel_from_area(brain_file, fs_info_file, parcel_area_file, parcel_file);

 [Input]
       brain_file : VBMEG BRAIN-MAT file(.brain.mat)
     fs_info_file : Freesurfer information file.
 parcel_area_file : Parcellation is done with the area file(.area.mat).
      parcel_file : Parcel file.

 Original code is vb_original_normal_statics.m
 2012/12/02 M.Fukushima
 2017/10/25 Y.Takeda
  Modified for parceling cortex according to area file

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

0001 function dmri_cortex_parcel_from_area(brain_file, fs_info_file, parcel_area_file, parcel_file)
0002 % Get indices of parcellated brain
0003 % [Usage]
0004 %    dmri_cortex_parcel_from_area(brain_file, fs_info_file, parcel_area_file, parcel_file);
0005 %
0006 % [Input]
0007 %       brain_file : VBMEG BRAIN-MAT file(.brain.mat)
0008 %     fs_info_file : Freesurfer information file.
0009 % parcel_area_file : Parcellation is done with the area file(.area.mat).
0010 %      parcel_file : Parcel file.
0011 %
0012 % Original code is vb_original_normal_statics.m
0013 % 2012/12/02 M.Fukushima
0014 % 2017/10/25 Y.Takeda
0015 %  Modified for parceling cortex according to area file
0016 %
0017 % Copyright (C) 2011, ATR All Rights Reserved.
0018 % License : New BSD License(see VBMEG_LICENSE.txt)
0019 
0020 %
0021 % --- Previous check
0022 %
0023 if nargin ~= 4
0024     error('Please check input arguments.');
0025 end
0026 if exist(brain_file, 'file') ~= 2
0027     error('brain_file not found.');
0028 end
0029 if exist(fs_info_file, 'file') ~= 2
0030     error('fs_info_file not found.');
0031 end
0032 
0033 %
0034 % --- Main Procedure
0035 %
0036 disp('Cortex parcellation..');
0037 start = tic;
0038 
0039 % Max radius of neighbor search
0040 Rmax = 0.004;    % 4mm
0041 
0042 %
0043 % --- Load original brain
0044 %
0045 load(fs_info_file);
0046 
0047 brain_parm.FS_left_file  = FS_wm.lh_smoothwm_asc_file;
0048 brain_parm.FS_right_file = FS_wm.rh_smoothwm_asc_file;
0049 
0050 vb_disp('Load original brain')
0051 [V0L,F0L,n0L,V0R,F0R,n0R] = vb_load_orig_brain_surf(brain_parm);
0052 %
0053 % --- SPM cordinate in [m] , Normal vector is outward.
0054 %
0055 vb_disp('Change coordinate to [m] ');
0056 
0057 V0L = V0L/1000;
0058 V0R = V0R/1000;
0059 
0060 %
0061 % --- Set normal direction outward
0062 %     normal vector for disconnected vertex is set to zero
0063 tic
0064 vb_disp_nonl('Set normal direction outward: ')
0065 [F0L, V0L, xxL] = vb_out_normal_surf(F0L,V0L);
0066 
0067 xxsig = sign( sum( n0L .* xxL , 2) );
0068 n0L = n0L .* repmat(xxsig, [1 3]);
0069 
0070 [F0R, V0R, xxR] = vb_out_normal_surf(F0R,V0R);
0071 
0072 xxsig = sign( sum( n0R .* xxR , 2) );
0073 n0R = n0R .* repmat(xxsig, [1 3]);
0074 vb_disp(sprintf('%f[sec]',toc));
0075 %vb_ptime(toc)
0076 
0077 % Original vertex index for cortex
0078 load(brain_file, 'subj');
0079 [tmp1, tmp2, BV_index] = vb_load_cortex_info(brain_file, 'subj');
0080 
0081 NV0L  = size(V0L,1);
0082 indxL = BV_index.Left;
0083 indxR = BV_index.Right;
0084 
0085 if min(indxR) > NV0L,
0086     indxR = indxR - NV0L;
0087 end
0088 
0089 %
0090 % Search neighbor points along cortex sheet
0091 % and calculate normal direction statistics
0092 %
0093 
0094 %
0095 % --- Make neighbor member list of original brain (Left)
0096 %
0097 tic
0098 vb_disp_nonl('Find neighbor index (Left): ')
0099 [xxDL,xxFL]   = vb_next_distance( F0L, V0L ); 
0100 vb_disp(sprintf('%f[sec]',toc));
0101 %vb_ptime(toc)
0102 
0103 vb_disp_nonl('Make neighbor member list of original brain (Left): ')
0104 tic
0105 Llist  = vb_find_near_member(xxFL, xxDL, indxL, Rmax);
0106 vb_disp(sprintf('%f[sec]',toc));
0107 %vb_ptime(toc)
0108 
0109 %
0110 % --- Make neighbor member list of original brain (Right)
0111 %
0112 tic
0113 vb_disp_nonl('Find neighbor index (Right): ')
0114 [xxDR,xxFR]   = vb_next_distance( F0R, V0R ); 
0115 vb_disp(sprintf('%f[sec]',toc));
0116 %vb_ptime(toc)
0117 
0118 vb_disp_nonl('Make neighbor member list of original brain (Right): ')
0119 tic
0120 Rlist  = vb_find_near_member(xxFR, xxDR, indxR, Rmax);
0121 vb_disp(sprintf('%f[sec]',toc));
0122 %vb_ptime(toc)
0123 
0124 %% Whole Brain Parcellation
0125 % Search radious for parcellation (as large as possible)
0126 Rmax = 0.040; % 4cm
0127 
0128 % Parcel cortex according to all the vertices in brain model
0129 vbmeg_area_ix=1:length(indxL)+length(indxR);
0130 
0131 Nroot_l = length(indxL);
0132 tmp = find(vbmeg_area_ix<=Nroot_l);
0133 half = tmp(end);
0134 
0135 % Parcellation
0136 [lh_id1] = vb_find_near_member...
0137   (xxFL, xxDL, indxL(vbmeg_area_ix(1:half)), Rmax);
0138 [rh_id1] = vb_find_near_member...
0139   (xxFR, xxDR, indxR(vbmeg_area_ix(half+1:end)-Nroot_l), Rmax);
0140 
0141 % Combine parcels according to area file
0142 load(parcel_area_file,'Area')
0143 Narea=length(Area);
0144 al=1;
0145 ar=1;
0146 for area=1:Narea
0147     ix=Area{area}.Iextract;
0148     ix1=[];
0149     if ix(1)<=Nroot_l% Left
0150         for i=1:length(ix)
0151             ix1=[ix1,lh_id1{ix(i)}];
0152         end
0153         lh_id2{al}=sort(ix1);
0154         area_num_l1(al)=area;
0155         al=al+1;
0156     else% Right
0157         ix=ix-Nroot_l;
0158         for i=1:length(ix)
0159             ix1=[ix1,rh_id1{ix(i)}];
0160         end
0161         rh_id2{ar}=sort(ix1);
0162         area_num_r1(ar)=area;
0163         ar=ar+1;
0164     end
0165 end
0166 
0167 % Make (lh/rh)_cortex_id1 which contains only cortex index (remove subcortex elements)
0168 lh_cortex_id1 = cell(length(lh_id2), 1);
0169 rh_cortex_id1 = cell(length(rh_id2), 1);
0170 for n=1:length(lh_id2)
0171     lh_cortex_id1{n} = lh_id2{n};
0172     [tmp, ind] = intersect(lh_id2{n}, FS_wm.lh_subcortex_index);
0173     lh_cortex_id1{n}(ind) = [];
0174 end
0175 for n=1:length(rh_id2)
0176     rh_cortex_id1{n} = rh_id2{n};
0177     [tmp, ind] = intersect(rh_id2{n}, FS_wm.rh_subcortex_index);
0178     rh_cortex_id1{n}(ind) = [];
0179 end
0180 
0181 % Extract parcels which have no or small subcortex
0182 th_p=0.1;% Threshold for portion of subcortex
0183 p_subcortex=zeros(length(lh_cortex_id1),1);
0184 for n=1:length(lh_cortex_id1)
0185     p_subcortex(n)=(length(lh_id2{n})-length(lh_cortex_id1{n}))/length(lh_id2{n});
0186 end
0187 ix=find(p_subcortex<th_p);
0188 lh_id=lh_id2(ix);
0189 lh_cortex_id=lh_cortex_id1(ix);
0190 area_num_l=area_num_l1(ix);
0191 p_subcortex=zeros(length(rh_cortex_id1),1);
0192 for n=1:length(rh_cortex_id1)
0193     p_subcortex(n)=(length(rh_id2{n})-length(rh_cortex_id1{n}))/length(rh_id2{n});
0194 end
0195 ix=find(p_subcortex<th_p);
0196 rh_id=rh_id2(ix);
0197 rh_cortex_id=rh_cortex_id1(ix);
0198 area_num_r=area_num_r1(ix);
0199 
0200 % Make vbmeg_area_ix
0201 Area_cortex=Area([area_num_l,area_num_r]);
0202 Narea=length(Area_cortex);
0203 V=vb_load_cortex(brain_file,'Inflate');
0204 vbmeg_area_ix=zeros(Narea,1);
0205 for area=1:Narea
0206     ix=Area_cortex{area}.Iextract;
0207     center=mean(V(ix,:),1);
0208     d=sum((V(ix,:)-repmat(center,length(ix),1)).^2,2);
0209     [~,b]=min(d);
0210     vbmeg_area_ix(area)=ix(b);
0211 end
0212 
0213 % Save index
0214 save(parcel_file, 'Area_cortex', 'vbmeg_area_ix', 'lh_id', 'rh_id', 'lh_cortex_id', 'rh_cortex_id', 'brain_file');
0215 
0216 % Create membership matrix(parcels to brain model(V))
0217 [membershipmat] = make_membershipmat(brain_file, parcel_file);
0218 save(parcel_file, '-append', 'membershipmat');
0219 
0220 disp(sprintf('\nSaved:%s', parcel_file));
0221 
0222 toc(start);