vb_solid_angle_check

 total solid angle for surface
   omega  = vb_solid_angle_check(V,F)
 omega : 重心から見た三角面の立体角総和 / (4*pi)
         閉局面に対しては omega = 1 となるべき
 V     : 境界三角面の頂点座標
 F     : 境界三角面頂点インデックス

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

0001 function omega  = vb_solid_angle_check(V,F)
0002 % total solid angle for surface
0003 %   omega  = vb_solid_angle_check(V,F)
0004 % omega : 重心から見た三角面の立体角総和 / (4*pi)
0005 %         閉局面に対しては omega = 1 となるべき
0006 % V     : 境界三角面の頂点座標
0007 % F     : 境界三角面頂点インデックス
0008 %
0009 % Copyright (C) 2011, ATR All Rights Reserved.
0010 % License : New BSD License(see VBMEG_LICENSE.txt)
0011 
0012 % 三角面の数　Number of triangle patch
0013 NF    = size(F,1);
0014 NV    = size(V,1);
0015 
0016 % 三角面の頂点　triangle vertex
0017 x1    = V(F(:,1),:);
0018 x2    = V(F(:,2),:);
0019 x3    = V(F(:,3),:);
0020 
0021 % 重心座標
0022 x0    = sum(V,1)/NV;
0023 xn    = x0( ones(NF,1) ,:);
0024 
0025 % xn から見た三角面頂点ベクトル　triangle vector
0026 xa    = x1 - xn;
0027 xb    = x2 - xn;
0028 xc    = x3 - xn;
0029 
0030 % 立体角計算
0031 abc    = sum(xa.*vb_cross2(xb,xc), 2);
0032 
0033 ra    = sqrt(sum(xa.^2,2));
0034 rb    = sqrt(sum(xb.^2,2));
0035 rc    = sqrt(sum(xc.^2,2));
0036 
0037 ya    = sum(xb.*xc,2);
0038 yb    = sum(xc.*xa,2);
0039 yc    = sum(xa.*xb,2);
0040 
0041 % 立体角総和
0042 omega = sum(atan(abc./(ra.*rb.*rc + ra.*ya + rb.*yb + rc.*yc)));
0043 omega = omega/(2*pi);