Home > vbmeg > functions > common > boundary > vb_surf_smooth_expand.m

vb_surf_smooth_expand

PURPOSE ^

Expand surface by morphology and smooth surface by spring model

SYNOPSIS ^

function [Vx, Fx, xx, B] = vb_surf_smooth_expand(V, F, xx, Para)

DESCRIPTION ^

 Expand surface by morphology and smooth surface by spring model
  [V_out, F_out, xx_out, B] = vb_surf_smooth_expand(V, F, xx, Para)
 --- Input & Output
  V(n, 1:3)  : vertex
  F(j, 1:3)  : patch index
 xx(n, 1:3)  : normal vector
 B(NX,NY,NZ) : mask image : inside of surface is filled out
 --- Parameter
 Para : parameter structure ( = default value )
 Para.Dim     : Dimension of mask image
 Para.Nvertex = 5000 : # of vertex
 Para.Nloop  = 200 : iteration number for fit boundary
 Para.Nlast  = 10  : iteration number for smoothing by spring model
 Para.tangent_rate   = 0.3  : spring constant
 Para.mask_ratio     = 0.5  : mask image force constant
 Para.mask_threshold = 0.3  : mask threthold
 Para.vstep  = 2   : voxel size of mask image [mm]
 Para.Radius = [ 6 -6 ] : morphological smoothing [mm]
             = [ 6 ]    : morphological expansion
             = []       : No morphological operation
 --- 膨張モデル
 Para.tangent_rate   = 0.3 :  バネ力係数
 Para.mask_ratio     = 0.5 : マスク強度力係数 ( 強度>閾値:外向き)
 Para.mask_threshold = 0.3 : マスク強度閾値   ( 強度<閾値:内向き)

 合力 =  tangent_rate * (平滑化バネ力)
       + mask_ratio   * (外向き法線方向) * (マスク強度 - 閾値)
 --- バネモデル平滑化

 合力 =  tangent_rate * (平滑化バネ力)

 2006/10/12 M. Sato
 2006/11/11 M. Sato
 M. Sato 2007-3-16

 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:

SOURCE CODE ^

0001 function    [Vx, Fx, xx, B] = vb_surf_smooth_expand(V, F, xx, Para)
0002 % Expand surface by morphology and smooth surface by spring model
0003 %  [V_out, F_out, xx_out, B] = vb_surf_smooth_expand(V, F, xx, Para)
0004 % --- Input & Output
0005 %  V(n, 1:3)  : vertex
0006 %  F(j, 1:3)  : patch index
0007 % xx(n, 1:3)  : normal vector
0008 % B(NX,NY,NZ) : mask image : inside of surface is filled out
0009 % --- Parameter
0010 % Para : parameter structure ( = default value )
0011 % Para.Dim     : Dimension of mask image
0012 % Para.Nvertex = 5000 : # of vertex
0013 % Para.Nloop  = 200 : iteration number for fit boundary
0014 % Para.Nlast  = 10  : iteration number for smoothing by spring model
0015 % Para.tangent_rate   = 0.3  : spring constant
0016 % Para.mask_ratio     = 0.5  : mask image force constant
0017 % Para.mask_threshold = 0.3  : mask threthold
0018 % Para.vstep  = 2   : voxel size of mask image [mm]
0019 % Para.Radius = [ 6 -6 ] : morphological smoothing [mm]
0020 %             = [ 6 ]    : morphological expansion
0021 %             = []       : No morphological operation
0022 % --- 膨張モデル
0023 % Para.tangent_rate   = 0.3 :  バネ力係数
0024 % Para.mask_ratio     = 0.5 : マスク強度力係数 ( 強度>閾値:外向き)
0025 % Para.mask_threshold = 0.3 : マスク強度閾値   ( 強度<閾値:内向き)
0026 %
0027 % 合力 =  tangent_rate * (平滑化バネ力)
0028 %       + mask_ratio   * (外向き法線方向) * (マスク強度 - 閾値)
0029 % --- バネモデル平滑化
0030 %
0031 % 合力 =  tangent_rate * (平滑化バネ力)
0032 %
0033 % 2006/10/12 M. Sato
0034 % 2006/11/11 M. Sato
0035 % M. Sato 2007-3-16
0036 %
0037 % Copyright (C) 2011, ATR All Rights Reserved.
0038 % License : New BSD License(see VBMEG_LICENSE.txt)
0039 
0040 if ~exist('Para','var'), Para = []; end
0041 
0042 if ~isfield(Para,'Dim') | isempty(Para.Dim)
0043     error('Size of mask image is not specified');
0044 else
0045     Dim = Para.Dim;
0046 end
0047 % --- subsampling step for mask image [mm]
0048 % 3D 画像処理のための間引きステップサイズ
0049 if ~isfield(Para,'vstep'), 
0050     step = 2; 
0051     Para.vstep = step;
0052 else
0053     step = Para.vstep;
0054 end;
0055 
0056 if ~isfield(Para,'Nloop'), Para.Nloop = 200; end
0057 if ~isfield(Para,'Nlast'), Para.Nlast = 20; end
0058 if ~isfield(Para,'Nvertex'), Para.Nvertex = 5000; end
0059 if ~isfield(Para,'tangent_rate'), Para.tangent_rate = 0.3; end
0060 if ~isfield(Para,'Radius'),Para.Radius  = [6 -6];end 
0061 if ~isfield(Para,'plot_mode'), Para.plot_mode = 0; end
0062 
0063 
0064 %%% DEBUG
0065 Debug_mode = 0;
0066 
0067 %
0068 % Closed surface check
0069 %
0070 disp('--- Closed surface check');
0071 omega = vb_solid_angle_check(V,F);
0072 fprintf('Solid angle = 4pi x %f\n',omega);
0073 if abs(omega-1)>1e-10, 
0074   disp('Surface is not closed.');
0075   return;
0076 else
0077   disp('Surface is closed.');
0078 end
0079 
0080 %
0081 %-------- 表面をマスクに変換
0082 %
0083 fprintf('--- Make mask image from surface\n')
0084 dmin  = step/2;    % 三角形分割長さ
0085 
0086 B = vb_surf_to_mask(V, F, step, dmin, Dim);
0087 
0088 %
0089 %------- 内部点を塗りつぶす
0090 %
0091 fprintf('--- Fill out inside the surface \n')
0092 
0093 filval = 1;        % 内部点マスクの値
0094 level  = 0.5;    % この値より小さい値を塗りつぶす
0095 
0096 X0  = fix(mean(V/step));    % 中心点
0097 B   = vb_flood_fill_3d(B, X0, filval, level);
0098 
0099 if Debug_mode == 2, return; end
0100 %
0101 % --- Apply morphology erosion
0102 %
0103 B = vb_morphology_operation(B, Para.Radius, step);
0104 
0105 if Debug_mode == 1, return; end
0106 %
0107 % --- 3次元データの平滑化
0108 %
0109 tic
0110 fprintf('smoothing\n')
0111 B = smooth3(B,'gaussian',3);
0112 vb_ptime(toc);
0113 
0114 %
0115 %------- マスクパターンの境界面に膨張
0116 %
0117 tic
0118 fprintf('--- Surface expansion\n')
0119 
0120 [Vx, Fx, xx] = vb_make_inner_sphere(V, Para.Nvertex);% 内接球面作成
0121 [Vx, Fx, xx] = vb_surf_smooth_fit(Vx,Fx,xx, B, Para);% 境界面に膨張
0122 
0123 vb_ptime(toc)
0124 
0125 if Para.plot_mode > 1
0126     dmax = 10;
0127     vb_plot_slice_surf( B, Vx/step, Fx, X0(1), 'x',[1 1],'r-',dmax);
0128     vb_plot_slice_surf( B, Vx/step, Fx, X0(2), 'y',[1 1],'r-',dmax);
0129     vb_plot_slice_surf( B, Vx/step, Fx, X0(3), 'z',[1 1],'r-',dmax);
0130 end
0131 
0132 %
0133 %------- バネモデル による平滑化
0134 %
0135 tic
0136 fprintf('--- Surface smoothing by spring model\n')
0137 
0138 Para.Nloop = Para.Nlast;     % バネモデル平滑化繰り返し回数
0139 [Vx, Fx, xx] = vb_surf_smooth(Vx,Fx,xx,Para);
0140 
0141 vb_ptime(toc)
0142 
0143 fprintf('# of patches : %d\n',size(Fx,1));
0144 fprintf('# of vertices: %d\n',size(Vx,1));
0145 
0146 if Para.plot_mode > 1
0147     dmax = 10;
0148     vb_plot_slice_surf( B, Vx/step, Fx, X0(1), 'x',[1 1],'r-',dmax);
0149     vb_plot_slice_surf( B, Vx/step, Fx, X0(2), 'y',[1 1],'r-',dmax);
0150     vb_plot_slice_surf( B, Vx/step, Fx, X0(3), 'z',[1 1],'r-',dmax);
0151 end
0152

Generated on Mon 22-May-2023 06:53:56 by m2html © 2005