「SDOI2016」平凡的骰子 - 计算几何+三维重心+二面角

题目大意

    这是一枚平凡的骰子。它是一个均质凸多面体，表面有$n$个端点，有$f$个面，每一面是一个凸多边形，且任意两面不共面。将这枚骰子抛向空中，骰子落地的时候不会发生二次弹跳（这是一种非常理想的情况）。你希望知道最终每一面着地的概率。每一面着地的概率可以用如下的方法计算：我们假设$O$为骰子的重心，并以$O$为球心，做半径为$1$的单位球面（记为$S$）。我们知道$S$的表面积即单位球的表面积，为$4\pi$，这里$\pi$为圆周率。对于骰子的某一面$C$来说，球面$S$上存在一块区域$T$满足：当下落时若骰子所受重力方向与$S$的交点落在$T$中，则$C$就是最终着地的一面。那么$C$着地的概率为区域$T$的面积除以$4\pi$。为了能更好地辅助计算球面上一块区域的面积，我们给出单位球面$S$上三角形的面积计算公式。考虑单位球面$S$上的三个两两相交的大圆，交点依次为$A$，$B$和$C$。则曲面三角形$ABC$的面积为$\text{Area}(ABC)=\alpha+\beta+\gamma-\pi$，其中$\alpha$，$\beta$和$\gamma$分别对应了三个二面角的大小。如下图所示。

    我们保证：每一面着地的时候，重心的垂心都恰好在这一面内。也就是说不会出现摆不稳的情况。

---

题目分析

这是一道结论题。
需要知道以下几个东西才能做。

1. 如何求凸多面体重心？
   答：将凸多面体剖分成若干个四面体，每个四面体计算体积，所有剖分出的向量带权（体积）的复合平均即为重心。
2. 如何四面体剖分？
   答：对于每个面进行三角形剖分，与剖分点连线即形成四面体。
3. 如何求剖分出的四面体体积？
   答：剖分出的三个向量作混合积除以$6$既是体积。
   混合积指任意两个做叉积后于另外一个点积。
   叉积得到底面面积，点积得到与高的乘积，除以$6$即为体积。
4. 如何求二面角？
   答：通过法向量的夹角计算，法向量可以通过叉积得到。
5. 如何求曲面三角形面积？
   答：$\text{Area}(ABC)=\alpha+\beta+\gamma-\pi$，题目已描述。
6. 如何求曲面多边形面积？
   答：可以剖分成曲面三角形，也可以对曲面多边形的内角求和减去$(\text{点数}-2)\pi$（根据高斯波捏公式可推导）。

然后就可以做了。
然而这之前以上六个我都不会。

---

代码

#include<algorithm>
#include<iostream>
#include<iomanip>
#include<cstring>
#include<cstdlib>
#include<climits>
#include<vector>
#include<cstdio>
#include<cmath>
#include<queue>

using namespace std;

inline const int Get_Int() {
	int num=0,bj=1;
	char x=getchar();
	while(!isdigit(x)) {
		if(x=='-')bj=-1;
		x=getchar();
	}
	while(isdigit(x)) {
		num=num*10+x-'0';
		x=getchar();
	}
	return num*bj;
}

const int maxn=105;
const double pi=acos(-1);

struct Point {
	double x,y,z;
	Point(double x=0,double y=0,double z=0):x(x),y(y),z(z) {}
	Point operator + (const Point& b) {
		return Point(x+b.x,y+b.y,z+b.z);
	}
	void operator += (const Point& b) {
		*this=*this+b;
	}
	Point operator - (const Point& b) {
		return Point(x-b.x,y-b.y,z-b.z);
	}
	void operator -= (const Point& b) {
		*this=*this-b;
	}
	Point operator * (double c) {
		return Point(x*c,y*c,z*c);
	}
	void operator *= (double c) {
		*this=*this*c;
	}
	Point operator / (double c) {
		return Point(x/c,y/c,z/c);
	}
	void operator /= (double c) {
		*this=*this/c;
	}
	double length() {
		return sqrt(x*x+y*y+z*z);
	}
} p[maxn];

typedef Point Vector;

double Dot(Vector a,Vector b) {
	return a.x*b.x+a.y*b.y+a.z*b.z;
}

Vector Cross(Vector a,Vector b) {
	return Vector(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x);
}

double Triple(Vector a,Vector b,Vector c) {
	return Dot(a,Cross(b,c));
}

double Cal(Vector a,Vector b,Vector c) {
	Vector u=Cross(a,b),v=Cross(a,c);
	return acos(Dot(u,v)/u.length()/v.length());
}

int n,f,F[maxn][maxn];

int main() {
	n=Get_Int();
	f=Get_Int();
	for(int i=1; i<=n; i++) {
		double x,y,z;
		scanf("%lf%lf%lf",&x,&y,&z);
		p[i]=Point(x,y,z);
	}
	for(int i=1; i<=f; i++) {
		F[i][0]=Get_Int();
		for(int j=1; j<=F[i][0]; j++)F[i][j]=Get_Int();
	}
	double vol=0;
	Vector ctr;
	for(int i=1; i<=f; i++)
		for(int j=2; j<F[i][0]; j++) {
			Vector ci=p[F[i][1]]+p[F[i][j]]+p[F[i][j+1]];
			double vi=Triple(p[F[i][1]],p[F[i][j]],p[F[i][j+1]]);
			vol+=vi;
			ctr+=ci*vi;
		}
	ctr/=vol*4;
	for(int i=1; i<=n; i++)p[i]-=ctr;
	for(int i=1; i<=f; i++) {
		double ans=0;
		for(int j=2; j<F[i][0]; j++)ans+=Cal(p[F[i][j]],p[F[i][j-1]],p[F[i][j+1]]);
		ans+=Cal(p[F[i][F[i][0]]],p[F[i][F[i][0]-1]],p[F[i][1]])+Cal(p[F[i][1]],p[F[i][F[i][0]]],p[F[i][2]]);
		printf("%0.7lf\n",(ans-(F[i][0]-2)*pi)/(4*pi));
	}
	return 0;
}