n维空间的多面体的有向测度和重心
缘起
在《三维凸包》中我们学习了如何求三维空间中的点集凸包，本文来论述二维、三维甚至高位几何体的测度和重心的计算. 所谓测度，对于二维，指的是面积，对于三维，指的是体积. 所谓重心，指的是空间中一个特殊的点，如果该物体是质量分布均匀的话（所谓质量分布均匀，指的是密度函数是常数函数），则该物体关于该点力矩平衡.
分析
概念讲清楚了，下面学习一下计算方法. 我们的学习是循序渐进的——从二维到高维，但是不管对于二维也好，还是高维也罢，我们总是先计算测度，再计算重心. 为什么是这个顺序呢? 因为重心是力矩平衡的点，而因为我们已经假设质量分布均匀了. 所以测度的分布就决定了重心的位置. 所以我们先要会计算测度才能进一步计算重心. 测度在本文不做特别说明的话指的就是有向测度.
三角形的面积和重心
这个在之前的学习中早就知道了，三角形的有向面积使用叉积可以方便的计算出来.
则三角形的有向面积是
其中，
是 A 在平面的坐标, 下同.
当然对于上图的话，有向面积为正. 至于三角形的重心 G， 我们通过简单的平面几何相似便知道
这里说一下，单个点给出的话，则表示从坐标原点出发的向量. 例如我们写 A，实际上表达的是向量
, 下同.
简单论证如下, 我们知道 G 是 三条中线的交点.
所以 EF 是该三角形的中位线. 所以根据 HGE 和 BGD 相似，便知道了一切. 即 HG 是 GD 的长度的一半. 而 AH = HD,所以 GD 是 AD 的 1/3, HG 是 AD 的 1/6，AH 是 AD 的 1/2，所以便有了上述三角形的重心公式. 值得注意的是，三角形的重心和质心重合. 所谓质心的意思是将 A、B、C 视作放置了三个质量相等的质点之后，这三个质点组成的质心系的重心位置. 所以三角形好像是其内部的质量完全分布在其三个顶点上一样.
平面多边形的面积和重心
计算平面多边形的面积有如下十分优美的 O(n) 伪代码, 这里 n 是多边形的顶点个数,
是多边形的 n 个顶点. 下同
ild area()
{
double ans = 0;
for (re i = 0; i < n - 1; i++)
{
ans += ps[i] / ps[i + 1];
}
ans += ps[n - 1] / ps[0];
return fabs(ans) / 2;
}
上面伪代码的含义就是平面上任取一点（此点不一定在多边形内部），然后将平面多边形进行三角剖分, 然后平面多边形的面积就等于剖分出来的三角形的面积之和.
即多边形 ABCDEFG 的面积就等于 三角形 OAB、OBC、OCD、ODE、OEF、OFG、OGA 的面积之和（再次强调，是有向面积）. 而且哈，这里必须要提一嘴. 就是多边形的有向面积可以极为方便的定义多边形的正方向. 一言以蔽之，就是如果计算出的有向面积 >0 ,则规定给出的顶点的顺序就是多边形的正向，否则就是多边形的逆向. 而数学中如果要刻画曲线（确切讲是简单闭曲线，所谓简单，指的是曲线没有自相交），用的是 Jordan 曲线定理（Jordan curve theorem）
在欧氏平面上，任意一条简单闭曲线J把平面分成两部分，使得在同一部分的任意两点，可以用一条不与J相交的弧相连;
在不同部分的两点若要相连，则连结的弧必须与J相交
Jordan 曲线定理是属于典型的数学中那种看起来容易，证起来难的定理. 直到1905年，维布伦才第一次给出了一个正确的证明.
有了 Jordan 曲线定理，我们便有了两个区域，其中那个测度有限的区域就称之为简单闭曲线的内部，另一个称之为外部. 那么简单闭曲线的正向指的是，你沿着此方向在闭曲线上行走，内部在你左手边. 但是这种定义显然对计算机是不友好的. 而上述计算有向面积的定义对于计算机是极为有好的.
会计算多边形的面积了，再来考虑多边形的重心. 结合上面的三角形的重心计算，一个自然的猜测是
但是很遗憾，反例太多了. 最直观的反例就是梯形.
如果按照上面的公式的话，则计算出来的重心就是 G，其中 E、F 分别是 AD 、 BC 的中点. G 是 EF 的中点. 但是你觉得在 G 处支起一根筷子，梯形的纸板能平衡住吗? 显然不能的，因为 梯形 ADIH 显然没有 梯形 HICB 面积大啊, 也就是前者没有后者重啊~ 所以在 G 处支起一根筷子，梯形纸板是显然不能保持住平衡的. 事实上，直观上我们感觉真实的重心应该在 G 的下方. 那么为什么会导致这个错误呢? 因为对于四边形，乃至多边形，只要不是三角形，则重心和质心并不相同. 上面求出的 G 其实是质心，而不是重心. 为什么非三角形的多边形的质心和重心不重合，而三角形却能做到这一点呢? 这是因为三角形的特殊性——三角形不需要指定这三个顶点的顺序就能唯一确定一个三角形，多边形则不能. 就拿五边形 ABCDE 为例. 如果就给你5个点的话，你是无法确定该五边形长啥样的. 例如
既可以长左边这样，又可以长右边这样.
那么重心该怎么求呢?
正确的姿势应该是首先将n个顶点的多边形（可以凸，可以凹）剖分成 n 个三角形. 例如下图
即将剖分出来的每个三角形都抽象成其对应的重心，例如上图中三角形 AOB 就抽象成了 G2，BOC 就抽象成了 G4 等等. 于是问题就规约为了计算质心系 {G1, G2, G3, G4, G5} 的质心. 但是，这里注意，质心系中的每个质心的质量是不一样的. 因为三角形的面积不一样. 所以我们需要赋予这些质心以权重，相应的权就是三角形的面积（再次强调，是有向面积）. 即多边形的重心计算公式如下
其中 A 是多边形的总的有向面积（也即 n 个剖出来的三角形的有向面积之和），
是每个三角形的有向面积，根据上面的学习，我们知道
注意，为了图方便，我们已经将上图中的 O 取为了坐标原点 （0，0）. 所以
其中，
这里值得提一嘴的是，抛开公式本身，我们不难理解多边形的面积有正负之分，也就是和顶点给出的顺逆时针方向有关，但是多边形的重心应该和顶点给出的顺逆时针方向是无关的. 事实上，上面的公式也印证了我们的看法.
三维多面体的体积和重心
有了前面多边形的面积和重心的学习，我们立刻知道了，要考虑三维多面体的体积（确切讲，是有向体积）和重心，同样是三角剖分，当然，既然到了三维空间，所谓的三角 指的就是四面体，而非三角形了. 所以三维多面体的有向体积等于剖分出来的四面体有向体积之和，而三维多面体的重心等于各个四面体的重心关于四面体有向体积的加权平均.
所以首先，我们应该知道如何计算四面体的体积. 学过解析几何的话，这就太简单了，就是三维向量的混合积嘛~
即四面体 ABCD 的有向体积为
其中上式最后那个等号就是简单的做一下行列式的行变换就行了.
大家可以对比一下二维三角形的有向面积S和这里的三维四面体的有向体积V的计算公式，
四不四发现了墙裂的美感?
关于三维多面体的重心，我们将在下面一般的 n 维空间多面体的体积和重心中做出一般性的论述.
n 维空间多面体的体积和重心
显然，我们需要考虑 n 维空间的多面体对应的三角剖分. 这里就不得不提及数学中单纯形的概念. 单纯形是二维三角形和三维四面体的一种泛化，一个 n 维单纯形是指包含 n + 1 个顶点的凸多面体.
令 n 维单纯形的顶点坐标为
那么，n 维单纯形的有向测度为
有了单纯形的有向测度公式，我们就来推导单纯形的重心公式. 显然，我们只需要计算单纯形的重心，然后做一下关于有向测度的加权平均就行了. 同样，我们从二维三角形开始找规律.
三角形的面积为
于是重心的 y 坐标为
这个公式通过微元法+杠杆原理是很容易得到的. 而根据相似三角形
于是
这说明三角形的重心在下 1/3 分点处. 推而广之，对于n维单纯形，
, 所以
即重心总在下
分点处. 所以我们便知道了，n 维单纯形的重心坐标为
而要进一步得到三维多面体的重心，我们自然就需要考虑 3 维多面体的四面体剖分. 受平面多边形的三角剖分启发，可以选定空间中任意一点 O 作为所有四面体的一个顶点——当然，你可以选择 O 为坐标原点，这样的好处是 4 阶行列式蜕化为 3 阶的行列式. 但是这样的话，剖分出来的是底面为平面多边形（可能不是三角形）的多棱锥. 例如下图是五棱锥 O-ABCDE
所以要进一步将平面多边形（上图中的 ABCDE） 做三角剖分. 例如上图中选择A点作为平面上的点，将 ABCDE 剖分为 ABC、ADC、ADE 这 3 个三角形. 于是，我们可以得到如下有向体积和重心的公式
假设三维多面体有 m 个面，每个面都是一个多边形，多边形的顶点个数是
， 第
个面的顶点集合为
，显然，对于不同的 i 之间, 点集{
}之间是可能有重复点的. 这里我们选择 O 为 三维空间的坐标原点.
则该三维多面体的有向体积为
其中
其中
，所以三维多面体的有向体积为
重心的话，
其中
其中（注意，我们选择 O 为三维空间的坐标原点）
综上所述，三维多面体的重心公式如下
其中
是
的坐标，
是
的坐标，
是
的坐标
关于上面公式中的符号的说明：
是多面体第 i （
）个面和坐标原点 O 形成的多棱锥的重心，再剖分下去, 得到
是把第 i 个面剖分为
个三角形之后， 第 j （
）个三角形
和坐标原点 O 构成的四面体的重心.
至此，就彻底解决了三维多面体的有向体积和重心问题. 一般对于比赛，至此基本够用了. 然鹅让我们的思绪再发散一下，考虑一般 n 维空间中的多面体（可凸可凹）的有向测度和重心问题.
自然的，我们的思路依旧是做单纯形剖分. 即将多面体剖分为坐标原点为多棱锥顶点，多面体的某个面为底面的多棱锥.
此时，多面体的面是 n - 1 维空间（更一般的，是 n - 1 维流形），而每个 n - 1 维面又由若干 n - 2 维面构成，如此不断分解，直至二维面由若干顶点构成. 所以聪明如你早就知道刻画整个不断剖分的过程的最好的数据结构是多叉树. 这棵树的叶子节点是 点. 整棵树的深度是 n - 1 （例如上面的n = 3维空间就剖分了 2 次, 第一次剖分是多面体剖分成了多棱锥，这本质上来讲也就是体-->面，第二次剖分是多棱锥剖分成了三棱锥（也就是四面体），这本质上来讲，也就是面-->点）
上面的图是针对 n = 3 画的, 第一次剖分其实就是选择了 O（坐标原点），第二次剖分其实就是选择了每个多棱锥的底面上的一个顶点，例如 P10、P20、P30、P40， 这里假定该三维多面体有四个面（m=4）, 所以树的第二层节点有4个子节点，然后假定 P10 所在面是5边形，所以 P10 有 4个子节点，类似的，P30 所在面是四边形，所以 P30 有 3个子节点... 等等. 而计算该 n 维空间的多面体的有向测度的算法就是在此棵树上进行 dfs , 例如我们找到一条路径 O-->P20-->P22, 则三棱锥的有向测度是 P20、P22、P23 （P23 是 P22 的右兄弟）构成的 n (上图是 n = 3 的例子)阶行列式的值. 特别的，对于 O --> P20 --> P24, 则P24 的右兄弟应该是 P21哈~ 所以这个 n 阶行列式的值的几何意义是这 n 个顶点和树的根节点 O 形成的 n 维单纯形的有向测度. 然后dfs遍历完整棵树之后，累加得到的有向测度就得到了整个n维多面体的有向测度. 然后各个单纯形的重心关于有向测度的加权平均得到的就是整个多面体的重心.
而n维单纯形的重心是极为好求的. 就是该单纯形的 n + 1 个顶点的算术平均.
需要指出的是，n 阶行列式的计算，如果按照行列式的定义去计算，复杂度是
的，如果使用高斯消元的话，可以优化到
, 而因为叶子节点个数是 O(V) 的，其中 V 是多面体的顶点个数，则整个 n 维空间多面体的有向测度+重心的算法复杂度是
的.
至此，n维空间的多面体的有向测度+重心问题已经得到了圆满的解决.
有了上面的学习，我们来看看这道经典题目 hdu 4273 Rescue
给你三维空间中的 n 个点（4 <= n <= 100）的三维坐标 x, y, z, 其中 |x|,|y|,|z|<=10000
请你计算出这 n 个点张成的凸包的重心到凸包的各个面的距离的最小值. 请精确到小数点3位.
【输入】
多样例. 样例的第一行是 n, 然后 n 行, 每行包含 x y z 三个整数. 输入保证这 n 个点是一个非退化的凸多面体的顶点.
【输出】
答案
【样例输入】
4
0 0 0
1 0 0
0 1 0
0 0 1
8
0 0 0
0 0 2
0 2 0
0 2 2
2 0 0
2 0 2
2 2 0
2 2 2
【样例输出】
0.144
1.000
【限制】
1s 32MB
此题的思路是非常明显的，就是先维护出三维凸包, 然后计算重心，最后计算出答案了.
//#define LOCAL
#pragma GCC optimize(2)
#pragma G++ optimize(2)
#pragma warning(disable:4996)
#pragma comment(linker, "/STACK:1024000000,1024000000")
#include <stdio.h>
#include <iostream>
#include <iomanip>
#include <string>
#include <ctype.h>
#include <string.h>
#include <fstream>
#include <sstream>
#include <math.h>
#include <map>
//#include <unordered_map>
#include <algorithm>
#include <vector>
#include <stack>
#include <queue>
#include <set>
#include <time.h>
#include <stdlib.h>
#include <bitset>
using namespace std;
//#define int unsigned long long
//#define int long long
#define re register int
#define ci const int
#define ui unsigned int
typedef pair<int, int> P;
#define FE(cur) for(re h = head[cur], to; ~h; h = g[h].nxt)
#define ilv inline void
#define ili inline int
#define ilc inline char
#define ild inline double
#define ilp inline P
#define LEN(cur) (hjt[cur].r - hjt[cur].l)
#define MID(cur) (hjt[cur].l + hjt[cur].r >> 1)
#define SQUARE(x) ((x) * (x))
typedef vector<int>::iterator vit;
typedef set<int>::iterator sit;
typedef map<int, int>::iterator mit;
const int inf = ~0u >> 1;
//const int inf = 0x3f3f3f3f;
const double PI = acos(-1.0), eps = 1e-8;
namespace fastio
{
const int BUF = 1 << 21;
char fr[BUF], fw[BUF], * pr1 = fr, * pr2 = fr; int pw;
ilc gc() { return pr1 == pr2 && (pr2 = (pr1 = fr) + fread(fr, 1, BUF, stdin), *pr2 = 0, pr1 == pr2) ? EOF : *pr1++; }
ilv flush() { fwrite(fw, 1, pw, stdout); pw = 0; }
ilv pc(char c) { if (pw >= BUF) flush(); fw[pw++] = c; }
ili read(int& x)
{
x = 0; int f = 1; char c = gc(); if (!~c) return EOF;
while (!isdigit(c)) { if (c == '-') f = -1; c = gc(); }
while (isdigit(c)) x = (x << 3) + (x << 1) + (c ^ 48), c = gc();
x *= f; return 1;
}
ili read(double& x)
{
int xx = 0; double f = 1.0, fraction = 1.0; char c = gc(); if (!~c) return EOF;
while (!isdigit(c)) { if (c == '-') f = -1.0; c = gc(); }
while (isdigit(c)) { xx = (xx << 3) + (xx << 1) + (c ^ 48), c = gc(); }
x = xx;
if (c ^ '.') { x = f * xx; return 1; }
c = gc();
while (isdigit(c)) x += (c ^ 48) * (fraction /= 10), c = gc();
x *= f; return 1;
}
ilv write(int x) { if (x < 0) pc('-'), x = -x; if (x > 9) write(x / 10); pc(x % 10 + 48); }
ilv writeln(int x) { write(x); pc(10); }
ili read(char* x)
{
char c = gc(); if (!~c) return EOF;
while (!isalpha(c) && !isdigit(c)) c = gc();
while (isalpha(c) || isdigit(c)) *x++ = c, c = gc();
*x = 0; return 1;
}
ili readln(char* x)
{
char c = gc(); if (!~c) return EOF;
while (c == 10) c = gc();
while (c >= 32 && c <= 126) *x++ = c, c = gc();
*x = 0; return 1;
}
ilv write(char* x) { while (*x) pc(*x++); }
ilv write(const char* x) { while (*x) pc(*x++); }
ilv writeln(char* x) { write(x); pc(10); }
ilv writeln(const char* x) { write(x); pc(10); }
ilv write(char c) { pc(c); }
ilv writeln(char c) { write(c); pc(10); }
} using namespace fastio;
#define cp const Point
#define cs const Surface
ci maxn = 105;
ili dbcmp(double x)
{
if (fabs(x) < eps)
{
return 0;
}
return x > eps ? 1 : -1;
}
ilv mmin(double& a, double b)
{
if (a > b)
{
a = b;
}
}
struct CH3D
{
struct Point
{
double x, y, z;
Point(double x = 0, double y = 0, double z = 0) :x(x), y(y), z(z) {}
Point operator + (cp& o) const
{
return Point(x + o.x, y + o.y, z + o.z);
}
Point operator - (cp& o) const
{
return Point(x - o.x, y - o.y, z - o.z);
}
double operator * (cp& o) const
{
return x * o.x + y * o.y + z * o.z;
}
Point operator / (cp& o) const
{
return Point(y * o.z - z * o.y, z * o.x - x * o.z, x * o.y - y * o.x);
}
double magnitude() const
{
return sqrt(SQUARE(x) + SQUARE(y) + SQUARE(z));
}
Point scalar(double lambda) const
{
return Point(x * lambda, y * lambda, z * lambda);
}
} ps[maxn];
int n;
struct Surface
{
int a, b, c, flag;
Surface(int a = 0, int b = 0, int c = 0) :a(a), b(b), c(c) { flag = 1; }
} surfaces[maxn << 3];
int num;
int g[maxn][maxn];
ild area(cp& a, cp& b, cp& c)
{
return ((b - a) / (c - a)).magnitude() / 2.0;
}
ild area(cs& sur)
{
return area(ps[sur.a], ps[sur.b], ps[sur.c]);
}
ild area()
{
double ans = 0;
if (n == 3)
{
return area(ps[0], ps[1], ps[2]);
}
for (re i = 0; i < num; i++)
{
ans += area(ps[surfaces[i].a], ps[surfaces[i].b], ps[surfaces[i].c]);
}
return ans;
}
ild volume(cp& a, cp& b, cp& c, cp& d)
{
return (((b - a) / (c - a)) * (d - a)) / 6.0;
}
ild volume(cp& p, cs& sur)
{
return volume(ps[sur.a], ps[sur.b], ps[sur.c], p);
}
ild volume()
{
double ans = 0;
Point o;
for (re i = 0; i < num; i++)
{
ans += volume(o, ps[surfaces[i].a], ps[surfaces[i].b], ps[surfaces[i].c]);
}
return ans;
}
ili ck(cp& p, cs& sur)
{
return dbcmp(volume(ps[sur.a], ps[sur.b], ps[sur.c], p)) == 1;
}
ili onsurface(cp& p, cs& sur)
{
return !dbcmp(volume(ps[sur.a], ps[sur.b], ps[sur.c], p));
}
ili same(cs& s, cs& t)
{
return onsurface(ps[s.a], t) && onsurface(ps[s.b], t) && onsurface(ps[s.c], t);
}
ili prework()
{
if (n < 4)
{
return 0;
}
int ok = 1;
for (re i = 1; i < n; i++)
{
if (dbcmp((ps[0] - ps[i]).magnitude()))
{
swap(ps[1], ps[i]);
ok = 0;
break;
}
}
if (ok)
{
return 0;
}
ok = 1;
for (re i = 2; i < n; i++)
{
if (dbcmp(area(ps[0], ps[1], ps[i])))
{
swap(ps[2], ps[i]);
ok = 0;
break;
}
}
if (ok)
{
return 0;
}
ok = 1;
for (re i = 3; i < n; i++)
{
if (dbcmp(volume(ps[0], ps[1], ps[2], ps[i])))
{
swap(ps[3], ps[i]);
ok = 0;
break;
}
}
return !ok;
}
ilv construct()
{
num = 0;
if (!prework())
{
return;
}
Surface add;
for (re i = 0; i < 4; i++)
{
Surface add((i + 1) % 4, (i + 2) % 4, (i + 3) % 4);
if (ck(ps[i], add))
{
swap(add.b, add.c);
}
g[add.a][add.b] = g[add.b][add.c] = g[add.c][add.a] = num;
surfaces[num++] = add;
}
for (re i = 4; i < n; i++)
{
for (re j = 0; j < num; j++)
{
if (surfaces[j].flag && ck(ps[i], surfaces[j]))
{
dfs(i, j);
break;
}
}
}
int tmp = num; num = 0;
for (re i = 0; i < tmp; i++)
{
if (surfaces[i].flag)
{
surfaces[num++] = surfaces[i];
}
}
}
void dfs(int p, int j)
{
surfaces[j].flag = 0;
rmv(p, surfaces[j].b, surfaces[j].a);
rmv(p, surfaces[j].c, surfaces[j].b);
rmv(p, surfaces[j].a, surfaces[j].c);
}
void rmv(int p, int a, int b)
{
int f = g[a][b];
if (surfaces[f].flag)
{
if (ck(ps[p], surfaces[f]))
{
dfs(p, f);
}
else
{
Surface add(b, a, p);
g[b][a] = g[a][p] = g[p][b] = num;
surfaces[num++] = add;
}
}
}
ild dis(cp& p, cs& sur)
{
return fabs(volume(p, sur)) / area(sur) * 3.0;
}
ild dis(cp& p)
{
double ans = inf;
for (int i = 0; i < num; i++)
{
mmin(ans, dis(p, surfaces[i]));
}
return ans;
}
Point gg()
{
Point ans, o = ps[0];
double v = 0, t;
for (re i = 0; i < num; i++)
{
cp& a = ps[surfaces[i].a];
cp& b = ps[surfaces[i].b];
cp& c = ps[surfaces[i].c];
t = volume(o, surfaces[i]);
ans = ans + (a + b + c + o).scalar(t / 4);
v += t;
}
return ans.scalar(1 / v);
}
}ch;
signed main()
{
#ifdef LOCAL
FILE* ALLIN = freopen("d:\\data.in", "r", stdin);
freopen("d:\\my.out", "w", stdout);
#endif
while (~scanf("%d", &ch.n))
// while (~read(ch.n))
{
// for (re i = 0; i < ch.n; i++) read(ch.ps[i].x), read(ch.ps[i].y), read(ch.ps[i].z);
for (re i = 0; i < ch.n; i++) scanf("%lf%lf%lf", &ch.ps[i].x, &ch.ps[i].y, &ch.ps[i].z);
ch.construct();
CH3D::Point centroid = ch.gg();
printf("%.3lf\n", ch.dis(centroid));
}
flush();
#ifdef LOCAL
fclose(ALLIN);
#endif
return 0;
}
ac情况
Accepted 31ms 1804kB C++
值得提醒一下的是，本题貌似用快读就会 T.