Home > functions > job > vb_job_head_3shell.m

vb_job_head_3shell

PURPOSE ^

Make 1 shell or 3 shell surface model for MEG/EEG

SYNOPSIS ^

function head_join = vb_job_head_3shell(proj_root, parm)

DESCRIPTION ^

 Make 1 shell or 3 shell surface model for MEG/EEG
     CSF (CSF + Gray + White matter) : 1 shell model surface (innermost)
     Skull (Bone)
     Scalp (Head Skin) (outermost)

   - CSF is imported from FreeSurfer/Curry surface 
     or made by using brain model and segmented gray matter file by SPM
   - Scalp surface is made from face and dilated CSF
   - Outer Skull surface is made by erosion from Scalp surface
 --- Usage
   head_join = vb_job_head_3shell(proj_root, parm)
   proj_root : Output file root directory
   parm      : Parameters explained below
   head_join : output head file name 
    (3 surfaces are joined in one file for 3 shell model)

 --- Setting for Output file  [relative path from 'proj_root']
 parm.head_file : Base name for head model

 --- Setting for Output surface
 parm.Nsurf   : Number of shell (=1 for MEG , =3 for EEG [default])
 parm.Nvertex : number of vertex in each surface 
              = [2000] (Default)
  if length(Nvertex) = 1, Nvertex for all surface are set to the same value

 --- Setting for MRI input file  (Absolute path)
   parm.analyze_file : MRI image file

 --- Setting for Scalp input file (Absolute path) [There are 2 cases]
    If the following files are specified at the same time,
    priority order is fs_scalp_file > face_file
 1.  parm.face_file     : face file made by vbmeg
 2.  parm.fs_scalp_file : FreeSurfer Scalp surface ('outer_skin_surface.asc')

 --- Setting for CSF input files  [There are 3 cases]
     (Absolute path except for 'brain_file')
   We recommend to use both FreeSurfer and Segment files:
      'freesurf_file' + 'gray_file' + 'brain_file'

 1. FreeSurfer file Case
   parm.freesurf_file : FreeSurfer CSF surface ('inner_skull_surface.asc')
 2. Segment file Case
   parm.gray_file  : Gray matter segment file made by SPM 'Segment'
   parm.brain_file : Brain model file (.brain.mat)
        ['brain_file' should be relative path from 'proj_root']

 3. Curry file Case
   parm.curry_file : Curry surface file

%% --- Advanced Optional Parameter
   !!! Default parameters are defined by 'vb_set_head_parm'

 parm.Scalp_mode 
     = 0 : Scalp surface is intersection of face and dilated CSF [Default]
     = 1 : Scalp surface is smoothed face & Skull is the same as mode=0
     = 2 : Scalp surface is smoothed face & Skull is erosion of face

 --- Optional setting for outer Skull surface   (Absolute path)
 parm.fs_skull_file : Outer Skull surface ('/outer_skull_surface.asc')

 --- Radius of Morphology operation [mm]
 If the obtained shell model is not good,
 you can adjust the following morphology radius parameters

 parm.Radius_scalp = 4 : Thickness of Scalp [mm]
 parm.Radius_skull = 2 : Min Thickness of Skull [mm]
 parm.Radius_fs  : smoothing radius for FreeSurfer surface
                 = [4 -4 ] Default (No size change)
 parm.Radius_csf : smoothing radius for Gray matter
                 = [6 6 -7 ] Default (Expand 6 mm)
 parm.Radius     : smoothing radius for CSF surface
                 = [4 -4 ] Default (No size change)
 
 parm.Nloop  = 200; Iteration number for fit boundary
 parm.Nlast_csf   = 10 ; Iteration number for CSF smoothing by spring model
                    As Nlast becomes larger, surface becomes smoother
                    but it shrinks inward
 parm.Nlast_scalp  = 0 ; Iteration number for Scalp
 parm.Nlast_skull  = 0 ; Iteration number for Skull
 (Gray matter parameters)
 parm.Glevel      = [0.8];  level of gray matter selection threshold
 parm.Radius_gray = [-2 2]; Radius to remove noise in gray matter image
 (Brain model parameters)
 parm.Radius_fill = [6];    Cortex fill radius

 Masa-aki Sato 2009-10-7
 Masa-aki Sato 2010-09-15  add FreeSurfer case
 Masa-aki Sato 2010-10-27  modified procedures

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function    head_join = vb_job_head_3shell(proj_root, parm)
0002 % Make 1 shell or 3 shell surface model for MEG/EEG
0003 %     CSF (CSF + Gray + White matter) : 1 shell model surface (innermost)
0004 %     Skull (Bone)
0005 %     Scalp (Head Skin) (outermost)
0006 %
0007 %   - CSF is imported from FreeSurfer/Curry surface
0008 %     or made by using brain model and segmented gray matter file by SPM
0009 %   - Scalp surface is made from face and dilated CSF
0010 %   - Outer Skull surface is made by erosion from Scalp surface
0011 % --- Usage
0012 %   head_join = vb_job_head_3shell(proj_root, parm)
0013 %   proj_root : Output file root directory
0014 %   parm      : Parameters explained below
0015 %   head_join : output head file name
0016 %    (3 surfaces are joined in one file for 3 shell model)
0017 %
0018 % --- Setting for Output file  [relative path from 'proj_root']
0019 % parm.head_file : Base name for head model
0020 %
0021 % --- Setting for Output surface
0022 % parm.Nsurf   : Number of shell (=1 for MEG , =3 for EEG [default])
0023 % parm.Nvertex : number of vertex in each surface
0024 %              = [2000] (Default)
0025 %  if length(Nvertex) = 1, Nvertex for all surface are set to the same value
0026 %
0027 % --- Setting for MRI input file  (Absolute path)
0028 %   parm.analyze_file : MRI image file
0029 %
0030 % --- Setting for Scalp input file (Absolute path) [There are 2 cases]
0031 %    If the following files are specified at the same time,
0032 %    priority order is fs_scalp_file > face_file
0033 % 1.  parm.face_file     : face file made by vbmeg
0034 % 2.  parm.fs_scalp_file : FreeSurfer Scalp surface ('outer_skin_surface.asc')
0035 %
0036 % --- Setting for CSF input files  [There are 3 cases]
0037 %     (Absolute path except for 'brain_file')
0038 %   We recommend to use both FreeSurfer and Segment files:
0039 %      'freesurf_file' + 'gray_file' + 'brain_file'
0040 %
0041 % 1. FreeSurfer file Case
0042 %   parm.freesurf_file : FreeSurfer CSF surface ('inner_skull_surface.asc')
0043 % 2. Segment file Case
0044 %   parm.gray_file  : Gray matter segment file made by SPM 'Segment'
0045 %   parm.brain_file : Brain model file (.brain.mat)
0046 %        ['brain_file' should be relative path from 'proj_root']
0047 %
0048 % 3. Curry file Case
0049 %   parm.curry_file : Curry surface file
0050 %
0051 %%% --- Advanced Optional Parameter
0052 %   !!! Default parameters are defined by 'vb_set_head_parm'
0053 %
0054 % parm.Scalp_mode
0055 %     = 0 : Scalp surface is intersection of face and dilated CSF [Default]
0056 %     = 1 : Scalp surface is smoothed face & Skull is the same as mode=0
0057 %     = 2 : Scalp surface is smoothed face & Skull is erosion of face
0058 %
0059 % --- Optional setting for outer Skull surface   (Absolute path)
0060 % parm.fs_skull_file : Outer Skull surface ('/outer_skull_surface.asc')
0061 %
0062 % --- Radius of Morphology operation [mm]
0063 % If the obtained shell model is not good,
0064 % you can adjust the following morphology radius parameters
0065 %
0066 % parm.Radius_scalp = 4 : Thickness of Scalp [mm]
0067 % parm.Radius_skull = 2 : Min Thickness of Skull [mm]
0068 % parm.Radius_fs  : smoothing radius for FreeSurfer surface
0069 %                 = [4 -4 ] Default (No size change)
0070 % parm.Radius_csf : smoothing radius for Gray matter
0071 %                 = [6 6 -7 ] Default (Expand 6 mm)
0072 % parm.Radius     : smoothing radius for CSF surface
0073 %                 = [4 -4 ] Default (No size change)
0074 %
0075 % parm.Nloop  = 200; Iteration number for fit boundary
0076 % parm.Nlast_csf   = 10 ; Iteration number for CSF smoothing by spring model
0077 %                    As Nlast becomes larger, surface becomes smoother
0078 %                    but it shrinks inward
0079 % parm.Nlast_scalp  = 0 ; Iteration number for Scalp
0080 % parm.Nlast_skull  = 0 ; Iteration number for Skull
0081 % (Gray matter parameters)
0082 % parm.Glevel      = [0.8];  level of gray matter selection threshold
0083 % parm.Radius_gray = [-2 2]; Radius to remove noise in gray matter image
0084 % (Brain model parameters)
0085 % parm.Radius_fill = [6];    Cortex fill radius
0086 %
0087 % Masa-aki Sato 2009-10-7
0088 % Masa-aki Sato 2010-09-15  add FreeSurfer case
0089 % Masa-aki Sato 2010-10-27  modified procedures
0090 %
0091 % Copyright (C) 2011, ATR All Rights Reserved.
0092 % License : New BSD License(see VBMEG_LICENSE.txt)
0093 
0094 proj_root = vb_rm_trailing_slash(proj_root);
0095 
0096 % Default Parameters
0097 head_parm = vb_set_head_parm;
0098 
0099 % Number of surface
0100 if ~isfield(parm,'Nsurf'), 
0101     Nsurf = head_parm.Nsurf;
0102 else
0103     Nsurf = parm.Nsurf; 
0104 end;
0105 
0106 % Number of vertex in one surface [2000]
0107 if ~isfield(parm,'Nvertex'), 
0108     switch    Nsurf
0109     case    1
0110         Nvertex = [6000]; 
0111     case    3
0112         Nvertex = head_parm.Nvertex; 
0113     end
0114 else
0115     Nvertex = parm.Nvertex;
0116 end;
0117 if length(Nvertex)==1, Nvertex = repmat(Nvertex,[1,Nsurf]); end;
0118 
0119 % Set Default Radius parameters
0120 if ~isfield(parm,'Radius_fs'),   parm.Radius_fs    = head_parm.Radius_fs ;end;
0121 if ~isfield(parm,'Radius_csf'),  parm.Radius_csf   = head_parm.Radius_csf;end;
0122 if ~isfield(parm,'Radius_scalp'),parm.Radius_scalp = head_parm.Radius_scalp;end
0123 if ~isfield(parm,'Radius_skull'),parm.Radius_skull = head_parm.Radius_skull;end
0124 if ~isfield(parm,'Skull_ratio'), parm.Skull_ratio  = head_parm.Skull_ratio ;end
0125 
0126 if ~isfield(parm,'Radius_final'),parm.Radius_final = head_parm.Radius_final;end
0127 if ~isfield(parm,'Radius'),      parm.Radius       = head_parm.Radius     ;end;
0128 if ~isfield(parm,'Radius_fill'), parm.Radius_fill  = head_parm.Radius_fill;end;
0129 if ~isfield(parm,'Radius_gray'), parm.Radius_gray  = head_parm.Radius_gray;end;
0130 if ~isfield(parm,'Glevel'),      parm.Glevel       = head_parm.Glevel     ;end;
0131 
0132 if ~isfield(parm,'Nloop'),       parm.Nloop       = head_parm.Nloop      ;end;
0133 if ~isfield(parm,'Nlast_csf'),   parm.Nlast_csf   = head_parm.Nlast_csf  ;end;
0134 if ~isfield(parm,'Nlast_scalp'), parm.Nlast_scalp = head_parm.Nlast_scalp;end;
0135 if ~isfield(parm,'Nlast_skull'), parm.Nlast_skull = head_parm.Nlast_skull;end;
0136 
0137 if ~isfield(parm,'Scalp_mode'), parm.Scalp_mode = head_parm.Scalp_mode;end;
0138 if ~isfield(parm,'vstep'),         parm.vstep = head_parm.vstep; end;
0139 
0140 step = parm.vstep;
0141 
0142 % Plot mode
0143 if ~isfield(parm,'plot_mode'), 
0144     plot_mode = 0;
0145 else
0146     plot_mode = parm.plot_mode; 
0147 end;
0148 
0149 % Brain & image file name
0150 analyzefile = parm.analyze_file;
0151 if isfield(parm, 'brain_file')
0152     brain_file = [proj_root '/' parm.brain_file];
0153     brain_parm = parm;
0154     brain_parm.brain_file = brain_file;
0155 else
0156     brain_file = [];
0157 end
0158 
0159 % Output head surface file name
0160 [headfile, basename] = change_name(parm.head_file);
0161 
0162 headmodel{1} = change_name([basename '_CSF'  ], Nvertex(1));
0163 headfile    = [proj_root '/' headfile];     % import surface file
0164 headsave{1} = [proj_root '/' headmodel{1}]; % CSF   surface file
0165 
0166 if Nsurf == 3
0167     % 3 shell model
0168     headmodel{2} = change_name([basename '_Skull'], Nvertex(2));
0169     headmodel{3} = change_name([basename '_Scalp'], Nvertex(3));
0170     head_join    = change_name([basename '_3shell'], sum(Nvertex));
0171     head_join   = [proj_root '/' head_join];    % 3 shell file
0172     headsave{2} = [proj_root '/' headmodel{2}]; % Skull surface file
0173     headsave{3} = [proj_root '/' headmodel{3}]; % Scalp surface file
0174 else
0175     % 1 shell model
0176     head_join   = headsave{1};
0177 end
0178 
0179 % Read MRI image size
0180 [Vdim, Vsize] = analyze_hdr_read(analyzefile);
0181 
0182 [V,F,Vs,Fs] = load_face(parm, Vdim, Vsize);
0183 
0184 %%% DEBUG %%%
0185 if plot_mode > 0
0186     debug_plot(V, F, analyzefile);
0187 end
0188 
0189 %
0190 % --- Make CSF surface
0191 %
0192 parm.Nvertex = Nvertex(1);
0193 parm.Nlast = parm.Nlast_csf;
0194 % Add margin to image for Morphology
0195 DW  = 20;
0196 %    Vhead = Vhead + step*DW;
0197 %    Vhead = Vhead - step*DW;
0198 
0199 if isfield(parm,'freesurf_file') && ~isempty(parm.freesurf_file)
0200 %
0201 % --- Import FreeSurfer surface file
0202     if ~exist(headfile,'file')
0203         [Vhead,Fhead] = vb_fs_load_surface(parm.freesurf_file);
0204         [Fhead,Vhead,XXhead] = vb_out_normal_surf(Fhead,Vhead);
0205         Vhead = Vhead * 0.001;
0206         
0207         disp(['--- Save CSF model file: ']);
0208         disp([headfile]);
0209         vb_fsave(headfile,'Vhead','Fhead','XXhead');
0210     else
0211         load(headfile,'Vhead','Fhead','XXhead');
0212     end
0213     
0214     %%% DEBUG %%%
0215     if plot_mode > 0
0216         debug_plot(Vhead, Fhead, analyzefile);
0217     end
0218     
0219     if ~isempty(brain_file)
0220         % Make mask from brain
0221         B = vb_make_head_surf_model(brain_parm);
0222         Dim  = size(B);
0223         
0224         Bcsf = zeros(Dim + 2*DW);
0225         Bcsf((1:Dim(1))+DW, (1:Dim(2))+DW, (1:Dim(3))+DW) = B;
0226         Dim = Dim + 2*DW;
0227     else
0228         Dim = ceil((Vdim .* Vsize)/step) + 2*DW;
0229         %Dim   = vb_mask_image_size(Vdim,Vsize,step);
0230     end
0231     
0232     % Make mask from FreeSurfer head surface
0233     Vhead = vb_spm_right_to_analyze_right_mm(Vhead, Vdim, Vsize);
0234     Vhead = Vhead + step*DW;
0235     
0236     % Make mask from FreeSurfer surface
0237     B = vb_surf_to_filled_mask(Vhead, Fhead, Dim, step);
0238     
0239     % Smooth & Shrink FreeSurfer surface
0240     B = vb_morphology_operation(B, parm.Radius_fs, step);
0241 
0242     % Merge cortex expand surface & FreeSurfer surface mask
0243     if ~isempty(brain_file)
0244         % Make mask from brain
0245         %Bcsf = vb_make_head_surf_model(brain_parm);
0246         
0247         if vb_matlab_version > 6,
0248             B = int8(B) + int8(Bcsf);
0249         else
0250             B = double(B) + double(Bcsf);
0251         end
0252     end
0253     
0254     % Make mask image from face surface and shrinks it by (Scalp+Skull) radius
0255     if ~isempty(Vs)
0256         Vs     = Vs + step*DW;
0257         Bface  = vb_face_to_mask(Vs, Fs, Dim, step);
0258         Radius = [- parm.Radius_scalp, - parm.Radius_skull];
0259         Bface  = vb_morphology_operation(Bface, Radius , step);
0260         
0261         % Intersection of CSF and shrinked face
0262         if vb_matlab_version > 6,
0263             B = int8(B) .* int8(Bface);
0264         else
0265             B = double(B) .* double(Bface);
0266         end
0267     end
0268     
0269     % Smooth surface with vertex number 'parm.Nvertex'
0270     %   ( Morphology by parm.Radius )
0271     [Vhead, Fhead, XXhead, B] = vb_mask_to_surf_expand(B, parm);
0272     Vhead = Vhead - step*DW;
0273     
0274     Vhead = vb_analyze_right_mm_to_spm_right(Vhead, Vdim, Vsize);
0275     
0276     disp(['--- Save CSF model file: ']);
0277     disp([headsave{1}]);
0278     vb_fsave(headsave{1},'Vhead','Fhead','XXhead');
0279     
0280     if Nsurf==1
0281         vb_plot_slice_head(headsave{1}, analyzefile, brain_file);
0282     end
0283 elseif isfield(parm,'curry_file') && ~isempty(parm.curry_file)
0284 %
0285 % ---     Import Curry surface file
0286     if ~exist(headfile,'file')
0287         vb_job_curry2vbmeg(parm.curry_file, headfile, analyzefile);
0288     end
0289     
0290     % Smooth Curry surface with vertex number 'parm.Nvertex'
0291     if Nsurf==3, parm.plot_mode = 0; end;
0292     
0293     vb_job_smooth_head(headfile, headsave{1}, parm);
0294     
0295 elseif    isfield(parm,'gray_file') && ~isempty(parm.gray_file)
0296 %
0297 % ---  Make mask from brain surface and gray image
0298     [B] = vb_make_head_surf_model(brain_parm);
0299     
0300     %Dim   = vb_mask_image_size(Vdim,Vsize,step);
0301     Dim   = size(B);
0302     
0303     % Make mask image from face surface and shrinks it by (Scalp+Skull) radius
0304     if ~isempty(Vs)
0305         Bface = vb_face_to_mask(Vs, Fs, Dim, step);
0306         Radius = [- parm.Radius_scalp, - parm.Radius_skull];
0307         Bface = vb_morphology_operation(Bface, Radius , step);
0308         
0309         % Intersection of CSF and shrinked face
0310         if vb_matlab_version > 6,
0311             B = int8(B) .* int8(Bface);
0312         else
0313             B = double(B) .* double(Bface);
0314         end
0315     end
0316 
0317     % Smooth surface with vertex number 'parm.Nvertex'
0318     % Morphology by parm.Radius
0319     [Vhead, Fhead, XXhead, B] = vb_mask_to_surf_expand(B, parm);
0320     Vhead = vb_analyze_right_mm_to_spm_right(Vhead, Vdim, Vsize);
0321 
0322     if isempty(XXhead), return; end;
0323     
0324     disp(['--- Save CSF model file: ']);
0325     disp([headsave{1}]);
0326     vb_fsave(headsave{1},'Vhead','Fhead','XXhead');
0327     
0328     if Nsurf==1,
0329         vb_plot_slice_head(headsave{1}, analyzefile, brain_file);
0330     end
0331 else
0332     error('No CSF input file');
0333 end
0334 
0335 if Nsurf==1, 
0336     % 1 shell case
0337     project_file_saving(parm);
0338     return; 
0339 end;
0340 
0341 if isempty(V), error('Scalp file is needed for 3 shell model'); end
0342 
0343 % Load CSF surface
0344 load(headsave{1},'Vhead','Fhead')
0345 
0346 %%% DEBUG %%%
0347 if plot_mode > 0
0348     debug_plot(Vhead, Fhead, analyzefile);
0349     if plot_mode == 1, return; end;
0350 end
0351 
0352 %
0353 % --- Make Scalp : intersection of face surface and dilated CSF surface
0354 %
0355 [Rmax,Rmin] = vb_distance_brain_face(Vhead,V);
0356 
0357 fprintf('Max distance from face to CSF = %5.1f\n',Rmax)
0358 fprintf('Min distance from face to CSF = %5.1f\n',Rmin)
0359 
0360 Rstep = round(Rmax/2);
0361 Radius = [Rstep  ceil(Rmax + 2 - Rstep)];
0362 fprintf('Expansion Radius from CSF to Scalp = %3.0f\n',sum(Radius))
0363 
0364 % --- Make intersection of dilated brain and face surface
0365 
0366 [B ,Bcsf, Bface] = ...
0367     vb_intersect_face_csf(Vhead,Fhead,V,F,Radius,step,Vdim,Vsize,DW);
0368 
0369 %%% DEBUG %%%
0370 if plot_mode > 1
0371     zindx = [40:20:200];
0372     zindx = ceil(zindx/step);
0373     vb_plot_slice( B, [], zindx, 1);
0374     vb_plot_slice( Bcsf, [], zindx, 1);
0375     vb_plot_slice( Bface, [], zindx, 1);
0376     if plot_mode == 2, return; end;
0377 end
0378 
0379 parm.Nlast = parm.Nlast_scalp;
0380 
0381 if parm.Scalp_mode == 0
0382     % --- Extract smooth surface from intersect_face image
0383     parm.Nvertex = Nvertex(3);
0384     parm.Radius  = parm.Radius_final;
0385     
0386     [Vhead, Fhead, XXhead, B] = vb_mask_to_surf_expand(B, parm);
0387     
0388     Vhead = Vhead - step*DW;
0389     Vhead = vb_analyze_right_mm_to_spm_right(Vhead,Vdim,Vsize);
0390     
0391     disp(['--- Save Scalp model file: ']);
0392     disp([headsave{3}]);
0393     vb_fsave(headsave{3},'Vhead','Fhead','XXhead');
0394 else
0395     % --- Extract closed face surface
0396     parm.Radius  = parm.Radius_final;
0397     parm.Nvertex = Nvertex(3);
0398     
0399     [Vhead, Fhead, XXhead, Bface] = vb_mask_to_surf_expand(Bface, parm);
0400     
0401     Vhead = Vhead - step*DW;
0402     Vhead = vb_analyze_right_mm_to_spm_right(Vhead,Vdim,Vsize);
0403     
0404     disp(['--- Save Face model file: ']);
0405     disp([headsave{3}]);
0406     vb_fsave(headsave{3},'Vhead','Fhead','XXhead');
0407     
0408     switch    parm.Scalp_mode 
0409     case    1
0410         % --- smooth Scalp surface
0411         parm.Radius  = parm.Radius_final;
0412         B = vb_morphology_operation(B, parm.Radius, step);
0413     case    2
0414         % --- smooth Face surface
0415         B = Bface;
0416     end
0417 end
0418 
0419 clear Bface
0420 
0421 %%% DEBUG %%%
0422 if plot_mode > 2
0423     debug_plot(Vhead, Fhead, analyzefile);
0424     if plot_mode == 3, return; end;
0425 end
0426 %
0427 % --- Make Outer Skull surface : erosion from Scalp surface
0428 %
0429 
0430 if isfield(parm,'fs_skull_file') && ~isempty(parm.fs_skull_file)
0431 %
0432 % --- intersection of Skull surface and dilated CSF surface
0433 %
0434     % Load Skull & CSF surface
0435     load(headsave{1},'Vhead','Fhead')
0436     [V,F] = vb_fs_load_surface(parm.fs_skull_file);
0437     V = V * 0.001;
0438     
0439     % CSF Expansion Radius
0440     Radius(2) = Radius(2) - parm.Radius_scalp;
0441     
0442     % --- Make intersection of dilated CSF and Skull surface mask
0443     [B ,Btmp, Bskl] = ...
0444         vb_intersect_face_csf(Vhead,Fhead,V,F,Radius,step,Vdim,Vsize,DW);
0445     
0446     clear Btmp
0447     switch    parm.Scalp_mode 
0448     case    {0,1}
0449         % --- smooth Skull surface
0450         parm.Radius  = parm.Radius_final;
0451         B = vb_morphology_operation(B, parm.Radius, step);
0452     case    2
0453         % --- smooth Face surface
0454         B = Bskl;
0455     end
0456 else    
0457     % --- Erosion from Scalp surface
0458     fprintf('Max width of Scalp = %5.1f\n',parm.Radius_scalp)
0459     
0460     B = vb_morphology_operation(B, - parm.Radius_scalp , step);
0461     %fprintf('# of Scalp mask %d\n',length(find(B(:)>0)))
0462 end
0463 % --- Set min dilation radius from CSF
0464 Rskull = Rmin * parm.Skull_ratio;
0465 fprintf('Min width of Skull = %5.1f\n',Rskull)
0466 
0467 Bcsf = vb_morphology_operation(Bcsf, Rskull, step);
0468 %fprintf('# of CSF   mask %d\n',length(find(Bcsf(:)>0)))
0469 
0470 % Make union of two images
0471 if vb_matlab_version > 6,
0472     B  = int8(B) + int8(Bcsf);
0473     ix = find( B(:) > 0 );
0474     B  = zeros(size(B),'int8');
0475     B(ix) = 1;
0476 else
0477     B  = double(B) + double(Bcsf);
0478     ix = find( B(:) > 0 );
0479     B  = zeros(size(B));
0480     B(ix) = 1;
0481 end
0482 
0483 %fprintf('# of Skull   mask %d\n',length(ix))
0484 
0485 % --- Extract smooth Skull surface from mask image [analyze_right_mm]
0486 parm.Radius  = parm.Radius_final;
0487 parm.Nvertex = Nvertex(2);
0488 parm.Nlast   = parm.Nlast_skull;
0489 
0490 [Vhead, Fhead, XXhead, B] = vb_mask_to_surf_expand(B, parm);
0491 
0492 Vhead = Vhead - step*DW;
0493 Vhead = vb_analyze_right_mm_to_spm_right(Vhead,Vdim,Vsize);
0494 
0495 disp(['--- Save Skull model file: ']);
0496 disp([headsave{2}]);
0497 vb_fsave(headsave{2},'Vhead','Fhead','XXhead');
0498 
0499 % join 3 surface into one file
0500 vb_head_join_files(headmodel, proj_root, head_join);
0501 
0502 %
0503 % --- Check 3 surface on MRI
0504 %
0505 
0506 vb_plot_slice_head(headsave, analyzefile, brain_file);
0507 
0508 %return
0509 
0510 % --- Check Scalp surface with face
0511 load(headsave{3},'Vhead','Fhead');
0512 
0513 if exist('surf_face','var') && isfield(surf_face,'V_reduce')
0514     V = surf_face.V_reduce;
0515     F = surf_face.F_reduce;
0516 end
0517 
0518 figure
0519 subplot(1,2,1)
0520 vb_plot_surf(V,F,[],[],1,1);
0521 hold on
0522 vb_plot_surf(Vhead,Fhead,'none','k');
0523 alpha 0.8;
0524 
0525 title('Scalp surface on face')
0526 
0527 view([ -210 30]);
0528 
0529 subplot(1,2,2)
0530 vb_plot_surf(Vhead,Fhead,[],[],1,1);
0531 
0532 title('Scalp surface')
0533 
0534 view([ -210 30]);
0535 
0536 project_file_saving(parm);
0537 
0538 return
0539 
0540 function    project_file_saving(parm)
0541 %
0542 % --- Save project file
0543 %
0544 proj_file = get_project_filename;
0545 if isempty(proj_file)
0546     return;
0547 end
0548 project_file_mgr('load', proj_file);
0549 project_file_mgr('add', 'head_parm', parm);
0550 
0551 
0552 return
0553 
0554 function    [fnew, basename] = change_name(fold, N);
0555 %
0556 % --- change_name
0557 %
0558 
0559 fext = '.head.mat';
0560 ix_ext = findstr(fold, fext);
0561 
0562 if ~isempty(ix_ext),
0563     basename = fold(1:ix_ext-1);
0564 else
0565     basename = fold;
0566 end
0567 
0568 if nargin == 2
0569     fnew = sprintf('%s_%d%s', basename, N, fext);
0570 else
0571     fnew = sprintf('%s%s', basename, fext);
0572 end
0573 
0574 return
0575 
0576 function    [V,F,Vs,Fs] = load_face(parm, Vdim, Vsize);
0577 
0578 % Load Face surface
0579 if isfield(parm,'fs_scalp_file') && ~isempty(parm.fs_scalp_file)
0580     [V,F] = vb_fs_load_surface(parm.fs_scalp_file);
0581     V = V * 0.001; % [mm] -> [m]
0582 elseif isfield(parm,'face_file') && ~isempty(parm.face_file)
0583     load(parm.face_file,'surf_face');
0584     V = surf_face.V;
0585     F = surf_face.F;
0586 else
0587     V  = []; F  = []; 
0588     Vs = []; Fs = []; 
0589     return
0590 end
0591 
0592 % Make surface closed
0593 [F,V] = vb_close_surf(F,V);
0594 
0595 % Coordinate transform to analyze_right_mm
0596 Vs = vb_spm_right_to_analyze_right_mm(V, Vdim, Vsize);
0597 Fs = F;
0598 
0599 % Min z-coordinate for face mask [mm]
0600 zmin = 2; 
0601 Vs(:,3) = max(Vs(:,3),zmin);
0602 
0603 return
0604 
0605 function    debug_plot(V, F, analyzefile)
0606     [B, Vdim, Vsize] = vb_load_analyze_to_right(analyzefile);
0607     Vana = vb_spm_right_to_analyze_right( V, Vdim, Vsize);
0608     
0609     zindx = 0.1*[2:2:8];
0610     xindx = 0.1*[2:2:8];
0611     
0612     Vmin = min(Vana);
0613     Vmax = max(Vana);
0614     zindx = ceil( zindx*(Vmax(3) - Vmin(3)) + Vmin(3));
0615     xindx = ceil( xindx*(Vmax(1) - Vmin(1)) + Vmin(1));
0616     zindx = zindx( zindx > 0 & zindx <= Vdim(3) );
0617     xindx = xindx( xindx > 0 & xindx <= Vdim(1) );
0618     
0619     vb_plot_slice_surf(B, Vana, F, zindx, 'z', [2 2]);
0620     vb_plot_slice_surf(B, Vana, F, xindx, 'x', [2 2]);
0621 
0622 return
0623 
0624 ix = find( Vs(:,3) >= zmin);
0625 Fs = vb_patch_select2(ix,F,size(Vs,1));

Generated on Tue 27-Aug-2013 11:46:04 by m2html © 2005