【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;
}