Make standard brain from ICBM152
(VBMEG public function)
[syntax]
make_std_brain_model(proj_root,brain_parm);
[input]
proj_root : <<string>> VBMEG project root directory.
brain_parm: <<struct>> Parameters for creating cortical surface model
--- fields of brain_parm
analyze_file: <<string>> MRI filename (.hdr)
brain_parm.FS_left_file = [fs_dir '/bem/lh.smoothwm.asc'];
brain_parm.FS_right_file = [fs_dir '/bem/rh.smoothwm.asc'];
brain_parm.FS_left_infl_file = [fs_dir '/bem/lh.inflated.asc'];
brain_parm.FS_right_infl_file = [fs_dir '/bem/rh.inflated.asc'];
brain_parm.FS_left_curv_file = [fs_dir '/bem/lh.curv.asc'];
brain_parm.FS_right_curv_file = [fs_dir '/bem/rh.curv.asc'];
brain_parm.FS_left_label_file = [fs_dir '/label/lh.cortex.label'];
brain_parm.FS_right_label_file = [fs_dir '/label/rh.cortex.label'];
brain_parm.FS_left_sphere_file = [fs_dir '/bem/lh.sphere.reg.asc'];
brain_parm.FS_right_sphere_file = [fs_dir '/bem/rh.sphere.reg.asc'];
brain_parm.FS_sphere_key = 'sphere.reg';
brain_file : <<string>> Cortical surface model file (.brain.mat).
act_file : <<string>> Cortical activity map file (.act.mat).
area_file : <<string>> Cortical area file (.area.mat).
Nvertex : <<int>> Total number of vertices for cortical surface model
Rmax : <<float>> Maximum radius for neighbor search
N_step : <<int>> Division number of z axis in matching between
cortical surface model and MRI points
display : <<int>> Display progress for this steps
---
[output]
The following files are created:
VBMEG cortical surface model file (.brain.mat)
VBMEG area file (.area.mat)
VBMEG activity map file (.act.mat)
The following variables are stored in .brain.mat file:
V : Cortical vertex point cordinate (SPM_Right_m)
F : Patch index structure
xx : Normal vector to cortical surface
xxA : area asigned for each vertex ([m^2])
BV_index : original vertex index in BV/FS corresponding to brain
Vinfo : Vertex dimension structure
xxF{i} : Next-neighbor index for the vertex-i
xxD{i} : Distance between next-neighbor and the vertex-i
xxT{n} : vertex index for n-th triangle
nextIX{i} : Neighbor index list for the vertex-i
nextDD{i} : Distance from the vertex-i
The following areas are automatically created and stored in .area.mat
file:
areakey = 'Cortex' : whole brain
areakey = 'cortex_region' : cortex region without corpus area
areakey = 'corpus_area' : corpus area
[history]
Ver.2.0 New version for group analysis
2015-12-27 M. Sato
Copyright (C) 2011, ATR All Rights Reserved.
License : New BSD License(see VBMEG_LICENSE.txt)0001 function [V,F,xx,BV_index,Vinfo] = ... 0002 make_std_brain_model(proj_root,brain_parm, DEBUG_MODE) 0003 % Make standard brain from ICBM152 0004 % (VBMEG public function) 0005 % 0006 % [syntax] 0007 % make_std_brain_model(proj_root,brain_parm); 0008 % 0009 % [input] 0010 % proj_root : <<string>> VBMEG project root directory. 0011 % brain_parm: <<struct>> Parameters for creating cortical surface model 0012 % --- fields of brain_parm 0013 % analyze_file: <<string>> MRI filename (.hdr) 0014 % 0015 % brain_parm.FS_left_file = [fs_dir '/bem/lh.smoothwm.asc']; 0016 % brain_parm.FS_right_file = [fs_dir '/bem/rh.smoothwm.asc']; 0017 % brain_parm.FS_left_infl_file = [fs_dir '/bem/lh.inflated.asc']; 0018 % brain_parm.FS_right_infl_file = [fs_dir '/bem/rh.inflated.asc']; 0019 % brain_parm.FS_left_curv_file = [fs_dir '/bem/lh.curv.asc']; 0020 % brain_parm.FS_right_curv_file = [fs_dir '/bem/rh.curv.asc']; 0021 % brain_parm.FS_left_label_file = [fs_dir '/label/lh.cortex.label']; 0022 % brain_parm.FS_right_label_file = [fs_dir '/label/rh.cortex.label']; 0023 % 0024 % brain_parm.FS_left_sphere_file = [fs_dir '/bem/lh.sphere.reg.asc']; 0025 % brain_parm.FS_right_sphere_file = [fs_dir '/bem/rh.sphere.reg.asc']; 0026 % 0027 % brain_parm.FS_sphere_key = 'sphere.reg'; 0028 % 0029 % brain_file : <<string>> Cortical surface model file (.brain.mat). 0030 % act_file : <<string>> Cortical activity map file (.act.mat). 0031 % area_file : <<string>> Cortical area file (.area.mat). 0032 % Nvertex : <<int>> Total number of vertices for cortical surface model 0033 % Rmax : <<float>> Maximum radius for neighbor search 0034 % N_step : <<int>> Division number of z axis in matching between 0035 % cortical surface model and MRI points 0036 % display : <<int>> Display progress for this steps 0037 % 0038 % --- 0039 % 0040 % [output] 0041 % The following files are created: 0042 % VBMEG cortical surface model file (.brain.mat) 0043 % VBMEG area file (.area.mat) 0044 % VBMEG activity map file (.act.mat) 0045 % 0046 % The following variables are stored in .brain.mat file: 0047 % V : Cortical vertex point cordinate (SPM_Right_m) 0048 % F : Patch index structure 0049 % xx : Normal vector to cortical surface 0050 % xxA : area asigned for each vertex ([m^2]) 0051 % 0052 % BV_index : original vertex index in BV/FS corresponding to brain 0053 % Vinfo : Vertex dimension structure 0054 % 0055 % xxF{i} : Next-neighbor index for the vertex-i 0056 % xxD{i} : Distance between next-neighbor and the vertex-i 0057 % xxT{n} : vertex index for n-th triangle 0058 % nextIX{i} : Neighbor index list for the vertex-i 0059 % nextDD{i} : Distance from the vertex-i 0060 % 0061 % The following areas are automatically created and stored in .area.mat 0062 % file: 0063 % areakey = 'Cortex' : whole brain 0064 % areakey = 'cortex_region' : cortex region without corpus area 0065 % areakey = 'corpus_area' : corpus area 0066 % 0067 % [history] 0068 % Ver.2.0 New version for group analysis 0069 % 2015-12-27 M. Sato 0070 % 0071 % Copyright (C) 2011, ATR All Rights Reserved. 0072 % License : New BSD License(see VBMEG_LICENSE.txt) 0073 0074 global vbmeg_inst 0075 const = vb_define_verbose; 0076 VERBOSE_LEVEL_NOTICE = const.VERBOSE_LEVEL_NOTICE; 0077 0078 % Set BV-SPM version (mode = 4 : std_brain mode) 0079 if ~exist('DEBUG_MODE','var'), DEBUG_MODE = 0; end; 0080 0081 % Define constant 0082 define = vbmeg_inst.const; 0083 0084 proj_root = vb_rm_trailing_slash(proj_root); 0085 if exist(proj_root, 'dir') ~= 7 0086 vb_mkdir(proj_root); 0087 end 0088 0089 % 0090 % Do not modify following lines 0091 % 0092 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0093 % 0094 0095 % Output file path 0096 if isempty(proj_root) 0097 brain_file = [brain_parm.brain_file]; 0098 area_file = [brain_parm.area_file ]; 0099 act_file = [brain_parm.act_file ]; 0100 else 0101 brain_file = [proj_root filesep brain_parm.brain_file]; 0102 area_file = [proj_root filesep brain_parm.area_file ]; 0103 act_file = [proj_root filesep brain_parm.act_file ]; 0104 end 0105 0106 % ID for MRI image data 0107 analyze_file = brain_parm.analyze_file; 0108 0109 [udir, fname] = fileparts(analyze_file); 0110 MRI_ID = fname; 0111 0112 flag_fs = (isfield(brain_parm, 'FS_left_sphere_file') ... 0113 && ~isempty(brain_parm.FS_left_sphere_file) ... 0114 && isfield(brain_parm, 'FS_right_sphere_file') ... 0115 && ~isempty(brain_parm.FS_right_sphere_file) ... 0116 && isfield(brain_parm, 'FS_sphere_key') ... 0117 && ~isempty(brain_parm.FS_sphere_key)); 0118 0119 if flag_fs==0, 0120 error('FreeSurfer file setting error'); 0121 end 0122 0123 Rmax = brain_parm.R_max; 0124 Display = brain_parm.display; 0125 0126 %------------------------------------------------------- 0127 % Make cortex coordinate from FreeSurfer 0128 %------------------------------------------------------- 0129 [V,F,xx,BV_index,Vinfo] = vb_make_brain_std(brain_parm); 0130 0131 if DEBUG_MODE==2, return; end; 0132 0133 %-------------------------------------------------------------------------- 0134 % Calculate MNI coordinate and Talairach coordinate 0135 %-------------------------------------------------------------------------- 0136 0137 Vmni = convert_spm_mni(V, analyze_file)*0.001; 0138 Vtal = mni2tal(Vmni); 0139 0140 % 0141 % Save cortical surface model 0142 % 0143 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0144 vb_disp(['Save brain model '], VERBOSE_LEVEL_NOTICE); 0145 vb_disp(sprintf(' filename = %s', brain_file), ... 0146 VERBOSE_LEVEL_NOTICE); 0147 vb_save([brain_file],... 0148 'F','V','xx','Vinfo','BV_index','MRI_ID','Vmni','Vtal'); 0149 0150 %---------------------------------------- 0151 % Search next-point index and distance 0152 %---------------------------------------- 0153 tic; 0154 vb_disp_nonl('Search next-point index and distance ', ... 0155 VERBOSE_LEVEL_NOTICE); 0156 0157 [xxD, xxF, xxT] = vb_next_distance( F.F3, V ); % unit : m 0158 0159 vb_disp(sprintf('%f[sec]',toc), VERBOSE_LEVEL_NOTICE); 0160 0161 % total number of vertex 0162 Ndipole = Vinfo.Ndipole; 0163 Vindx = 1:Ndipole; 0164 0165 %---------------------------------------------- 0166 % Calculate area assigned to the vertex 0167 %---------------------------------------------- 0168 tic; 0169 [xxA] = vb_calc_patch_area(V, F.F3, xxT); 0170 vb_disp(sprintf('%f[sec]',toc), VERBOSE_LEVEL_NOTICE); 0171 0172 % 0173 % Save next neighbor index 0174 % 0175 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0176 vb_disp(['Save next neighbor index'], VERBOSE_LEVEL_NOTICE); 0177 vb_save([brain_file], 'xxF','xxD','xxA'); 0178 0179 %--------------------------------------------------- 0180 % Find original vertex near the cortical vertex 0181 %--------------------------------------------------- 0182 normal_stat = vb_original_normal_statics(proj_root,brain_parm,BV_index); 0183 0184 % normal average over neighbor points of original brain 0185 xx = vb_calc_normal_average(normal_stat.neighbor_org, normal_stat.normal_org); 0186 0187 vb_disp('Save neighbor statics in the original brain', ... 0188 VERBOSE_LEVEL_NOTICE); 0189 vb_save(brain_file,'normal_stat','xx'); 0190 0191 %-------------------------------------------------------------------------- 0192 % Save FreeSurfer sphere coordinate 0193 %-------------------------------------------------------------------------- 0194 parm.brain_file = brain_parm.brain_file; 0195 parm.brain_sphR = brain_parm.FS_right_sphere_file; 0196 parm.brain_sphL = brain_parm.FS_left_sphere_file; 0197 parm.key = brain_parm.FS_sphere_key; 0198 vb_job_brain_add_sphcoord(proj_root, parm); 0199 0200 %-------------------------------------------------------------------------- 0201 % Make two area files derived from two anatomical divisons 0202 %-------------------------------------------------------------------------- 0203 vb_atlas_label(brain_file, 'aal'); 0204 vb_remove_island_atlas(brain_file, 'aal'); 0205 0206 %vb_atlas_label(brain_file, 'brodmann'); 0207 %vb_remove_island_atlas(brain_file, 'brodmann'); 0208 0209 if DEBUG_MODE==1, return; end; 0210 0211 %---------------------------------------------- 0212 % Search neighbor points along cortex sheet 0213 %---------------------------------------------- 0214 0215 [nextIX , nextDD] = vb_find_neighbor_all(Rmax, xxF, xxD, Vindx, Display); 0216 0217 vb_disp('Save neighbor index of cortical vertex', ... 0218 VERBOSE_LEVEL_NOTICE); 0219 vb_save([brain_file],'nextIX','nextDD'); 0220 0221 %-------------------------------------- 0222 % Make new area and act files 0223 %-------------------------------------- 0224 act_new.key = 'Uniform'; 0225 act_new.xxP = ones(Ndipole,1); 0226 act_new.comment = 'homogeneous fMRI constraint'; 0227 0228 vb_add_act([act_file], act_new, MRI_ID, OFF); 0229 0230 % Make Area file 0231 Vindx = [Vinfo.cortexL; Vinfo.cortexR]; 0232 AreaNew.key = 'Cortex'; 0233 AreaNew.Iextract = Vindx; %[1:Ndipole]'; 0234 0235 vb_add_area([area_file], AreaNew, MRI_ID, OFF); 0236 0237 %----------------------------------------- 0238 % Check brain model by 3D & MRI image 0239 %----------------------------------------- 0240 vb_check_brain_model(proj_root,brain_parm); 0241 0242 0243 %-------------------------- 0244 % Import inflate model 0245 %-------------------------- 0246 if has_inflate_model(brain_parm) 0247 vb_job_inflate(proj_root, brain_parm); 0248 end 0249 0250 %-------------------------- 0251 % Save project file 0252 %-------------------------- 0253 proj_file = get_project_filename; 0254 if ~isempty(proj_file) 0255 project_file_mgr('load', proj_file); 0256 project_file_mgr('add', 'brain_parm', brain_parm); 0257 end 0258 0259 return; 0260 0261 function [result] = has_inflate_model(brain_parm) 0262 % Check if brain_parm contains valid fields for inflate model. 0263 % result : return bit pattern 0264 % = 01(1 as int): BV inflate data is set. 0265 % = 10(2 as int): FS inflate data is set. 0266 % = 11(3 as int): Both inflate model were set 0267 result = 0; 0268 BV_BIT = 1; FS_BIT = 2; 0269 % bit all on 0270 result = bitset(result, BV_BIT, 1); 0271 result = bitset(result, FS_BIT, 1); 0272 0273 BV_fields = {'BV_left_infl_file', 'BV_right_infl_file'}; 0274 FS_fields = {'FS_left_infl_file', 'FS_right_infl_file', ... 0275 'FS_left_curv_file', 'FS_right_curv_file'}; 0276 for k=1:length(BV_fields) 0277 if ~isfield(brain_parm, BV_fields{k}) || ... 0278 exist(brain_parm.(BV_fields{k}), 'file') ~= 2 0279 % bit off 0280 result = bitset(result, BV_BIT, 0); 0281 break; 0282 end 0283 end 0284 for k=1:length(FS_fields) 0285 if ~isfield(brain_parm, FS_fields{k}) || ... 0286 exist(brain_parm.(FS_fields{k}), 'file') ~= 2 0287 % bit off 0288 result = bitset(result, FS_BIT, 0); 0289 break; 0290 end 0291 end 0292 0293 function [Vmni, origin] = convert_spm_mni(Vspm, fname) 0294 % Transform SPM-Right-m coordinate to MNI mm coordinate 0295 % [Vmni, origin] = convert_spm_mni(Vspm, fname) 0296 % Vmni ; MNI [mm] coordinate (Npoint x 3) 0297 % fname : file name (MNI_152) of T1 image file 0298 % Vspm : SPM-Right [m] coordinate 0299 % origin: origin of MNI cordinate in SPM-Right [m] coordinate 0300 % Vspm = Vmni + origin; (Vmni = 0 at origin) 0301 0302 [Trans , dim] = get_coord_trans_mat(fname); 0303 0304 % Trans : transform matrix from voxcel to MNI mm coordinate 0305 % [x y z 1] = [i j k 1] * Trans 0306 % Xmni = Vox * Trans 0307 % Vox = Xmni * inv(Trans) 0308 0309 % Center position of 3D image in MNI mm 0310 center = vb_affine_trans(dim/2 , Trans); 0311 0312 % Xspm = 0 <-> Xmni = center 0313 % Xspm = -center <-> Xmni = 0 0314 % Xspm = Xmni - center 0315 % Xmni = Xspm + center 0316 0317 % Get MNI mm coordinate 0318 Vmni = vb_repadd(Vspm * 1000 , center) ; 0319 0320 % origin of MNI cordinate in SPM-Right coordinate 0321 %origin = -center ; 0322 origin = -center * 0.001;