Make cortex data from MNI standard brain for group analysis
1. Map MNI standard brain to subject brain
2. Find nearest point from subject cortical surface
extracted by FreeSurfer
[V,F,xx,ORG_index,Vinfo] = vb_make_brain_mni(brain_parm,reg_mode);
--- Input
brain_parm : structure with following field
--- fields of brain_parm
.FS_left_file = [fs_dir '/bem/lh.smoothwm.asc'];
.FS_right_file = [fs_dir '/bem/rh.smoothwm.asc'];
.FS_left_label_file = [fs_dir '/label/lh.cortex.label'];
.FS_right_label_file = [fs_dir '/label/rh.cortex.label'];
if reg_mode=='SPM'
.spm_normalization_file: SPM normalization file
(_sn.mat)
if reg_mode=='FS'
.FS_left_sphere_file
.FS_right_sphere_file
.FS_sphere_key
--- Output
V : Cortical vertex point cordinate (SPM_Right_m) [Nvertex, 3]
頂点位置 （ 頂点の数, 3）
xx : Normal vector to cortical surface [Nvertex, 3]
頂点の法線方向と位置（ 頂点の数, 3）
F : Patch index structure
面(3角形)を構成する３つの頂点番号 ( 面の数, 3)
.F3R : Right cortex
.F3L : Left cortex
.F3 : Left + Right
BV_index : original vertex index in BV corresponding to brain
.Left : Left brain
.Right : Right brain
Vinfo : Vertex dimension structure
.Ndipole : # of vertex
.NdipoleL : # of vertex in Left cortex
.Npatch : # of patch
.Coord = 'SPM_Right_m';
[history]
Ver.2.0 New version for group analysis
2015-12-27 M. Sato
2017-01-30 Y. Takeda vb_find_nearest_point->vb_find_nearest_point_no_overlap
2017-03-06 M. Sato create Vext, Fext.
Copyright (C) 2011, ATR All Rights Reserved.
License : New BSD License(see VBMEG_LICENSE.txt)0001 function [V,F,xx,BV_index,Vinfo,Vext,Fext] = ... 0002 vb_make_brain_mni(brain_parm,reg_mode) 0003 % Make cortex data from MNI standard brain for group analysis 0004 % 1. Map MNI standard brain to subject brain 0005 % 2. Find nearest point from subject cortical surface 0006 % extracted by FreeSurfer 0007 % 0008 % [V,F,xx,ORG_index,Vinfo] = vb_make_brain_mni(brain_parm,reg_mode); 0009 % 0010 %--- Input 0011 % 0012 % brain_parm : structure with following field 0013 % --- fields of brain_parm 0014 % .FS_left_file = [fs_dir '/bem/lh.smoothwm.asc']; 0015 % .FS_right_file = [fs_dir '/bem/rh.smoothwm.asc']; 0016 % 0017 % .FS_left_label_file = [fs_dir '/label/lh.cortex.label']; 0018 % .FS_right_label_file = [fs_dir '/label/rh.cortex.label']; 0019 % 0020 % if reg_mode=='SPM' 0021 % .spm_normalization_file: SPM normalization file 0022 % (_sn.mat) 0023 % if reg_mode=='FS' 0024 % .FS_left_sphere_file 0025 % .FS_right_sphere_file 0026 % .FS_sphere_key 0027 % 0028 % 0029 %--- Output 0030 % 0031 % V : Cortical vertex point cordinate (SPM_Right_m) [Nvertex, 3] 0032 % 頂点位置 （ 頂点の数, 3） 0033 % xx : Normal vector to cortical surface [Nvertex, 3] 0034 % 頂点の法線方向と位置（ 頂点の数, 3） 0035 % F : Patch index structure 0036 % 面(3角形)を構成する３つの頂点番号 ( 面の数, 3) 0037 % .F3R : Right cortex 0038 % .F3L : Left cortex 0039 % .F3 : Left + Right 0040 % BV_index : original vertex index in BV corresponding to brain 0041 % .Left : Left brain 0042 % .Right : Right brain 0043 % Vinfo : Vertex dimension structure 0044 % .Ndipole : # of vertex 0045 % .NdipoleL : # of vertex in Left cortex 0046 % .Npatch : # of patch 0047 % .Coord = 'SPM_Right_m'; 0048 % 0049 % [history] 0050 % Ver.2.0 New version for group analysis 0051 % 2015-12-27 M. Sato 0052 % 2017-01-30 Y. Takeda vb_find_nearest_point->vb_find_nearest_point_no_overlap 0053 % 2017-03-06 M. Sato create Vext, Fext. 0054 % 0055 % 0056 % Copyright (C) 2011, ATR All Rights Reserved. 0057 % License : New BSD License(see VBMEG_LICENSE.txt) 0058 0059 global vbmeg_inst; 0060 const = vb_define_verbose; 0061 0062 % Max distance for nearest point search 0063 Rmax = 30.0; % 30 mm 0064 0065 EPS = 1e-10; 0066 0067 Nvertex = brain_parm.Nvertex; 0068 0069 vb_disp('Load original subject brain'); 0070 0071 %--------------------------------------------------- 0072 % Load Subject original cortical surface (SPM_Right mm cordinate) 0073 %--------------------------------------------------- 0074 [V0L,F0L,n0L,V0R,F0R,n0R] = vb_load_orig_brain_surf(brain_parm); 0075 % V0L/R vertex coordinate (SPM_Right mm) 0076 % n0L/R unit vector 0077 % FOL/R triangle patch 0078 0079 % Extract cortex area index of original surface 0080 vb_disp('Extract cortex area index'); 0081 0082 [Vindex] = vb_extract_cortex(V0L,F0L,V0R,F0R,brain_parm); 0083 % Cortex area index of original FreeSurfer surface 0084 % Vindex.cortexL; 0085 % Vindex.cortexR + NdipoleL0; 0086 0087 %--------------------------------------------------- 0088 % Load standard brain cortical surface 0089 %--------------------------------------------------- 0090 [std_parm, std_brain_dir] = vb_set_icbm152(Nvertex); 0091 0092 std_brain = [std_brain_dir filesep std_parm.brain_file]; 0093 0094 [Vstd, Fmni] = vb_load_cortex(std_brain); 0095 % Fmni : triangle patch for MNI cordinate 0096 0097 load(std_brain,'Vinfo'); 0098 0099 Ndipole = Vinfo.Ndipole; 0100 NdipoleL = Vinfo.NdipoleL; 0101 NdipoleR = Ndipole - NdipoleL; 0102 0103 %--------------------------------------------------- 0104 % Map MNI-coordinate to individual's cortical surface 0105 %--------------------------------------------------- 0106 0107 switch reg_mode 0108 case 'SPM' 0109 vb_disp('SPM registration Mode'); 0110 0111 % reference coordinate: subject original SPM-R mm coordinate 0112 VLref = V0L; 0113 VRref = V0R; 0114 0115 snfile = brain_parm.spm_normalization_file; 0116 0117 % target coordinate: MNI cordinate in standard brain 0118 Vmni = vb_load_cortex(std_brain,'MNI'); 0119 0120 % Back transformation to subject brain 0121 % Use mm-coordinates 0122 sn = load((snfile)); 0123 [Vsubj] = vb_unsn(Vmni'*1000,sn); 0124 0125 % target coordinate: MNI cordinate 0126 % transformed to SPM-R mm-coordinates in subject brain 0127 VL = Vsubj(:,1:NdipoleL)'; 0128 VR = Vsubj(:,(NdipoleL+1):Ndipole)'; 0129 0130 case 'FS' 0131 vb_disp('FS spherical coordinate registration Mode'); 0132 0133 % reference coordinate: subject original spherical coordinate in FreeSurfer 0134 VLref = vb_fs_load_surface(brain_parm.FS_left_sphere_file ); 0135 VRref = vb_fs_load_surface(brain_parm.FS_right_sphere_file); 0136 0137 % target coordinate: spherical coordinate of MNI standard brain 0138 Vsph = vb_load_cortex(std_brain,'sphere.reg'); 0139 0140 VL = Vsph(1:NdipoleL,:); 0141 VR = Vsph((NdipoleL+1):Ndipole,:); 0142 0143 otherwise 0144 error('Registration mode error'); 0145 end 0146 0147 vb_disp('--- Mapping MNI-coordinate onto a subject-brain'); 0148 %--------------------------------------------------- 0149 % Find nearest point from 'V' to 'VRref' 0150 %--------------------------------------------------- 0151 0152 % vertex number of reference original FreeSurfer surface 0153 NLref = size(VLref,1); 0154 NRref = size(VRref,1); 0155 0156 % Cortex area index of original FreeSurfer surface 0157 ref_corL = Vindex.cortexL; 0158 ref_corR = Vindex.cortexR; 0159 ref_nonL = vb_setdiff2([1:NLref], ref_corL); 0160 ref_nonR = vb_setdiff2([1:NRref], ref_corR); 0161 0162 % Cortex area index of MNI brain model 0163 corL = Vinfo.cortexL; 0164 corR = Vinfo.cortexR - NdipoleL; 0165 nonL = vb_setdiff2([1:NdipoleL], corL); 0166 nonR = vb_setdiff2([1:NdipoleR], corR); 0167 0168 % Cortex area vertex of original FreeSurfer surface 0169 VLref_cor = VLref(ref_corL,:); 0170 VRref_cor = VRref(ref_corR,:); 0171 VLref_non = VLref(ref_nonL,:); 0172 VRref_non = VRref(ref_nonR,:); 0173 0174 % Cortex area vertex of brain model 0175 VL_cor = VL(corL,:); 0176 VR_cor = VR(corR,:); 0177 VL_non = VL(nonL,:); 0178 VR_non = VR(nonR,:); 0179 0180 %--> Find nearest point from 'V' to 'Vref' 0181 [IL_cor, ddL_cor] = vb_find_nearest_point_no_overlap(VLref_cor, VL_cor,Rmax); 0182 [IR_cor, ddR_cor] = vb_find_nearest_point_no_overlap(VRref_cor, VR_cor,Rmax); 0183 [IL_non, ddL_non] = vb_find_nearest_point_no_overlap(VLref_non, VL_non,Rmax); 0184 [IR_non, ddR_non] = vb_find_nearest_point_no_overlap(VRref_non, VR_non,Rmax); 0185 0186 % index of brain model mapped from MNI brain 0187 IndxL = zeros(NdipoleL,1); 0188 IndxR = zeros(NdipoleR,1); 0189 0190 IndxL(corL) = ref_corL(IL_cor); 0191 IndxR(corR) = ref_corR(IR_cor); 0192 IndxL(nonL) = ref_nonL(IL_non); 0193 IndxR(nonR) = ref_nonR(IR_non); 0194 0195 ddL = zeros(NdipoleL,1); 0196 ddR = zeros(NdipoleR,1); 0197 0198 ddL(corL) = ddL_cor; 0199 ddR(corR) = ddR_cor; 0200 ddL(nonL) = ddL_non; 0201 ddR(nonR) = ddR_non; 0202 0203 dd = [ddL; ddR]; 0204 0205 % % Cortex area index of original FreeSurfer surface 0206 % NLref = size(VLref,1); 0207 % NRref = size(VRref,1); 0208 fprintf('# of subj_orig total vertex = %d (L), %d (R)\n', ... 0209 NLref,NRref) 0210 fprintf('# of subj_orig cortical vertex = %d (L), %d (R)\n',... 0211 length(ref_corL),length(ref_corR)) 0212 fprintf('# of target cortical vertex = %d (L), %d (R)\n',... 0213 length(corL),length(corR)) 0214 fprintf('# of target non-cortical vertex = %d (L), %d (R)\n\n',... 0215 length(nonL),length(nonR)) 0216 0217 fprintf('Max (Mean) distance in cortex = %4.1f (%4.1f) mm \n', ... 0218 max([ddL_cor(:); ddR_cor(:)]),mean([ddL_cor(:); ddR_cor(:)])) 0219 fprintf('Max (Mean) distance in non-cortex = %4.1f (%4.1f) mm \n', ... 0220 max([ddL_non(:); ddR_non(:)]),mean([ddL_non(:); ddR_non(:)])) 0221 0222 %--------------------------------------------------- 0223 % --- SPM cordinate in subject brain [m] 0224 %--------------------------------------------------- 0225 VL = V0L(IndxL,:); 0226 VR = V0R(IndxR,:); 0227 V = [VL; VR]/1000; 0228 0229 %--------------------------------------------------- 0230 % --- Normal vector obtained from original FS-model 0231 %--------------------------------------------------- 0232 n3L = n0L(IndxL,:); 0233 n3R = n0R(IndxR,:); 0234 xx = [n3L ; n3R]; 0235 [xx] = check_normal_vector(xx); 0236 0237 %--------------------------------------------------- 0238 % --- Get patch information for V 0239 %--------------------------------------------------- 0240 % 0241 % --- Reference coordinate: 0242 % Reduced Subject SPM-R mm cordinate 0243 % 0244 [FLs,VLs] = vb_reducepatch( F0L, V0L, Nvertex ); 0245 [FRs,VRs] = vb_reducepatch( F0R, V0R, Nvertex ); 0246 0247 % 0248 % ---- Check closed surface by solid_angle/(2*pi) = 1 0249 % 0250 [FLs, VLs, xxLs] = check_surface(FLs,VLs); 0251 [FRs, VRs, xxRs] = check_surface(FRs,VRs); 0252 0253 %%--------------------------------------------------- 0254 % --- Find nearest point of VL from Reduced subject cordinate VLs 0255 %--------------------------------------------------- 0256 0257 %--> Find nearest point from 'VL','VR' to 'VLs','VRs' 0258 [IL, ddL,VLs,FLs] = vb_find_no_overlap_divide(VLs, VL, FLs,Rmax); 0259 [IR, ddR,VRs,FRs] = vb_find_no_overlap_divide(VRs, VR, FRs,Rmax); 0260 0261 [Vext,FL,FR] = vb_calc_subj_patch(IL,IR,VLs,FLs,VRs,FRs); 0262 0263 % Add extra vertex to define complete reduced patch 0264 Vext = [V; Vext]; 0265 0266 % reduced patch for extended vertex (Vext) 0267 Fext.F3 = [FL ; FR]; 0268 Fext.F3R = FR; 0269 Fext.F3L = FL; 0270 0271 % reduced patch for model vertex (V) 0272 % select patch whose vertices consist of 'V' 0273 %FL = inner_patch_select([1:NdipoleL], FL); 0274 %FR = inner_patch_select([(1:NdipoleR) + NdipoleL], FR); 0275 %F.F3 = [FL ; FR]; 0276 %F.F3R = FR; 0277 %F.F3L = FL; 0278 % 0279 %F.NdipoleL = NdipoleL; 0280 %F.Fmni = Fmni; 0281 0282 F = []; 0283 0284 %--------------------------------------------------- 0285 % --- Vinfo : Dimensional info 0286 %--------------------------------------------------- 0287 0288 NdipoleL0 = size(VLref,1); 0289 Ndipole0 = size(VRref,1) + NdipoleL0; 0290 0291 %--------------------------------------------------- 0292 % Original vertex index corresponding to extended vertex (Vext) 0293 %--------------------------------------------------- 0294 %--> Find nearest point from 'Vref' for 'V' 0295 %[II, dd] = vb_find_nearest_point_no_overlap([V0L;V0R], Vext); 0296 %BV_index.Vall = II; 0297 0298 BV_index.Left = IndxL; 0299 BV_index.Right = IndxR + NdipoleL0; 0300 0301 % Cortex area index of original FreeSurfer surface 0302 BV_index.cortexL = Vindex.cortexL; 0303 BV_index.cortexR = Vindex.cortexR + NdipoleL0; 0304 0305 % Dimensional info 0306 % number of all vertex in 'VLs' 0307 NLall = size(VLs, 1); 0308 NRall = size(VRs, 1); 0309 0310 Vinfo.Ndipole = Ndipole; 0311 Vinfo.NdipoleL = NdipoleL; 0312 Vinfo.Ndipole_extL = Ndipole; 0313 Vinfo.Ndipole_extR = Ndipole + NLall; 0314 Vinfo.Ndipole0 = Ndipole0; 0315 Vinfo.NdipoleL0 = NdipoleL0; 0316 Vinfo.Coord = 'SPM_Right_m'; 0317 Vinfo.dd_match = dd; 0318 0319 return 0320 0321 function Fnew = inner_patch_select(Vix,F) 0322 % select patch whose vertices consist given vertex points 0323 0324 Nmax = max( max(Vix(:)) , max(F(:)) ); 0325 0326 Itrans = zeros(Nmax,1); 0327 Itrans(Vix) = 1:length(Vix); 0328 0329 Fnew = Itrans(F); 0330 0331 ix = find( Fnew(:,1).*Fnew(:,2).*Fnew(:,3) > 0); 0332 0333 Fnew = Fnew(ix,:); 0334 0335 return 0336 0337 function [FL, VL, xxL] = check_surface(FL,VL) 0338 0339 EPS = 1e-10; 0340 0341 [FL, VL, xxL] = vb_out_normal_surf(FL,VL); 0342 0343 % 0344 % ---- Check closed surface by solid_angle/(2*pi) = 1 0345 % 0346 omega = vb_solid_angle_check(VL,FL); 0347 vb_disp(sprintf('solid_angle/(2*pi) = %5.3f',omega)); 0348 if abs(omega - 1) > EPS, 0349 vb_disp('Surface is not closed '); 0350 end 0351 0352 return 0353 0354 function [xx] = check_normal_vector(xx) 0355 % Normalize normal vectors 0356 xxn = zeros(size(xx)); 0357 nn = sum(xx.^2 ,2); 0358 ix = find(nn > 0); 0359 xxn(ix,:) = xx(ix,:)./repmat(sqrt(nn(ix)) ,[1 3]); 0360 xx = xxn; 0361 0362 % Check zero normal vector 0363 if length(ix) == size(xx,1), 0364 vb_disp('--- All normal vector norms are 1'); 0365 else 0366 warning('There are zero normal vector'); 0367 end 0368 0369 return