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

最新推荐文章于 2024-02-03 15:27:49 发布

guankeliang

最新推荐文章于 2024-02-03 15:27:49 发布

阅读量1.7k

点赞数 1

分类专栏：数学

数学专栏收录该内容

1 篇文章 0 订阅

订阅专栏

问题

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

多面体有关的问题之所以重要，看看图就明白了，这也是多面体及其三角网格：
这里写图片描述
竖起来是超人,躺下去就是金缕玉衣:

多面体数据结构的约定

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

假设多面体的所有面都已经化作三角形网格。这个实现起来应该不难，如果没有细化的要求，随便把多面体的某个面分解成一个（凸时）或多个凸多边形或至多一个顶点内凹的凹多边形（非凸时），然后凸多边形任取一个顶点与不相邻顶点连线，凹多边形则取凹顶点与不相邻顶点连线，与每条边的两个端点组成一个三角形;
这里写图片描述
需要细化的情况下,直接用Delaunay triangulation吧还是，这样，就有了顶点，以及由三个两两相连顶点组成的面（三角形，顶点相连即局部拓扑关系）。数据存储方式是这样的：

顶点， Vn×3 的矩阵, 每一行保存一个顶点的三维欧式坐标;
面， Fm×3 的矩阵, 每一行对应于一个面，按照逆时针顺序保存着组成这个面的三个点在顶点矩阵 V 中的序号

上图中6个顶点4个面的情况下,数据是这么存储的:
$V 6 \times 3 = ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ x 0 x 1 x 2 x 3 x 4 x 5 y 0 y 1 y 2 y 3 y 4 y 5 z 0 z 1 z 2 z 3 z 4 z 5 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟,$
$F 4 \times 3 = ⎛ ⎝ ⎜ ⎜ ⎜ 012312345555 ⎞ ⎠ ⎟ ⎟ ⎟$

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

数学原理

算法的类型:

一个简单的算法的原理

算法的数学原理：http://www2.imperial.ac.uk/~rn/centroid.pdf

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

\int \int ◯ S F \cdot n d S = ∭ D \nabla \cdot F d V

三维实心物体质心坐标的确定方法

ρ(x,y,z)≡1

x ¯ = M y z m = 1 m ∭ D x \cdot ρ (x, y, z) d V = 1 m ∭ D x d V

y ¯ = M x z m = 1 m ∭ D y \cdot ρ (x, y, z) d V = 1 m ∭ D y d V

z ¯ = M x y m = 1 m ∭ D z \cdot ρ (x, y, z) d V = 1 m ∭ D z d V

上面这些还算容易理解。为了求质心坐标，先要计算多面体体积，找一个满足 ∇⋅F=1 的向量场 F 即可， ⟨x,0,0⟩ , ⟨0,y,0⟩ , ⟨0,0,z⟩ , 13⟨x,y,z⟩ 都行，取最后一个。从而：

m = V = ∭ V d V = 1 3 \int \int ◯ S (⟨ x, y, z ⟩ \cdot n) d S = 1 6 \sum A i a i n i = 1 6 \sum A i b i n i

其中

ni 是多面体第

i 个三角形面

Ai(ai,bi,ci) 法向量

(bi−ai)⊗(ci−bi) ，它的模是小三角形面积的

2 倍，——这是

13 变

16 的原因；实际上不但三个顶点

ai,bi,ci ，小三角形所在平面上所有点对应的向量，与这个法向量的内积都是常数。 到这里都还算容易理解，而且每一步都是精确相等，没有近似计算的必要。

算法的Matlab实现:

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

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

Stackoverflow1:http://stackoverflow.com/questions/9325303/centroid-of-convex-polyhedron/

Stackoverflow2：http://stackoverflow.com/questions/14186837/calculation-of-centroid-volume-of-a-polyhedron-when-the-vertices-are-given

代码实现

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)说明前者根本就是多余的；

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

  <span class="hljs-comment">% default values</span>
  robust = false;
  <span class="hljs-comment">% Map of parameter names to variable names</span>
  params_to_variables = <span class="hljs-transposed_variable">containers.</span>Map( <span class="hljs-cell">{<span class="hljs-string">'Robust'</span>}</span>, <span class="hljs-cell">{<span class="hljs-string">'robust'</span>}</span>);
  v = <span class="hljs-number">1</span>;
  <span class="hljs-keyword">while</span> v <= <span class="hljs-built_in">numel</span>(varargin)
    param_name = varargin<span class="hljs-cell">{v}</span>;
    <span class="hljs-keyword">if</span> isKey(params_to_variables,param_name)
      assert(v+<span class="hljs-number">1</span><=<span class="hljs-built_in">numel</span>(varargin));
      v = v+<span class="hljs-number">1</span>;
      <span class="hljs-comment">% Trick: use feval on anonymous function to use assignin to this workspace </span>
      feval(@()assignin(<span class="hljs-string">'caller'</span>,params_to_variables(param_name),varargin<span class="hljs-cell">{v}</span>));
    <span class="hljs-keyword">else</span>
      error(<span class="hljs-string">'Unsupported parameter: %s'</span>,varargin<span class="hljs-cell">{v}</span>);
    <span class="hljs-keyword">end</span>
    v=v+<span class="hljs-number">1</span>;
  <span class="hljs-keyword">end</span>

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

    <span class="hljs-comment">% Rename corners</span>
    A = V(F(:,<span class="hljs-number">1</span>,:),:);
    B = V(F(:,<span class="hljs-number">2</span>,:),:);
    C = V(F(:,<span class="hljs-number">3</span>,:),:);
    <span class="hljs-comment">% Needs to be **unnormalized** normals</span>
    N = <span class="hljs-built_in">cross</span>(B-A,C-A,<span class="hljs-number">2</span>);
    <span class="hljs-comment">% total volume via divergence theorem: ∫ 1</span>
    vol = sum(sum(<span class="hljs-transposed_variable">A.</span>*N))/<span class="hljs-number">6</span>;
    <span class="hljs-comment">% centroid via divergence theorem and midpoint quadrature: ∫ x</span>
    C = <span class="hljs-number">1</span>/(<span class="hljs-number">2</span>*vol)*(<span class="hljs-number">1</span>/<span class="hljs-number">24</span>* sum(<span class="hljs-transposed_variable">N.</span>*((A+B).^<span class="hljs-number">2</span> + (B+C).^<span class="hljs-number">2</span> + (C+A).^<span class="hljs-number">2</span>)));
  <span class="hljs-keyword">end</span>
<span class="hljs-keyword">end</span></code>