P4724 【模板】三維凸包
阿新 • • 發佈:2019-01-01
\(\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; }