Extract face surface from MRI image
[F, V] = vb_face_extract(imagefile, Radius, step)
[F, V] = vb_face_extract(imagefile, Radius, step, pmax, Bval, hstep)
---- Output
V : Face vertex coordinate: (x,y,z) [Nvertex x 3]
F : Triangle patch index [Npatch x 3]
---- Input
imagefile : MRI Analyze/NIFTI image file name
----- Optional input
Radius = [ -2 2 6 -6]; [mm]
: Radius of Morfology operation
R > 0 : dilation
R < 0 : erosion
step : subsampling step size [mm]
pmax : Prob. to detrmine Threshold value for face extraction (0.998)
Bval : Threshold value for background separation
This value is automatically determined in the default mode
hstep : step size to round the value of image intensity
--------------------------------------------------
MRI構造画像から顔表面画像の抽出
V : 顔表面頂点 face vertex
F : その三角面パッチ patch index
Masa-aki Sato 2010-10-18
Image intensity is rounded by deviding 'hstep'
before making mask to avoid noisy fluctuation
2011-06-20 taku-y
[minor] Progress message was added.
Copyright (C) 2011, ATR All Rights Reserved.
License : New BSD License(see VBMEG_LICENSE.txt)0001 function [F, V] = vb_face_extract(imagefile, Radius, step, pmax, Bval, hstep) 0002 % Extract face surface from MRI image 0003 % [F, V] = vb_face_extract(imagefile, Radius, step) 0004 % [F, V] = vb_face_extract(imagefile, Radius, step, pmax, Bval, hstep) 0005 % 0006 % ---- Output 0007 % V : Face vertex coordinate: (x,y,z) [Nvertex x 3] 0008 % F : Triangle patch index [Npatch x 3] 0009 % ---- Input 0010 % imagefile : MRI Analyze/NIFTI image file name 0011 % 0012 % ----- Optional input 0013 % Radius = [ -2 2 6 -6]; [mm] 0014 % : Radius of Morfology operation 0015 % R > 0 : dilation 0016 % R < 0 : erosion 0017 % 0018 % step : subsampling step size [mm] 0019 % pmax : Prob. to detrmine Threshold value for face extraction (0.998) 0020 % Bval : Threshold value for background separation 0021 % This value is automatically determined in the default mode 0022 % hstep : step size to round the value of image intensity 0023 % -------------------------------------------------- 0024 % MRI構造画像から顔表面画像の抽出 0025 % V : 顔表面頂点 face vertex 0026 % F : その三角面パッチ patch index 0027 % 0028 % Masa-aki Sato 2010-10-18 0029 % Image intensity is rounded by deviding 'hstep' 0030 % before making mask to avoid noisy fluctuation 0031 % 2011-06-20 taku-y 0032 % [minor] Progress message was added. 0033 % 0034 % Copyright (C) 2011, ATR All Rights Reserved. 0035 % License : New BSD License(see VBMEG_LICENSE.txt) 0036 0037 if nargin < 2, Radius = [ -2 2 ]; end; 0038 if nargin < 3, step = 2; end; 0039 if nargin < 4, pmax = 0.998; end; 0040 if nargin < 5, Bval = []; end; 0041 if nargin < 6, hstep = 1; end; 0042 0043 h = waitbar(0,'Wait for face extraction'); 0044 vb_disp_nonl(' 0 % processed'); 0045 0046 % 0047 % Load 3D MRI image to analyze_right_hand coordinate 0048 % 0049 % Vdim : Voxel dimension of Analyze image 0050 % Vsize : Voxel size of Analyze image 0051 % --- Analyze righthand voxel coordinate 0052 % X: Left(1) -> Right (Vdim(1)) 0053 % Y: Back(1) -> Front (Vdim(2)) 0054 % Z: Bottom(1) -> Top (Vdim(3)) 0055 [B, Vdim, Vsize] = vb_load_analyze_to_right(imagefile); 0056 0057 waitbar(0.1); 0058 for ii=1:15; vb_disp_nonl(sprintf('\b')); end 0059 vb_disp_nonl(' 10 % processed'); 0060 0061 % Calculate background noise threshold value 0062 % by estimating backgound distribution 0063 if isempty(Bval) 0064 [Bval, hstep] = vb_back_threshold(B, pmax); 0065 end 0066 0067 waitbar(0.2); 0068 for ii=1:15; vb_disp_nonl(sprintf('\b')); end 0069 vb_disp_nonl(' 20 % processed'); 0070 0071 B = round(B/hstep); 0072 Bval = round(Bval/hstep); 0073 0074 % 閾値より大きい値を持つボクセルを頭部として抽出 0075 % Make mask image from MRI data 0076 Bmask = vb_image_to_mask(B, Bval, step); 0077 0078 if vb_matlab_version > 6, Bmask = int8(Bmask); end; 0079 0080 % Expand Z-axis to avoid top head reaches to image boundary 0081 Zex = ceil(10/step); % 10 mm 0082 [Nx,Ny,Nz]=size(Bmask); 0083 0084 % vb_disp_nonl(sprintf('\nExpand image')); 0085 0086 if vb_matlab_version > 6, 0087 B = zeros(Nx,Ny,Nz + Zex, 'int8'); 0088 B(1:Nx,1:Ny,1:Nz) = Bmask; 0089 else 0090 B = zeros(Nx,Ny,Nz + Zex); 0091 B(1:Nx,1:Ny,1:Nz) = Bmask; 0092 end 0093 0094 waitbar(0.4); 0095 for ii=1:15; vb_disp_nonl(sprintf('\b')); end 0096 vb_disp_nonl(' 40 % processed'); 0097 0098 % Apply morphology operations 0099 % Dilation/erosion are done consecutively according to Radius 0100 B = vb_morphology_operation(B, Radius, step); 0101 0102 waitbar(0.8); 0103 for ii=1:15; vb_disp_nonl(sprintf('\b')); end 0104 vb_disp_nonl(' 80 % processed'); 0105 0106 % face extraction 0107 % V : 顔表面頂点 face vertex 0108 % F : その三角面パッチ patch index 0109 0110 [F, V] = vb_surf_extract(B, step); 0111 0112 waitbar(0.95); 0113 for ii=1:15; vb_disp_nonl(sprintf('\b')); end 0114 vb_disp_nonl(' 95 % processed'); 0115 0116 % Change Analyze Left-hand voxcel coordinate to Right-hand SPM (m) coord. 0117 % Transformation from Analyze to SPM is not changed 0118 % despite of mask image expansion along Z-axis 0119 % since SPM is defined in original image 0120 V = vb_analyze_right_to_spm_right(V,Vdim,Vsize); 0121 0122 close(h); 0123 drawnow;