三维计算几何

三维凸包

算法简述

三维凸包一个较为容易实现的算法是增量法。先将所有的点打乱顺序，然后选择四个不共面的点组成一个四面体，如果找不到说明凸包不存在。然后遍历剩余的点，不断更新凸包。对遍历到的点做如下处理。

1. 如果点在凸包内，则不更新。
2. 如果点在凸包外，那么找到所有原凸包上所有分隔了对于这个点可见面和不可见面的边，以这样的边的两个点和新的点创建新的面加入凸包中。

模板代码

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

/* Macros */
/******************************************************************************/
#define nxt(i) ((i+1)%n)
#define nxt2(i, x) ((i+1)%((x).size()))
#define prv(i) ((i+(x).size()-1)%n)
#define prv2(i, x) ((i+(x).size()-1)%((x).size()))
#define sz(x) (int((x).size()))
#define setpre(x) do{cout<<setprecision(x)<<setiosflags(ios::fixed);}while(0)

/* Real number tools */
/******************************************************************************/
const double PI = acos(-1.0);
const double eps = 1e-8;
double mysqrt(double x) {
    return x <= 0.0 ? 0.0 : sqrt(x);
}
double sq(double x) {
    return x*x;
}
int sgn(double x) {
    return x < -eps ? -1 : x > eps ? 1 : 0;
}

/* 3d Point */
/******************************************************************************/
struct Pt3 {
    double x, y, z;
    Pt3() { }
    Pt3(double x, double y, double z) : x(x), y(y), z(z) { }
};
typedef const Pt3 cPt3;
typedef cPt3 & cPt3r;

Pt3 operator + (cPt3r a, cPt3r b) { return Pt3(a.x+b.x, a.y+b.y, a.z+b.z); }
Pt3 operator - (cPt3r a, cPt3r b) { return Pt3(a.x-b.x, a.y-b.y, a.z-b.z); }
Pt3 operator * (cPt3r a, double A) { return Pt3(a.x*A, a.y*A, a.z*A); }
Pt3 operator * (double A, cPt3r a) { return Pt3(a.x*A, a.y*A, a.z*A); }
Pt3 operator / (cPt3r a, double A) { return Pt3(a.x/A, a.y/A, a.z/A); }
bool operator == (cPt3r a, cPt3r b) {
    return !sgn(a.x-b.x) && !sgn(a.y-b.y) && !sgn(a.z-b.z);
}
istream& operator >> (istream& sm, Pt3 &pt) {
    sm >> pt.x >> pt.y >> pt.z; return sm;
}
ostream & operator << (ostream& sm, cPt3r pt) {
    sm << "(" << pt.x << ", " << pt.y << ", " << pt.z << ")"; return sm;
}
double len(cPt3r p) { return mysqrt(sq(p.x) + sq(p.y) + sq(p.z)); }
double dist(cPt3r a, cPt3r b) { return len(a-b); }
Pt3 unit(cPt3r p) { return p / len(p); }
Pt3 det(cPt3r a, cPt3r b) {
    return Pt3(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 dot(cPt3r a, cPt3r b) {
    return a.x*b.x + a.y*b.y + a.z*b.z;
}
double mix(cPt3r a, cPt3r b, cPt3r c) {
    return dot(a, det(b, c));
}
/* 3d Line & Segment */
/******************************************************************************/
struct Ln3 {
    Pt3 a, b;
    Ln3() { }
    Ln3(cPt3r a, cPt3r b) : a(a), b(b) { }
};
typedef const Ln3 cLn3;
typedef cLn3 & cLn3r;

bool ptonln(cPt3r a, cPt3r b, cPt3r c) {
    return sgn(len(det(a-b, b-c))) <= 0;
}

/* 3d Plane */
/******************************************************************************/
struct Pl {
    Pt3 a, b, c;
    Pl() { }
    Pl(cPt3r a, cPt3r b, cPt3r c) : a(a), b(b), c(c) { }
};
typedef const Pl cPl;
typedef cPl & cPlr;

Pt3 nvec(cPlr pl) {
    return det(pl.a-pl.b, pl.b-pl.c);
}

/* Solution  */
/******************************************************************************/
bool cmp(cPt3r a, cPt3r b) {
    if (sgn(a.x-b.x)) return sgn(a.x-b.x) < 0;
    if (sgn(a.y-b.y)) return sgn(a.y-b.y) < 0;
    if (sgn(a.z-b.z)) return sgn(a.z-b.z) < 0;
    return false;
}

struct Face {
    int a, b, c;
    Face() { }
    Face(int a, int b, int c) : a(a), b(b), c(c) { }
};

void convex3d(vector<Pt3> &p, vector<Pl> &out) {
    sort(p.begin(), p.end(), cmp);
    p.erase(unique(p.begin(), p.end()), p.end());
    random_shuffle(p.begin(), p.end());
    vector<Face> face;
    for (int i = 2; i < sz(p); ++i) {
        if (ptonln(p[0], p[1], p[i])) continue;
        swap(p[i], p[2]);
        for (int j = i + 1; j < sz(p); ++j)
            if (sgn(mix(p[1]-p[0], p[2]-p[1], p[j]-p[0])) != 0) {
                swap(p[j], p[3]);
                face.push_back(Face(0, 1, 2));
                face.push_back(Face(0, 2, 1));
                goto found;
            }
    }
found:
    vector<vector<int> > mark(sz(p), vector<int>(sz(p), 0));
    for (int v = 3; v < sz(p); ++v) {
        vector<Face> tmp;
        for (int i = 0; i < sz(face); ++i) {
            int a = face[i].a, b = face[i].b, c = face[i].c;
            if (sgn(mix(p[a]-p[v], p[b]-p[v], p[c]-p[v])) < 0) {
                mark[a][b] = mark[b][a] = v;
                mark[b][c] = mark[c][b] = v;
                mark[c][a] = mark[a][c] = v;
            }else{
                tmp.push_back(face[i]);
            }
        }
        face = tmp;
        for (int i = 0; i < sz(tmp); ++i) {
            int a = face[i].a, b = face[i].b, c = face[i].c;
            if (mark[a][b] == v) face.push_back(Face(b, a, v));
            if (mark[b][c] == v) face.push_back(Face(c, b, v));
            if (mark[c][a] == v) face.push_back(Face(a, c, v));
        }
    }
    out.clear();
    for (int i = 0; i < sz(face); ++i)
        out.push_back(Pl(p[face[i].a], p[face[i].b], p[face[i].c]));
}

vector<Pt3> p;
vector<Pl> out;

int main() {
    int n;
    cin >> n;
    for (int i = 0; i < n; ++i) {
        Pt3 pt;
        cin >> pt;
        p.push_back(pt);
    }
    convex3d(p, out);
    double area = 0.0;
    for (int i = 0; i < sz(out); ++i)
        area += len(det(out[i].a-out[i].b, out[i].b-out[i].c));
    setpre(3);
    cout << area / 2.0 << "\n";
    return 0;
}

经典题目

• POJ 3528 Ultimate Weapon