mesh_3shell_script Assumes preprocessing of MRI with FreeSurfer http://surfer.nmr.mgh.harvard.edu/ FSL tools (BET & FAST) http://fmrib.ox.ac.uk/ An example FSL script can be found at the end of this script
0001 0002 % mesh_3shell_script 0003 % 0004 % Assumes preprocessing of MRI with 0005 % FreeSurfer http://surfer.nmr.mgh.harvard.edu/ 0006 % FSL tools (BET & FAST) http://fmrib.ox.ac.uk/ 0007 % 0008 % An example FSL script can be found at the end of this script 0009 % 0010 0011 % Licence: GNU GPL, no implied or express warranties 0012 % History: 08/2002, Darren.Weber_at_radiology.ucsf.edu 0013 % 08/2003, Darren.Weber_at_radiology.ucsf.edu 0014 % added more control variables and the FSL script 0015 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0016 0017 0018 % the plot option is exclusive of the following commands, 0019 % it should be done after initial attempt to get scalp etc. 0020 plot = 0; 0021 0022 % the initial results may not be very pretty, these tools 0023 % are not ideal! It usually takes some time to tweak the 0024 % shrinkwrap results and then finally use bem_shells_correct to 0025 % get reasonable BEM surfaces. All development work needs to 0026 % be done on modifying the shrinkwrap method (it's not great!). 0027 doScalp = 0; 0028 doIskull = 0; 0029 doOskull = 0; 0030 doCortex = 0; 0031 0032 % only use this after initial processing above, use the plot 0033 % option to check the initial results. This option will use 0034 % bem_shells_correct, which can do a neat job of correcting 0035 % the surfaces for intersections. 0036 doCorrection = 0; 0037 0038 doOutputEMSE = 0; 0039 doOutputBrainStorm = 0; 0040 0041 %IMG.path = 'D:\freesurfer\subjects\ucsf_morgan\mri\analyze'; 0042 IMG.path = '/data/ucsf/mri/ucsf_morgan/mri/analyze'; 0043 IMG.meshpath = '/data/ucsf/mri/ucsf_morgan/surf'; 0044 0045 IMG.scalp = 'morgan_orig_axial_las_smooth'; 0046 IMG.oskull = 'morgan_orig_axial_las_skull'; 0047 IMG.iskull = 'morgan_orig_axial_las_brain'; 0048 0049 IMG.intensity.scalp = 60; 0050 IMG.tolerance.scalp = 20; 0051 IMG.intensity.oskull = 15; 0052 IMG.tolerance.oskull = 10; 0053 IMG.intensity.iskull = 80; 0054 IMG.tolerance.iskull = 40; 0055 0056 if ~exist('p','var'), 0057 [p] = eeg_toolbox_defaults; 0058 end 0059 0060 if plot, 0061 0062 % plot scalp surface (index 4), in red 0063 if ~isempty(p.mesh.data.vertices{4}), 0064 patch('vertices',p.mesh.data.vertices{4},'faces',p.mesh.data.faces{4},... 0065 'FaceColor',[.6 0 0],'Edgecolor','none','FaceAlpha',.2); 0066 end 0067 % plot outer skull surface (index 3), in green 0068 if ~isempty(p.mesh.data.vertices{3}), 0069 patch('vertices',p.mesh.data.vertices{3},'faces',p.mesh.data.faces{3},... 0070 'FaceColor',[.0 .6 .0],'Edgecolor',[.8 .8 .8],'FaceAlpha',.4); 0071 end 0072 % plot inner skull surface (index 2), in blue 0073 if ~isempty(p.mesh.data.vertices{2}), 0074 patch('vertices',p.mesh.data.vertices{2},'faces',p.mesh.data.faces{2},... 0075 'FaceColor',[.0 .0 .6],'Edgecolor',[.8 .8 .8],'FaceAlpha',.6); 0076 end 0077 % plot cortex surface (index 1), in black 0078 if ~isempty(p.mesh.data.vertices{1}), 0079 patch('vertices',p.mesh.data.vertices{1},'faces',p.mesh.data.faces{1},... 0080 'FaceColor',[.0 .0 .0],'Edgecolor',[.8 .8 .8],'FaceAlpha',1); 0081 end 0082 0083 camlight; mouse_rotate; 0084 axis tight; hold on 0085 0086 % FV = ISkull; 0087 % patch('vertices',FV.vertices,'faces',FV.faces,'FaceColor',[.6 .0 .0],'Edgecolor','none','FaceAlpha',.2); axis tight; hold on 0088 % x = FV.vertices(1,1); 0089 % y = FV.vertices(1,2); 0090 % z = FV.vertices(1,3); 0091 % plot3(x,y,z,'ro') 0092 % 0093 return 0094 end 0095 0096 0097 0098 % TESSELATE SKULL/SCALP SURFACES 0099 0100 0101 % Image files to process 0102 0103 % These images were created with the freesurfer mri_convert command, eg: 0104 % mri_convert -oid -1 0 0 -ojd 0 1 0 -okd 0 0 1 orig mri\analyze\morgan_orig_axial_las.img 0105 % These Analyze files were then processed with the FSL tools to find the 0106 % skull. 0107 0108 cd(IMG.path) 0109 0110 if doScalp, 0111 0112 % --- Tesselate Scalp --- 0113 fprintf('\n\n**************************************\nProcessing Scalp\n\n'); 0114 Nsurf = 4; 0115 intensity = IMG.intensity.scalp; 0116 tolerance = IMG.tolerance.scalp; 0117 avw = avw_img_read(IMG.scalp); 0118 % [FV, Edges] = avw_shrinkwrap(avw,FV,smooth,vthresh,interpVal,... 0119 % fitval,fittol,fititer,fitchange,fitvattr) 0120 scalp = avw_shrinkwrap(avw,[],0,0,[],intensity,tolerance,50,0.5,0.4); 0121 if isfield(scalp,'vertices'), 0122 p.mesh.data.meshtype{Nsurf} = 'scalp'; 0123 p.mesh.data.vertices{Nsurf} = scalp.vertices ./ 1000; % convert to meters 0124 p.mesh.data.faces{Nsurf} = scalp.faces; 0125 end 0126 end 0127 0128 0129 if doIskull, 0130 0131 % --- Tesselate Inner skull 0132 fprintf('\n\n**************************************\nProcessing "Inner Skull"\n\n'); 0133 Nsurf = 2; 0134 intensity = IMG.intensity.iskull; 0135 tolerance = IMG.tolerance.iskull; 0136 avw = avw_img_read(IMG.iskull); 0137 % [FV, Edges] = avw_shrinkwrap(avw,FV,smooth,vthresh,interpVal,... 0138 % fitval,fittol,fititer,fitchange,fitvattr) 0139 ISkull = avw_shrinkwrap(avw,[],0,0,[],intensity,tolerance,50,0.5,0.4); 0140 if isfield(ISkull,'vertices'), 0141 p.mesh.data.meshtype{Nsurf} = 'inner_skull'; 0142 p.mesh.data.vertices{Nsurf} = ISkull.vertices ./ 1000; % convert to meters 0143 p.mesh.data.faces{Nsurf} = ISkull.faces; 0144 end 0145 end 0146 0147 0148 if doOskull, 0149 0150 % --- Tesselate Outer skull 0151 fprintf('\n\n**************************************\nProcessing "Outer Skull"\n\n'); 0152 Nsurf = 3; 0153 intensity = IMG.intensity.oskull; 0154 tolerance = IMG.tolerance.oskull; 0155 avw = avw_img_read(IMG.oskull); 0156 % [FV, Edges] = avw_shrinkwrap(avw,FV,smooth,vthresh,interpVal,... 0157 % fitval,fittol,fititer,fitchange,fitvattr) 0158 OSkull = avw_shrinkwrap(avw,[],0,0,0,intensity,tolerance,50,0.05,0.1); 0159 if isfield(OSkull,'vertices'), 0160 p.mesh.data.meshtype{Nsurf} = 'outer_skull'; 0161 p.mesh.data.vertices{Nsurf} = OSkull.vertices ./ 1000; % convert to meters 0162 p.mesh.data.faces{Nsurf} = OSkull.faces; 0163 %patch('vertices',OSkull.vertices,'faces',OSkull.faces,... 0164 % 'FaceColor',[.6 .0 .0],'Edgecolor',[.8 .8 .8],... 0165 % 'FaceAlpha',.2); axis tight; hold on 0166 end 0167 end 0168 0169 if doCortex, 0170 0171 % LOAD FREESURFER CORTEX 0172 % Freesurfer surface to process 0173 % This freesurfer surface contains the whole brain surface, after 0174 % altering the freesurfer pons/corpus callosum cutting positions. 0175 p.mesh.path = IMG.meshpath; 0176 p.mesh.file = 'rh.pial'; 0177 p.mesh.type = 'fs_surf'; 0178 0179 [p] = mesh_open(p); 0180 0181 end 0182 0183 if doCorrection, 0184 0185 dist.four2three = 6; % scalp to outer skull 0186 dist.four2two = 8; % scalp to inner skull 0187 dist.three2two = 5; % outer skull to inner skull 0188 dist.two2one = 2; % inner skull to cortex 0189 0190 cortex = 1; % process cortex (1 = yes, 0 = no) 0191 [p] = mesh_bem_correct(p, [], dist, 1, cortex); 0192 0193 end 0194 0195 0196 if doOutputEMSE, 0197 mesh_write_emse(p); 0198 end 0199 if doOutputBrainStorm, 0200 mesh_write_brainstorm(p); 0201 end 0202 0203 0204 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0205 % Example FSL script 0206 % Copy this text into a text file and chmod a+x that file. 0207 % Remove the % matlab comments and everything before the next line, 0208 % including these instructions. 0209 % #! /bin/bash 0210 % 0211 % 0212 % imagepath="." 0213 % imagetype="_orig_axial_las" 0214 % 0215 % sub="subject"; 0216 % img=`echo $sub$imagetype`; 0217 % 0218 % printf "\nprocessing $img\n"; 0219 % 0220 % printf "\n5mm Guassian spatial smoothing: ${img}_smooth.*\n"; 0221 % /usr/local/fsl/bin/ip ${img} ${img}_smooth 0 -s 5 0222 % 0223 % printf "\nUsing BET for brain extraction: ${img}_brain.*\n"; 0224 % bet ${img} ${img}_brain -s 0225 % 0226 % # Binarise the bet skull 0227 % #avwmaths ${img}_brain_skull -bin BETskull 0228 % # Dilate bet skull and add a constant 0229 % #avwmaths_32R BETskull -dil -add 10000 BETskull 0230 % 0231 % printf "\nSkull extraction and 2mm Guassian spatial smoothing: ${img}_skull.*\n"; 0232 % /usr/local/fsl/bin/ip ${img}_brain_skull ${img}_skull 0 -s 2 0233 % 0234 % printf "\nSpatial normalization\n"; 0235 % flirt -in ${img}_brain.hdr -ref /usr/local/fsl/etc/standard/avg152T1_brain.hdr -out ${img}_brain_reg.hdr -omat %{img}_brain_reg.mat -bins 256 -cost corratio -searchrx -180 180 -searchry -180 180 -searchrz -180 180 -dof 12 -interp trilinear 0236 % convert_xfm -matonly -omat ${img}_brain_reginv.mat -inverse ${img}_brain_reg.mat 0237 % 0238 % printf "\nBrain segmentation\n"; 0239 % fast -t 1 -c 3 -op -e -ov -od $img ${img}_brain 0240 % 0241 % exit