P4724 【模板】三维凸包

(color{#0066ff}{题目描述})

给出空间中n个点，求凸包表面积。

(color{#0066ff}{输入格式})

第一行一个整数n，表示点数。

接下来n行，每行三个实数x，y，z描述坐标。

(color{#0066ff}{输出格式})

输出凸包表面积，保留3位小数。

(color{#0066ff}{输入样例})

4 
0 0 0
1 0 0
0 1 0
0 0 1

(color{#0066ff}{输出样例})

2.366

(color{#0066ff}{数据范围与提示})

n≤2000

(color{#0066ff}{题解})

增量法

把每个面，分成正面，反面

先选出3个点（构成一个面）

每次加点

先把当前的点能看见的面全部删除（最后的凸包一定不存在能被某个点看见的面）

然后枚举前面的面中的某两个点，与当前点构成新面，成立则加入

最后的就是凸包的面‘

为了防止共面共线问题，可以在精度允许范围内微调一下坐标

#include<bits/stdc++.h>
using namespace std;
#define LL long long
LL in() {
	char ch; int x = 0, f = 1;
	while(!isdigit(ch = getchar()))(ch == '-') && (f = -f);
	for(x = ch ^ 48; isdigit(ch = getchar()); x = (x << 1) + (x << 3) + (ch ^ 48));
	return x * f;
}
const double eps = 1e-9;
const int maxn = 2050;
struct node {
	double x, y, z;
	node(double x = 0, double y = 0, double z = 0): x(x), y(y), z(z) {}
	node operator - (const node &b) const {
		return node(x - b.x, y - b.y, z - b.z);
	}
	node operator ^ (const node &b) const {
		return node(y * b.z - z * b.y, z * b.x - x * b.z, x * b.y - y * b.x);
	}
	double operator * (const node &b) const {
		return x * b.x + y * b.y + z * b.z;
	}
	void init() {
		x = x + ((double)rand() / (double)RAND_MAX - 0.5) * eps * 10;
		y = y + ((double)rand() / (double)RAND_MAX - 0.5) * eps * 10;
		z = z + ((double)rand() / (double)RAND_MAX - 0.5) * eps * 10;
	}
	double mo() {
		return sqrt(*this * *this);
	}
	bool jud() {
		return fabs(x) <= eps && fabs(y) <= eps && fabs(z) <= eps;
	}
}e[maxn];

struct plane {
	int v[3];
	plane(int a = 0, int b = 0, int c = 0) { v[0] = a, v[1] = b, v[2] = c; }
	int &operator [] (const int &b) {
		return v[b];
	}
	node F() const {
		return ((e[v[1]] - e[v[0]]) ^ (e[v[2]] - e[v[0]]));
	}
	bool cansee(node x) const {
		return (x - e[v[0]]) * F() > 0;
	}
};
int n, cnt;
bool vis[maxn][maxn];
void init() {
	n = in();
	for(int i = 1; i <= n; i++) {
		node o;
		scanf("%lf%lf%lf", &o.x, &o.y, &o.z);
		for(int j = 1; j <= cnt; j++) if((e[j] - o).jud()) goto cant;
		e[++cnt] = o;
		cant:;
	}
	n = cnt;
	for(int i = 1; i <= n; i++) e[i].init();
}
double D() {
	double ans = 0;
	using std::vector;
	vector<plane> c;
	c.push_back(plane(1, 2, 3));
	c.push_back(plane(3, 2, 1));
	for(int i = 4; i <= n; i++) {
		vector<plane> q;
		for(int j = 0; j < (int)c.size(); j++) {
			plane t = c[j];
			bool flag = t.cansee(e[i]);
			if(!flag) q.push_back(c[j]);
			for(int k = 0; k < 3; k++)
				vis[t[k]][t[(k + 1) % 3]] = flag;
		}
		for(int j = 0; j < (int)c.size(); j++) 
			for(int k = 0; k < 3; k++) {
				int a = c[j][k], b = c[j][(k + 1) % 3];
				if(vis[a][b] != vis[b][a] && vis[a][b]) 
					q.push_back(plane(a, b, i));
			}
		c = q;
	}
	for(int i = 0; i < (int)c.size(); i++) ans += c[i].F().mo();
	return ans;
}
	
			
int main() {
	init();
	printf("%.3f", D() / 2.0);
	return 0;
}