vb_surf_smooth_data

 smoothing surface by spring model and 3D-image data
  [V, F, xxn] = vb_surf_smooth_data(V,F,xx,vg,Para)

 Ver 1.0  by M. Sato  2004-2-10

 ポリゴンモデルをマスクパターンの境界面に膨らませ
 バネ力による平滑化を行う
 平滑化バネ力＋外向き力(内側)／内向き力(外側)

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

 Para.tangent_rate         : 接線方向力係数
 Para.tangent_normal_ratio : 法線方向力係数
 Para.mask_ratio           : 画像強度力係数 ( 強度>閾値:外向き)
 Para.mask_threshold       : 画像強度閾値   ( 強度<閾値:内向き)
 Para.vsize                : ボクセルサイズ
 Para.normal_subtract;
 Para.Nloop

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

0001 function    [V, F, xxn] = vb_surf_smooth_data(V,F,xx,vg,Para)
0002 % smoothing surface by spring model and 3D-image data
0003 %  [V, F, xxn] = vb_surf_smooth_data(V,F,xx,vg,Para)
0004 %
0005 % Ver 1.0  by M. Sato  2004-2-10
0006 %
0007 % ポリゴンモデルをマスクパターンの境界面に膨らませ
0008 % バネ力による平滑化を行う
0009 % 平滑化バネ力＋外向き力(内側)／内向き力(外側)
0010 %
0011 %  V(n, 1:3)  : 頂点の位置
0012 %  F(j, 1:3)  : 三角面３頂点のインデックス
0013 % xx(n, 1:3)  : 頂点の法線方向(nx,ny,nz)
0014 % xxn(n, 1:3) : 頂点の法線方向(nx,ny,nz)
0015 % vg(NX,NY,NZ): 3D-ボクセル・マスクパターン
0016 %
0017 % Para.tangent_rate         : 接線方向力係数
0018 % Para.tangent_normal_ratio : 法線方向力係数
0019 % Para.mask_ratio           : 画像強度力係数 ( 強度>閾値:外向き)
0020 % Para.mask_threshold       : 画像強度閾値   ( 強度<閾値:内向き)
0021 % Para.vsize                : ボクセルサイズ
0022 % Para.normal_subtract;
0023 % Para.Nloop
0024 %
0025 % Copyright (C) 2011, ATR All Rights Reserved.
0026 % License : New BSD License(see VBMEG_LICENSE.txt)
0027 
0028 c1       = Para.tangent_rate;
0029 c2       = Para.tangent_rate*Para.tangent_normal_ratio;
0030 c3       = Para.normal_subtract;
0031 c4       = Para.mask_ratio;
0032 g0       = Para.mask_threshold;
0033 
0034 Nloop  = Para.Nloop;
0035 Ndisp  = Para.Ndisp;
0036 cmode  = Para.cmode;
0037 vangle = Para.vangle;
0038 NX       = Para.NX;
0039 NY       = Para.NY;
0040 nfig   = 1;
0041 
0042 Npoint = size(V,1);          % number of dipoles
0043 Npatch = size(F,1);          % number of patch
0044 gsize  = size(vg);
0045 
0046 % ３角面の法線ベクトル
0047 xxf  = zeros(Npatch,3);
0048 % ３角面面積
0049 nns  = zeros(Npatch,1);
0050 % ３角面頂点の差分ベクトル
0051 VV1  = zeros(Npatch,3);
0052 VV2  = zeros(Npatch,3);
0053 VV3  = zeros(Npatch,3);
0054 
0055 % 頂点法線 = 頂点に隣接する三角面法線の平均
0056 xxn  = zeros(Npoint,3);    
0057 xxs  = zeros(Npoint,1);
0058 % 頂点差分ベクトル
0059 fd   = zeros(Npoint,3);
0060 % 各頂点の近傍点数
0061 Nv   = zeros(Npoint,1);
0062 % 頂点差分ベクトルと法線の内積
0063 nnf  = zeros(Npoint,1);
0064 % 接線方向の力
0065 ft     = zeros(Npoint,3);
0066 % 法線方向の力
0067 fn     = zeros(Npoint,3);
0068 % 外向き(内側)・内向き(外側)の力
0069 fg     = zeros(Npoint,3);
0070 
0071 % V 座標のボクセルインデックス
0072 jv     = zeros(Npoint,3);
0073 jj     = zeros(Npoint,1);
0074 % 内側マスク強度
0075 gg     = zeros(Npoint,1);
0076 
0077 % 三角面３頂点のインデックス
0078 F1     = F(:,1);
0079 F2     = F(:,2);
0080 F3     = F(:,3);
0081 
0082 % ３角面の法線ベクトル
0083 xxf  = vb_cross2( V(F2,:)-V(F1,:) , V(F3,:)-V(F1,:) );
0084 
0085 % 三角面３頂点の法線ベクトル平均
0086 xxk  = xx(F1,:)+xx(F2,:)+xx(F3,:);
0087 
0088 xdot = sum( xxk .* xxf ,2);
0089 ix     = find(sign(xdot) < 0 );
0090 
0091 % ３角面の法線ベクトルと
0092 % 頂点の法線ベクトル 'xx' (外向き)の向きをそろえる
0093 F(ix,2) = F3(ix);
0094 F(ix,3) = F2(ix);
0095 
0096 F2        = F(:,2);
0097 F3        = F(:,3);
0098 
0099 % 各頂点の近傍点数
0100 for n=1:Npatch,
0101     inx = F(n,:)';
0102     Nv(inx) = Nv(inx) + 2; 
0103 end;
0104 
0105 % バネ力平滑化＋外向き力(内側)＋内向き力(外側)
0106 
0107 for i=1:Nloop,
0108 
0109     % 三角面３頂点
0110     V1=V(F1,:);
0111     V2=V(F2,:);
0112     V3=V(F3,:);
0113     
0114     % ３角面の法線ベクトル
0115     xxf   = vb_cross2( V2 - V1, V3 - V1 );
0116 
0117     % Normalization
0118     nns   = sqrt(sum(xxf.^2,2));
0119     xxf   = xxf./nns(:,ones(1,3));        
0120     
0121     % 三角面各頂点の差分ベクトル
0122     VV1   = V2 + V3 - 2*V1;
0123     VV2   = V3 + V1 - 2*V2;
0124     VV3   = V1 + V2 - 2*V3;
0125 
0126     % 頂点法線 = 頂点に隣接する三角面法線の平均
0127     xxn   = zeros(Npoint,3);    
0128     % 頂点差分ベクトル(近傍和)
0129     fd      = zeros(Npoint,3);
0130     
0131     for n=1:Npatch,
0132         % 三角面３頂点インデックス
0133         j1=F1(n);
0134         j2=F2(n);
0135         j3=F3(n);
0136         
0137         xxn(j1,:) = xxn(j1,:) + xxf(n,:);
0138         xxn(j2,:) = xxn(j2,:) + xxf(n,:);
0139         xxn(j3,:) = xxn(j3,:) + xxf(n,:);
0140         
0141         % 頂点差分ベクトル(近傍和)
0142         fd(j1,:) = fd(j1,:) + VV1(n,:); 
0143         fd(j2,:) = fd(j2,:) + VV2(n,:); 
0144         fd(j3,:) = fd(j3,:) + VV3(n,:); 
0145     end;
0146     
0147     % 法線の正規化
0148     xxs   = sqrt(sum(xxn.^2,2));
0149     xxn   = xxn./xxs(:,ones(1,3));
0150 
0151     % 近傍点数で正規化
0152     fd     = fd./Nv(:,ones(1,3));
0153 
0154     % 頂点差分ベクトルと法線の内積
0155     nnf  = sum(xxn.*fd,2);
0156     
0157     % 頂点差分ベクトルの法線射影
0158     fn   = xxn.*nnf(:,ones(1,3));
0159     
0160     % 頂点差分ベクトルの接線射影
0161     % DEBUG
0162     ft  = fd - fn; 
0163     % DEBUG
0164     
0165     % 法線方向の平均力を引く
0166     fns = sum(fn,1)/Npoint;
0167     
0168     % V 座標をボクセルインデックスに変換
0169     jv    = floor(V)+1;
0170     jj    = sub2ind(gsize,jv(:,1),jv(:,2),jv(:,3));
0171 
0172     % 内側マスク強度
0173     gg    = tanh(vg(jj)-g0);
0174 
0175     fg    = xxn.*gg(:,ones(1,3));
0176     
0177     V    = V + c1*ft + c2*(fn - c3*fns(ones(Npoint,1),:)) + c4*fg;
0178 
0179     if rem(i,Ndisp)==0,
0180         subplot(NY,NX,nfig); nfig=nfig+1;
0181         vb_plot_head_V(V,F,cmode);
0182         view(vangle);
0183         tlabel = sprintf('Iteration = %d',i);
0184         title(tlabel);
0185         axis equal
0186     end;
0187 end
0188