【HDU 3662】3D Convex Hull

http://acm.hdu.edu.cn/showproblem.php?pid=3662
求给定空间中的点的三维凸包上有多少个面。
用增量法，不断加入点，把新加的点能看到的面都删掉，不能看到的面与能看到的面的棱与新点相连构成一个新的三角形面。
这样的面全都是三角形，注意最后统计答案时要把重合的面算成一个。
时间复杂度\(O(n^2)\)。

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;

const int N = 303;
const double eps = 1e-6;

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 &A) const {
		return Point(x + A.x, y + A.y, z + A.z);
	}
	Point operator - (const Point &A) const {
		return Point(x - A.x, y - A.y, z - A.z);
	}
	double operator * (const Point &A) const {
		return x * A.x + y * A.y + z * A.z;
	}
	Point operator ^ (const Point &A) const {
		return Point(y * A.z - z * A.y, z * A.x - x * A.z, x * A.y - y * A.x);
	}
	double sqrlen() {
		return x * x + y * y + z * z;
	}
} P[N];

struct Face {
	int a, b, c; bool ex;
	Face(int _a = 0, int _b = 0, int _c = 0, bool _ex = false) : a(_a), b(_b), c(_c), ex(_ex) {}
} F[N * N];

int n, ftot, LeftFace[N][N];

void insFace(int a, int b, int c, int n1, int n2, int n3) {
	F[++ftot] = (Face) {a, b, c, true};
	LeftFace[a][b] = LeftFace[b][c] = LeftFace[c][a] = ftot;
	LeftFace[b][a] = n1;
	LeftFace[c][b] = n2;
	LeftFace[a][c] = n3;
}

bool visible(int f, int p) {
	Point a = P[F[f].b] - P[F[f].a], b = P[F[f].c] - P[F[f].a];
	return (P[p] - P[F[f].a]) * (a ^ b) > eps;
}

int st, to[N], ps[N], pt[N], ptot = 0, pf[N];

void dfs(int x, int s, int t, int p) {
	if (F[x].ex == false) return;
	
	if (visible(x, p))
		F[x].ex = false;
	else {
		to[st = s] = t;
		return;
	}
	
	dfs(LeftFace[F[x].b][F[x].a], F[x].a, F[x].b, p);
	dfs(LeftFace[F[x].c][F[x].b], F[x].b, F[x].c, p);
	dfs(LeftFace[F[x].a][F[x].c], F[x].c, F[x].a, p);
}

Point ff;
void dfs2(int x) {
	if (F[x].ex == false) return;
	
	Point now = (P[F[x].b] - P[F[x].a]) ^ (P[F[x].c] - P[F[x].a]);
	if (fabs(now * ff - sqrt(now.sqrlen() * ff.sqrlen())) < 1e-6)
		F[x].ex = false;
	else
		return;
	
	dfs2(LeftFace[F[x].b][F[x].a]);
	dfs2(LeftFace[F[x].c][F[x].b]);
	dfs2(LeftFace[F[x].a][F[x].c]);
}

int main() {
	while (~scanf("%d", &n)) {
		for (int i = 1; i <= n; ++i) scanf("%lf%lf%lf", &P[i].x, &P[i].y, &P[i].z);
		ftot = 0;
		Point a, b, c, d, e;
		a = P[2] - P[1];
		int tmp, id1, id2;
		for (tmp = 3; tmp <= n; ++tmp) {
			b = P[tmp] - P[1];
			d = a ^ b;
			if (d.sqrlen() < eps) continue;
			id1 = tmp; break;
		}
		for (++tmp; tmp <= n; ++tmp) {
			c = P[tmp] - P[1];
			if (fabs(d * c) < eps) continue;
			id2 = tmp; break;
		}
		
		if (d * c < 0) swap(id1, id2);
		insFace(1, 2, id2, 3, 4, 2);
		insFace(1, id2, id1, 1, 4, 3);
		insFace(1, id1, 2, 2, 4, 1);
		insFace(2, id1, id2, 3, 2, 1);
		
		for (int i = 3; i <= n; ++i) {
			if (i == id1 || i == id2) continue;
			for (int j = 1; j <= ftot; ++j)
				if (F[j].ex && visible(j, i)) {
					dfs(j, 0, 0, i);
					ptot = 0;
					int tmps = st, tmpt = to[st], ppff = ftot;
					do {
						++ptot;
						ps[ptot] = tmps; pt[ptot] = tmpt;
						pf[ptot] = ++ppff;
						tmps = tmpt; tmpt = to[tmpt];
					} while (tmps != st);
					
					for (int k = 1, pre, suc; k <= ptot; ++k) {
						pre = k - 1; suc = k + 1;
						if (pre == 0) pre = ptot;
						if (suc > ptot) suc = 1;
						pre = pf[pre]; suc = pf[suc];
						insFace(pt[k], i, ps[k], suc, pre, LeftFace[pt[k]][ps[k]]);
					}
					
					break;
				}
		}
		
		int ans = 0;
		for (int i = 1; i <= ftot; ++i)
			if (F[i].ex) {
				++ans;
				ff = (P[F[i].b] - P[F[i].a]) ^ (P[F[i].c] - P[F[i].a]);
				dfs2(i);
			}
		printf("%d\n", ans);
	}
	
	return 0;
}