用高斯定理求多面体的质心

问题

任意给定一个未必凸的、均质多面体的顶点坐标、顶点和面的局部拓扑关系，求的多面体的质心坐标的算法如何设计？

多面体有关的问题之所以重要，看看图就明白了，这也是多面体及其三角网格：

竖起来是超人,躺下去就是金缕玉衣:

多面体数据结构的约定

算法总要基于特定的数据结构。关键应该是多面体以何种形式给出的问题。算法不纠结于多面体本身，以及对多面体进行三角形网格化的算法的问题。所以，一个多面体的数据结构是这样的：

假设多面体的所有面都已经化作三角形网格。这个实现起来应该不难，如果没有细化的要求，随便把多面体的某个面分解成一个（凸时）或多个凸多边形或至多一个顶点内凹的凹多边形（非凸时），然后凸多边形任取一个顶点与不相邻顶点连线，凹多边形则取凹顶点与不相邻顶点连线，与每条边的两个端点组成一个三角形;

需要细化的情况下,直接用Delaunay triangulation吧还是，这样，就有了顶点，以及由三个两两相连顶点组成的面（三角形，顶点相连即局部拓扑关系）。数据存储方式是这样的：

• 顶点，$V_{n\times 3}$ 的矩阵, 每一行保存一个顶点的三维欧式坐标;
• 面，$F_{m\times 3}$ 的矩阵, 每一行对应于一个面，按照逆时针顺序保存着组成这个面的三个点在顶点矩阵 $V$ 中的序号
  
  上图中6个顶点4个面的情况下,数据是这么存储的:
  $$V_{6\times 3}=\left( \begin{array}{ccc} x_0 & y_0 & z_0 \\ x_1 & y_1 & z_1 \\ x_2 & y_2 & z_2 \\ x_3 & y_3 & z_3 \\ x_4 & y_4 & z_4 \\ x_5 & y_5 & z_5 \\ \end{array}\right),$$
  $$F_{4\times 3}=\left( \begin{array}{ccc} 0 & 1 & 5 \\ 1 & 2 & 5 \\ 2 & 3 & 5 \\ 3 & 4 & 5 \\ \end{array} \right)$$

后面的Matlab代码也是如此，不过代码引用中出现 F(:,1,:)似乎是有问题的. 在Matlab中会忽略这个问题,但是这样写多余而且在对 $V$ 和 $F$ 数据结构不知情的情况下容易引起误解。

数学原理

算法的类型:

一个简单的算法的原理

算法的数学原理

Divergence theorem
${\bf F}=\langle X,Y,Z\rangle$ 是向量场, 则在光滑有向曲面 $S$ 为边界所包围的简单连通的三维封闭区域 $D$ 上

$$\mathop{{\int\!\!\!\!\!\int}\mkern-24mu \bigcirc}\limits_S \;{\bf{F}}\cdot {\bf n}\;{\rm d}S=\iiint\limits_D \nabla \cdot {\bf F}\;{\rm d}V$$
三维实心物体质心坐标的确定方法 $\rho(x,y,z)\equiv 1$
$${\bar x}=\dfrac{M_{yz}}{m}=\dfrac{1}{m}\iiint\limits_D x\cdot \rho(x,y,z){\rm d}V= \dfrac{1}{m}\iiint\limits_D x\;{\rm d}V$$
$${\bar y}=\dfrac{M_{xz}}{m}=\dfrac{1}{m}\iiint\limits_D y\cdot \rho(x,y,z){\rm d}V= \dfrac{1}{m}\iiint\limits_D y\;{\rm d}V$$
$${\bar z}=\dfrac{M_{xy}}{m}=\dfrac{1}{m}\iiint\limits_D z\cdot \rho(x,y,z){\rm d}V= \dfrac{1}{m}\iiint\limits_D z\;{\rm d}V$$

上面这些还算容易理解。为了求质心坐标，先要计算多面体体积，找一个满足 $\nabla\cdot {\bf F}=1$ 的向量场 $\bf F$ 即可， $\langle x,0,0\rangle$, $\langle 0,y,0\rangle$, $\langle 0,0,z\rangle$, $\tfrac{1}{3}\langle x,y,z\rangle$ 都行，取最后一个。从而：

$$m=V=\iiint\limits_V {\rm d}V=\dfrac{1}{3}\mathop{{\int\!\!\!\!\!\int}\mkern-24mu \bigcirc}\limits_S\;\left(\langle x,y,z\rangle\cdot {\bf n}\right)\;{\rm d}S=\dfrac{1}{6}\sum\limits_{A_i}{\bf a}_i\;{\bf n}_i=\dfrac{1}{6}\sum\limits_{A_i}{\bf {\bf b}}_i\;{\bf n}_i$$
其中 ${\bf n}_i$ 是多面体第 $i$ 个三角形面 $A_i \left({\bf a}_i,{\bf b}_i,{\bf c}_i\right)$ 法向量 $\left({\bf b}_i-{\bf a}_i\right)\otimes \left({\bf c}_i-{\bf b}_i\right)$，它的模是小三角形面积的 $2$ 倍，——这是$\dfrac{1}{3}$ 变 $\dfrac{1}{6}$ 的原因；实际上不但三个顶点 ${\bf a}_i,{\bf b}_i,{\bf c}_i$ ，小三角形所在平面上所有点对应的向量，与这个法向量的内积都是常数 。 到这里都还算容易理解，而且每一步都是精确相等，没有近似计算的必要。

算法的Matlab实现:

[假设已经对面作了三级剖分] :

Matlab代码:

Stackoverflow1:

Stackoverflow2

代码实现

http://www.alecjacobson.com/weblog/?p=3854

I ran into a concise algorithm for compute the centroid (center of mass) of an arbitrary closed 3D polyhedron. If you have a closed triangle mesh with vertices $V$ and faces in $F$ then here’s an implementation in matlab:

下面的算法基于这样的数据结构和前提： 封闭的多面体表面已经被划分为三角形网格，网格的顶点存贮在矩阵 $V$ 中，用来表达多面体三角形网格中每个小三角面的面与顶点之间的局部拓扑关系则保存在矩阵 $F$ 中。

从代码看毛病， 前面用 F(:,1,:), 后面则 F(:,1)说明前者根本就是多余的；

% Rename corners
A = V(F(:,1,:),:);
B = V(F(:,2,:),:);
C = V(F(:,3,:),:);
% Needs to be **unnormalized** normals
N = cross(B-A,C-A,2);
% total volume via divergence theorem: ∫ 1
vol = sum(sum(V(F(:,1),:).*N))/6;
% centroid via divergence theorem and midpoint quadrature: ∫ x
C = 1/(2*vol)*(1/24* sum(N.*((A+B).^2 + (B+C).^2 + (C+A).^2)));

这里的文件centroid.m在同一作者发布的Matlab工具箱 gptoolbox 中包含(其中的稳健算法果然把卷绕数(广义)用上了):

function [C,vol] = centroid(V,F,varargin)
  % CENTROID Compute the centroid of a closed polyhedron boudned by (V,F)
  %
  % C = centroid(V,F)
  % [C,vol] = centroid(V,F,'ParameterName',ParameterValue, ...)
  %
  % Inputs:
  %   V  #V by 3 list of mesh vertex positions
  %   F  #F by 3 list of triangle mesh indices
  %   Optional:
  %     'Robust' followed by whether to use more robust but costlier method for
  %       nearly closed input. {false}
  % Outputs:
  %   C  3-vector of centroid location
  %   vol  total volume of polyhedron
  %

  % default values
  robust = false;
  % Map of parameter names to variable names
  params_to_variables = containers.Map( {'Robust'}, {'robust'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace 
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  if robust
    [TV,TT] = cdt(V,F);
    BC = barycenter(TV,TT);
    w = winding_number(V,F,barycenter(TV,TT))/(4*pi);
    TT = TT(abs(w)>0.5,:);
    BC = BC(abs(w)>0.5,:);
    vol = volume(TV,TT);
    C = sum(bsxfun(@times,vol,BC))/sum(vol);
    vol = sum(vol);
  else
    % "Calculating the volume and centroid of a polyhedron in 3d" [Nuernberg 2013]
    % [http://www2.imperial.ac.uk/~rn/centroid.pdf](http://www2.imperial.ac.uk/~rn/centroid.pdf)

    % Rename corners
    A = V(F(:,1,:),:);
    B = V(F(:,2,:),:);
    C = V(F(:,3,:),:);
    % Needs to be **unnormalized** normals
    N = cross(B-A,C-A,2);
    % total volume via divergence theorem: ∫ 1
    vol = sum(sum(A.*N))/6;
    % centroid via divergence theorem and midpoint quadrature: ∫ x
    C = 1/(2*vol)*(1/24* sum(N.*((A+B).^2 + (B+C).^2 + (C+A).^2)));
  end
end

这是打个草稿,更新中
更新着更新着就跑题了…