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