1. 程式人生 > >P4724 【模板】三維凸包

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