三维凸包的求解如果有四点共面，那么就比较繁琐了。
增量法，每次添加一个点，然后更新凸包。
关于增量法的图示看这里:
http://www.cse.unsw.edu.au/~lambert/java/3d/hull.html
现在的主要问题就是，一直所有面的法向量，如何更新凸包。
法向量点乘一个向量，可以得到这个向量在法向量上的偏移，如果这个偏移是正的，说明这个向量相关的点在面的上方，负的是下方，所以只要找到一条边关联的两个面的点积是一正一负即可。
注意这里的所有面都是由三个点组成，且是有序的，每个边只被记录两次，一次顺时针，一次逆时针。
这里不能使用int，数值太大会越界。
1
2 const int N = 1024;
3 struct F {
4 int v[4];
5 F(){}
6 F(int a,int b,int c){ v[0] = a, v[1] = b, v[2] = c, v[3] = v[0];}
7 };
8
9 struct point_t {
10 double x, y, z;
11 point_t (){}
12 point_t (double a, double b, double c){x = a, y = b, z = c;}
13 };
14 point_t operator +(point_t a, point_t b) { return point_t (a.x + b.x, a.y + b.y, a.z + b.z);}
15 point_t operator -(point_t a, point_t b) { return point_t (a.x - b.x, a.y - b.y, a.z - b.z);}
16
17 double dot_mul(point_t a, point_t b) { return a.x * b.x + a.y * b.y + a.z * b.z;}
18 point_t cross_mul(point_t a,point_t b)
19 {
20 return point_t (a.y*b.z - b.y*a.z ,
21 a.z*b.x - b.z*a.x ,
22 a.x*b.y - b.x*a.y);
23 }
24
25 double area(point_t a, point_t b)
26 {
27 point_t v = cross_mul(a, b);
28 return sqrt(0.0 + v.x * v.x + v.y * v.y + v.z * v.z) / 2;
29 }
30
31 point_t p[N];
32 int vis[N][N], n;
33
34 int main()
35 {
36 int i, j, k;
37 scanf("%d", &n);
38 for (i = 0;i < n;i++) {
39 scanf("%lf %lf %lf", &p[i].x, &p[i].y, &p[i].z);
40 }
41 vector<F> cur;
42 cur.pb(F(0, 1, 2)), cur.pb(F(2, 1, 0));
43 for (i = 3;i < n;i++) {
44 vector<F> next;
45 for (j = 0;j < cur.size();j++) {
46 F f = cur[j];
47 point_t v = cross_mul(p[f.v[1]] - p[f.v[0]],
48 p[f.v[2]] - p[f.v[1]]);
49 int val = dot_mul(v, p[i] - p[f.v[1]]) < 0 ? -1 : 1;
50 if (val < 0) {
51 next.pb(f);
52 }
53 for (k = 0;k < 3;k++) {
54 if (vis[f.v[k+1]][f.v[k]] == 0) {
55 vis[f.v[k]][f.v[k+1]] = val;
56 }else { //http://www.cppblog.com/schindlerlee
57 if (vis[f.v[k+1]][f.v[k]] != val) {
58 if (val > 0) {
59 next.pb(F(f.v[k], f.v[k+1], i));
60 }else {
61 next.pb(F(f.v[k+1], f.v[k], i));
62 }
63 }
64 vis[f.v[k+1]][f.v[k]] = 0;
65 }
66 }
67 }
68 cur = next;
69 }
70 double ans = 0;
71 for (i = 0;i < cur.size();i++) {
72 F f = cur[i];
73 ans += area(p[f.v[0]] - p[f.v[1]], p[f.v[2]] - p[f.v[1]]);
74 }
75 printf("%.3f\n", ans);
76 return 0;
77 }