vb_surf_smooth_fit

 Fit surface to the boundary of the mask
  [V, F, xx] = vb_surf_smooth_fit(V,F,xx,B,Para)
 --- Input
  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
               surface is expanded to fit to the boundary of mask region
 Para.Nloop  : iteration number
 Para.tangent_rate   : spring constant
 Para.mask_ratio     : mask image force constant
 Para.mask_threshold : mask threthold

 --- Output
  V(n, 1:3)  : vertex
  F(j, 1:3)  : patch index
 xx(n, 1:3)  : normal vector
 --- Mehod
 ポリゴンモデルをマスクパターンの境界面に膨張、バネ力による平滑化を行う

  V(n, 1:3)  : 頂点の位置
  F(j, 1:3)  : 三角面３頂点のインデックス
 xx(n, 1:3)  : 頂点の初期法線方向(nx,ny,nz)
 xxn(n, 1:3) : 頂点の最終法線方向(nx,ny,nz)
 B(NX,NY,NZ) : 3D-ボクセル・マスクパターン

 Para.Nloop  : 繰り返し回数
 Para.vstep  : ボクセルサイズ

 Para.tangent_rate   : バネ強度
 Para.mask_ratio     : マスク強度力係数 ( 強度>閾値:外向き)
 Para.mask_threshold : マスク強度閾値   ( 強度<閾値:内向き)

 change coordinate to index
     V = [-1/2,-1/2,-1/2] - [1/2,1/2,1/2]
 <=> J = [1,1,1]
   J   = floor(V + 0.5) + 1;

 合力 =  tangent_rate * (平滑化バネ力)
       + mask_ratio   * (外向き法線方向) * (マスク強度 - 閾値)
 
 Ver 1.0  by M. Sato  2004-2-5
 Ver 2.0  by M. Sato  2006-11-11
 M. Sato 2007-3-16

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

0001 function    [V, F, xxn] = vb_surf_smooth_fit(V,F,xx,B,Para)
0002 % Fit surface to the boundary of the mask
0003 %  [V, F, xx] = vb_surf_smooth_fit(V,F,xx,B,Para)
0004 % --- Input
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 %               surface is expanded to fit to the boundary of mask region
0010 % Para.Nloop  : iteration number
0011 % Para.tangent_rate   : spring constant
0012 % Para.mask_ratio     : mask image force constant
0013 % Para.mask_threshold : mask threthold
0014 %
0015 % --- Output
0016 %  V(n, 1:3)  : vertex
0017 %  F(j, 1:3)  : patch index
0018 % xx(n, 1:3)  : normal vector
0019 % --- Mehod
0020 % ポリゴンモデルをマスクパターンの境界面に膨張、バネ力による平滑化を行う
0021 %
0022 %  V(n, 1:3)  : 頂点の位置
0023 %  F(j, 1:3)  : 三角面３頂点のインデックス
0024 % xx(n, 1:3)  : 頂点の初期法線方向(nx,ny,nz)
0025 % xxn(n, 1:3) : 頂点の最終法線方向(nx,ny,nz)
0026 % B(NX,NY,NZ) : 3D-ボクセル・マスクパターン
0027 %
0028 % Para.Nloop  : 繰り返し回数
0029 % Para.vstep  : ボクセルサイズ
0030 %
0031 % Para.tangent_rate   : バネ強度
0032 % Para.mask_ratio     : マスク強度力係数 ( 強度>閾値:外向き)
0033 % Para.mask_threshold : マスク強度閾値   ( 強度<閾値:内向き)
0034 %
0035 % change coordinate to index
0036 %     V = [-1/2,-1/2,-1/2] - [1/2,1/2,1/2]
0037 % <=> J = [1,1,1]
0038 %   J   = floor(V + 0.5) + 1;
0039 %
0040 % 合力 =  tangent_rate * (平滑化バネ力)
0041 %       + mask_ratio   * (外向き法線方向) * (マスク強度 - 閾値)
0042 %
0043 % Ver 1.0  by M. Sato  2004-2-5
0044 % Ver 2.0  by M. Sato  2006-11-11
0045 % M. Sato 2007-3-16
0046 %
0047 % Copyright (C) 2011, ATR All Rights Reserved.
0048 % License : New BSD License(see VBMEG_LICENSE.txt)
0049 
0050 if ~exist('Para','var'), Para = []; end
0051 if ~isfield(Para,'Nloop'), Para.Nloop = 200; end
0052 if ~isfield(Para,'Ndisp'), Para.Ndisp = Para.Nloop+1; end
0053 if ~isfield(Para,'vstep'), Para.vstep = 1; end
0054 if ~isfield(Para,'tangent_rate'),   Para.tangent_rate   = 0.3; end
0055 if ~isfield(Para,'mask_ratio'),     Para.mask_ratio     = 0.5; end
0056 if ~isfield(Para,'mask_threshold'), Para.mask_threshold = 0.3; end
0057 
0058 tangent_rate = Para.tangent_rate;
0059 image_rate   = Para.mask_ratio;
0060 threshold    = Para.mask_threshold;
0061 step         = Para.vstep;
0062 
0063 Nloop         = Para.Nloop;
0064 Npoint         = size(V,1);      % number of vertex
0065 Npatch         = size(F,1);      % number of patch
0066 [NX,NY,NZ]     = size(B);        % image size
0067 
0068 % Plot parameter
0069 Ndisp = Para.Ndisp;
0070 Nfig  = fix(Nloop/Ndisp);
0071 
0072 if Nfig < 2,
0073     NYfig = 1;
0074 else
0075     NYfig = 2;
0076 end
0077 NXfig = max(ceil(Nfig/NYfig), 1);
0078 nfig  = NXfig*NYfig + 1;
0079 
0080 fclr = [];
0081 eclr = [];
0082 light_mode = 1;
0083 vangle = [-70 20];
0084 
0085 % ３角面の法線ベクトル
0086 xxf  = zeros(Npatch,3);
0087 % ３角面面積
0088 nns  = zeros(Npatch,1);
0089 % ３角面頂点の差分ベクトル
0090 VV1  = zeros(Npatch,3);
0091 VV2  = zeros(Npatch,3);
0092 VV3  = zeros(Npatch,3);
0093 
0094 % 頂点法線 = 頂点に隣接する三角面法線の平均
0095 xxn  = zeros(Npoint,3);    
0096 xxs  = zeros(Npoint,1);
0097 % 頂点差分ベクトル
0098 fd   = zeros(Npoint,3);
0099 % 各頂点の近傍点数
0100 Nv   = zeros(Npoint,1);
0101 % 頂点差分ベクトルと法線の内積
0102 nnf  = zeros(Npoint,1);
0103 % 接線方向の力
0104 ft     = zeros(Npoint,3);
0105 % 法線方向の力
0106 fn     = zeros(Npoint,3);
0107 % 外向き(内側)・内向き(外側)のイメージ強度による力
0108 fg     = zeros(Npoint,3);
0109 
0110 % V 座標のボクセルインデックス
0111 jv     = zeros(Npoint,3);
0112 jj     = zeros(Npoint,1);
0113 % 内側マスク強度
0114 gg     = zeros(Npoint,1);
0115 
0116 % 三角面３頂点のインデックス
0117 F1     = F(:,1);
0118 F2     = F(:,2);
0119 F3     = F(:,3);
0120 
0121 % ３角面の法線ベクトル
0122 xxf  = vb_cross2( V(F2,:)-V(F1,:) , V(F3,:)-V(F1,:) );
0123 
0124 % 三角面３頂点の法線ベクトル平均
0125 xxk  = xx(F1,:)+xx(F2,:)+xx(F3,:);
0126 
0127 xdot = sum( xxk .* xxf ,2);
0128 ix     = find(sign(xdot) < 0 );
0129 
0130 % ３角面の法線ベクトルと
0131 % 頂点の法線ベクトル 'xx' (外向き)の向きをそろえる
0132 F(ix,2) = F3(ix);
0133 F(ix,3) = F2(ix);
0134 
0135 F2        = F(:,2);
0136 F3        = F(:,3);
0137 
0138 % 各頂点の近傍点数
0139 for n=1:Npatch,
0140     inx = F(n,:)';
0141     Nv(inx) = Nv(inx) + 2; 
0142 end;
0143 
0144 % バネ力平滑化＋外向き力(内側)＋内向き力(外側)
0145 
0146 for i=1:Nloop,
0147 
0148     % 三角面３頂点
0149     V1=V(F1,:);
0150     V2=V(F2,:);
0151     V3=V(F3,:);
0152     
0153     % ３角面の法線ベクトル
0154     xxf   = vb_cross2( V2 - V1, V3 - V1 );
0155 
0156     % Normalization
0157     nns   = sqrt(sum(xxf.^2,2));
0158     xxf   = xxf./nns(:,ones(1,3));        
0159     
0160     % 三角面各頂点の差分ベクトル
0161     VV1   = V2 + V3 - 2*V1;
0162     VV2   = V3 + V1 - 2*V2;
0163     VV3   = V1 + V2 - 2*V3;
0164 
0165     % 頂点法線 = 頂点に隣接する三角面法線の平均
0166     xxn   = zeros(Npoint,3);    
0167     % 頂点差分ベクトル(近傍和)
0168     fd      = zeros(Npoint,3);
0169     
0170     for n=1:Npatch,
0171         % 三角面３頂点インデックス
0172         j1=F1(n);
0173         j2=F2(n);
0174         j3=F3(n);
0175         
0176         xxn(j1,:) = xxn(j1,:) + xxf(n,:);
0177         xxn(j2,:) = xxn(j2,:) + xxf(n,:);
0178         xxn(j3,:) = xxn(j3,:) + xxf(n,:);
0179         
0180         % 頂点差分ベクトル(近傍和)
0181         fd(j1,:) = fd(j1,:) + VV1(n,:); 
0182         fd(j2,:) = fd(j2,:) + VV2(n,:); 
0183         fd(j3,:) = fd(j3,:) + VV3(n,:); 
0184     end;
0185     
0186     % 法線の正規化
0187     xxs   = sqrt(sum(xxn.^2,2));
0188     xxn   = xxn./xxs(:,ones(1,3));
0189 
0190     % 頂点差分ベクトルを近傍点数で正規化
0191     fd     = fd./Nv(:,ones(1,3));
0192     
0193     % V 座標をボクセルインデックスに変換
0194     jv    = floor(V/step + 0.5) + 1;
0195     % Check limit
0196     ixv = find(jv(:,1) > 0 & jv(:,1) <= NX ...
0197              & jv(:,2) > 0 & jv(:,2) <= NY ...
0198              & jv(:,3) > 0 & jv(:,3) <= NZ);
0199 
0200     %---- sub2ind ----
0201     jj  = jv(ixv,1) + NX*( (jv(ixv,2) - 1) + NY*(jv(ixv,3) - 1) );
0202 
0203     % 内側マスク強度
0204     gg    = zeros(Npoint,1);
0205     gg(ixv)    = B(jj) - threshold;
0206 %    gg(ixv) = max(B(jj) - threshold, 0);
0207 
0208     fg    = xxn.*gg(:,ones(1,3));
0209     
0210     V    = V + tangent_rate*fd + image_rate*fg;
0211 
0212     if rem(i,Ndisp)==0,
0213         if nfig > NYfig*NXfig
0214             figure;
0215             nfig=1;
0216         else
0217             nfig=nfig+1;
0218         end
0219         subplot(NYfig,NXfig,nfig); 
0220         vb_plot_surf(V,F,fclr,eclr,light_mode);
0221         view(vangle);
0222         tlabel = sprintf('Iteration = %d',i);
0223         title(tlabel);
0224     end;
0225 end
0226