Blog of RuSun

\begin {array}{c} \mathfrak {One Problem Is Difficult} \\\\ \mathfrak {Because You Don't Know} \\\\ \mathfrak {Why It Is Diffucult} \end {array}

三维凸包模板

用于求三维凸包的面积。

查看代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <bitset>
using namespace std;
const int N = 2e3 + 10;
const double eps = 1e-8, pi = acos(-1);
double rand_eps()
{
return ((double)rand() / RAND_MAX - 0.5) * eps;
}
struct Point
{
double x, y, z;
void input()
{
scanf("%lf%lf%lf", &x, &y, &z);
}
void rd()
{
x += rand_eps(), y += rand_eps(), z += rand_eps();
}
Point operator-(Point t) const
{
return {x - t.x, y - t.y, z - t.z};
}
Point operator+(Point t) const
{
return {x + t.x, y + t.y, z + t.z};
}
double operator^(Point t) const
{
return x * t.x + y * t.y + z * t.z;
}
Point operator&(Point t) const
{
return {y * t.z - z * t.y, z * t.x - x * t.z, x * t.y - y * t.x};
}
double len()
{
return sqrt(x * x + y * y + z * z);
}
} p[N];
struct Plain
{
int v[3];
Point normal()
{
return (p[v[1]] - p[v[0]]) & (p[v[2]] - p[v[0]]);
}
double area()
{
return normal().len() / 2;
}
bool above(Point a)
{
return ((a - p[v[0]]) ^ normal()) >= 0;
}
} pl[N], np[N];
int ConvexHull(int n)
{
static bitset<N> g[N];
int cnt = 0, ncnt = 0;
pl[++cnt] = {1, 2, 3};
pl[++cnt] = {3, 2, 1};
for (int i = 4; i <= n; ++i)
{
ncnt = 0;
for (int j = 1; j <= cnt; ++j)
{
bool t = pl[j].above(p[i]);
if (!t)
np[++ncnt] = pl[j];
for (int k = 0; k < 3; ++k)
g[pl[j].v[k]][pl[j].v[(k + 1) % 3]] = t;
}
for (int j = 1; j <= cnt; ++j)
for (int k = 0; k < 3; ++k)
{
int a = pl[j].v[k], b = pl[j].v[(k + 1) % 3];
if (g[a][b] && !g[b][a])
np[++ncnt] = {a, b, i};
}
cnt = ncnt;
for (int j = 1; j <= cnt; ++j)
pl[j] = np[j];
}
return cnt;
}
int n;
int main()
{
scanf("%d", &n);
for (int i = 1; i <= n; ++i)
p[i].input();
for (int i = 1; i <= n; ++i)
p[i].rd();
int m = ConvexHull(n);
double res = 0;
for (int i = 1; i <= m; ++i)
res += pl[i].area();
printf("%.3lf\n", res);
return 0;
}