vb_surf_smooth

 smoothing surface by spring model
  [V_out, F_out, xx_out] = vb_surf_smooth(V, F, xx, Para)

  V(n, 1:3)  : vertex on the surface
  F(j, 1:3)  : patch index for surface
 xx(n, 1:3)  : normal vector at each vertex
 
 Para.Nloop  : loop iteration number
 Para.tangent_rate : coefficient of spring force
 Para.normal_mode = 0: force = spring force
 Para.normal_mode = 1: force = spring force + compensatory expansion force 
 --------------------------------------------------
 ポリゴンモデルをバネ力による平滑化を行う

  V(n, 1:3)  : 頂点の位置
  F(j, 1:3)  : 三角面３頂点のインデックス
 xx(n, 1:3)  : 頂点の法線方向(nx,ny,nz)
 
 Para.Nloop  : 繰り返し回数
 Para.tangent_rate : 接線方向力係数
 Para.normal_mode = 0: 平滑化バネ力
 Para.normal_mode = 1: 平滑化バネ力 + 平均半径を保つ法線方向力 (default)

 平滑化バネ力 =  tangent_rate * (近傍点の平均座標との差)

 Ver 1.0  by M. Sato  2004-2-10
 modefied M. Sato 2006-10-15
 modefied M. Sato 2010-10-3

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

0001 function    [V, F, xxn] = vb_surf_smooth(V,F,xx,Para)
0002 % smoothing surface by spring model
0003 %  [V_out, F_out, xx_out] = vb_surf_smooth(V, F, xx, Para)
0004 %
0005 %  V(n, 1:3)  : vertex on the surface
0006 %  F(j, 1:3)  : patch index for surface
0007 % xx(n, 1:3)  : normal vector at each vertex
0008 %
0009 % Para.Nloop  : loop iteration number
0010 % Para.tangent_rate : coefficient of spring force
0011 % Para.normal_mode = 0: force = spring force
0012 % Para.normal_mode = 1: force = spring force + compensatory expansion force
0013 % --------------------------------------------------
0014 % ポリゴンモデルをバネ力による平滑化を行う
0015 %
0016 %  V(n, 1:3)  : 頂点の位置
0017 %  F(j, 1:3)  : 三角面３頂点のインデックス
0018 % xx(n, 1:3)  : 頂点の法線方向(nx,ny,nz)
0019 %
0020 % Para.Nloop  : 繰り返し回数
0021 % Para.tangent_rate : 接線方向力係数
0022 % Para.normal_mode = 0: 平滑化バネ力
0023 % Para.normal_mode = 1: 平滑化バネ力 + 平均半径を保つ法線方向力 (default)
0024 %
0025 % 平滑化バネ力 =  tangent_rate * (近傍点の平均座標との差)
0026 %
0027 % Ver 1.0  by M. Sato  2004-2-10
0028 % modefied M. Sato 2006-10-15
0029 % modefied M. Sato 2010-10-3
0030 %
0031 % Copyright (C) 2011, ATR All Rights Reserved.
0032 % License : New BSD License(see VBMEG_LICENSE.txt)
0033 
0034 % normal_mode = 0: 平滑化バネ力
0035 % normal_mode = 1: 平滑化バネ力 + 平均半径を保つ法線方向力
0036 normal_mode = 1;
0037 
0038 if ~exist('Para','var'), Para = []; end
0039 if ~isfield(Para,'Nloop'), Para.Nloop = 50; end
0040 if ~isfield(Para,'tangent_rate'), Para.tangent_rate = 0.3; end
0041 if isfield(Para,'normal_mode'), normal_mode = Para.normal_mode; end;
0042 
0043 tangent_rate = Para.tangent_rate;
0044 
0045 Nloop  = Para.Nloop;
0046 
0047 Npoint = size(V,1);          % number of dipoles
0048 Npatch = size(F,1);          % number of patch
0049 
0050 % ３角面頂点ベクトル
0051 V1  = zeros(Npatch,3);
0052 V2  = zeros(Npatch,3);
0053 V3  = zeros(Npatch,3);
0054 % ３角面頂点の差分ベクトル
0055 VV1 = zeros(Npatch,3);
0056 VV2 = zeros(Npatch,3);
0057 VV3 = zeros(Npatch,3);
0058 
0059 % 頂点差分ベクトル
0060 fd   = zeros(Npoint,3);
0061 % 各頂点の近傍点数
0062 Nv   = zeros(Npoint,1);
0063 % 頂点差分ベクトルと法線の内積
0064 nnf  = zeros(Npoint,1);
0065 
0066 % 三角面３頂点のインデックス
0067 F1     = F(:,1);
0068 F2     = F(:,2);
0069 F3     = F(:,3);
0070 
0071 % ３角面の法線ベクトル
0072 xxf  = vb_cross2( V(F2,:)-V(F1,:) , V(F3,:)-V(F1,:) );
0073 
0074 % 三角面３頂点の法線ベクトル平均
0075 xxk  = xx(F1,:)+xx(F2,:)+xx(F3,:);
0076 
0077 xdot = sum( xxk .* xxf ,2);
0078 ix     = find(sign(xdot) < 0 );
0079 
0080 % ３角面の法線ベクトルと
0081 % 頂点の法線ベクトル 'xx' (外向き)の向きをそろえる
0082 F(ix,2) = F3(ix);
0083 F(ix,3) = F2(ix);
0084 
0085 F2        = F(:,2);
0086 F3        = F(:,3);
0087 
0088 % 各頂点の近傍点数
0089 for n=1:Npatch,
0090     inx = F(n,:)';
0091     Nv(inx) = Nv(inx) + 2; 
0092 end;
0093 
0094 % 平滑化バネ力＋外向き力
0095 
0096 for i=1:Nloop,
0097 
0098     % 三角面３頂点
0099     V1=V(F1,:);
0100     V2=V(F2,:);
0101     V3=V(F3,:);
0102     
0103     % ３角面の法線ベクトル
0104     xxf   = vb_cross2( V2 - V1, V3 - V1 );
0105 
0106     % Normalization
0107     xxf   = vb_repmultiply(xxf, 1./sqrt(sum(xxf.^2,2)));
0108     
0109     % 三角面各頂点の差分ベクトル
0110     VV1   = V2 + V3 - 2*V1;
0111     VV2   = V3 + V1 - 2*V2;
0112     VV3   = V1 + V2 - 2*V3;
0113 
0114     % 頂点法線 = 頂点に隣接する三角面法線の平均
0115     xxn   = zeros(Npoint,3);    
0116     % 頂点差分ベクトル(近傍和)
0117     fd      = zeros(Npoint,3);
0118     
0119     for n=1:Npatch,
0120         % 三角面３頂点インデックス
0121         j1=F1(n);
0122         j2=F2(n);
0123         j3=F3(n);
0124         
0125         % ３角面の法線ベクトルの和
0126         xxn(j1,:) = xxn(j1,:) + xxf(n,:);
0127         xxn(j2,:) = xxn(j2,:) + xxf(n,:);
0128         xxn(j3,:) = xxn(j3,:) + xxf(n,:);
0129         
0130         % 頂点差分ベクトル(近傍和)
0131         fd(j1,:) = fd(j1,:) + VV1(n,:); 
0132         fd(j2,:) = fd(j2,:) + VV2(n,:); 
0133         fd(j3,:) = fd(j3,:) + VV3(n,:); 
0134     end;
0135 
0136     % 法線の正規化
0137     xxn = vb_repmultiply(xxn, 1./sqrt(sum(xxn.^2,2)));
0138 
0139     % 近傍点数で正規化
0140     fd  = vb_repmultiply(fd, 1./Nv);
0141     
0142     % 各点の変化ベクトルの法線方向射影平均
0143     ds  = mean( sum(fd .* xxn, 2) );
0144     % dd  = mean( sqrt(sum(fd .^2, 2)) ) % mean displacement
0145     
0146     if normal_mode==1
0147         % 平滑化バネ力 + 平均半径を保つ法線方向力
0148         V  = V + tangent_rate * (fd - ds*xxn);
0149     else
0150         % 平滑化バネ力
0151         V  = V + tangent_rate * fd ;
0152     end
0153 end
0154 
0155 if nargout < 3, return; end;
0156 
0157 % 三角面３頂点
0158 V1=V(F1,:);
0159 V2=V(F2,:);
0160 V3=V(F3,:);
0161 
0162 % ３角面の法線ベクトル
0163 xxf   = vb_cross2( V2 - V1, V3 - V1 );
0164 
0165 % Normalization
0166 xxf   = vb_repmultiply(xxf, 1./sqrt(sum(xxf.^2,2)));
0167 
0168 % 頂点法線 = 頂点に隣接する三角面法線の平均
0169 xxn   = zeros(Npoint,3);    
0170 
0171 for n=1:Npatch,
0172     % 三角面３頂点インデックス
0173     j1=F1(n);
0174     j2=F2(n);
0175     j3=F3(n);
0176     
0177     xxn(j1,:) = xxn(j1,:) + xxf(n,:);
0178     xxn(j2,:) = xxn(j2,:) + xxf(n,:);
0179     xxn(j3,:) = xxn(j3,:) + xxf(n,:);
0180 end;
0181 
0182 % 法線の正規化
0183 xxn = vb_repmultiply(xxn, 1./sqrt(sum(xxn.^2,2)));