Description of vb_job_head

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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
