直接从别人那复制的代码，惭愧啊，看懂了为自己储存模板。
View Code
1 #include<iostream> 2 #include<cmath> 3 #include<cstring> 4 #include<cstdio> 5 #define N 505 6 #define eps 0.000001 7 using namespace std; 8 struct Point 9 { 10 double x,y,z; 11 Point(){} 12 Point(double _x,double _y,double _z) 13 { 14 x=_x; 15 y=_y; 16 z=_z; 17 } 18 Point operator-(Point t1)//向量减法 19 { 20 return Point(x-t1.x,y-t1.y,z-t1.z); 21 } 22 Point operator*(Point t2)//叉积 23 { 24 return Point(y*t2.z-t2.y*z,z*t2.x-x*t2.z,x*t2.y-y*t2.x); 25 } 26 double operator^(Point t3)//点积 27 { 28 return x*t3.x+y*t3.y+z*t3.z; 29 } 30 }; 31 struct Plane 32 { 33 int a,b,c;//a,b,c为三个点的编号，a,b,c要满足从凸包外面看成右手系 34 bool in;//表示该平面是否在凸包内 35 }; 36 void Swap(Point &a,Point &b) 37 { 38 Point c; 39 c=a; 40 a=b; 41 b=c; 42 } 43 Point point[N]; 44 Plane plane[N*10]; 45 int edge[N][N]; 46 int plen;//计算过的面的个数 47 void dfs(int p,int t); 48 double vol(Point p,Plane f)//p与平面abc组成的四面体的有向体积的倍 49 { 50 Point a=point[f.a]-p,b=point[f.b]-p,c=point[f.c]-p; 51 return (a*b)^c; 52 } 53 double vlen(Point a)//求向量a的模 54 { 55 return sqrt(a.x*a.x+a.y*a.y+a.z*a.z); 56 } 57 void deal(int p,int t1,int t2) 58 { 59 int t=edge[t1][t2];//搜索与该边相邻的另外一个平面 60 if(plane[t].in) 61 { 62 if(vol(point[p],plane[t])>eps) 63 dfs(p,t); 64 else 65 { 66 Plane add; 67 add.a=t2,add.b=t1,add.c=p,add.in=true;//这里注意顺序，就可以不用Swap了，add.a,add.b,add.c要成右手系 68 edge[add.a][add.b]=edge[add.b][add.c]=edge[add.c][add.a]=plen; 69 plane[plen++]=add; 70 } 71 } 72 } 73 void dfs(int p,int t)//递归搜索所有应该从凸包内删除的面 74 { 75 plane[t].in=false; 76 deal(p,plane[t].b,plane[t].a);//注意，a和b的顺序刚好跟下面的相反，为的就是搜索与边(point[plane[t].a],point[plane[t].b])相邻的另外一个平面 77 deal(p,plane[t].c,plane[t].b); 78 deal(p,plane[t].a,plane[t].c); 79 } 80 int del(int n)//增量法，有n个点，返回计算过的平面个数，若无法构成凸包，则返回-1 81 { 82 if(n<4)//如果点数小于，则无法构成凸包，若已保证点数大于或等于，可略去 83 return -1; 84 /******************这一段用来保证前四点不共面，若已保证，可略去 85 bool allTheSamePoint=true; 86 for(int i=1;i<n;i++)//保证前两点不共点 87 { 88 if(vlen(point[i]-point[0])>eps) 89 { 90 Swap(point[1],point[i]); 91 allTheSamePoint=false; 92 break; 93 } 94 } 95 if(allTheSamePoint) 96 return -1; 97 bool allTheSameLine=true; 98 for(int i=2;i<n;i++)//保证前三点不共线 99 { 100 if(vlen((point[1]-point[0])*(point[i]-point[1]))>eps) 101 { 102 Swap(point[2],point[i]); 103 allTheSameLine=false; 104 break; 105 } 106 } 107 if(allTheSameLine) 108 return -1; 109 bool allTheSamePlane=true; 110 for(int i=3;i<n;i++)//保证前四点不共面 111 { 112 if(fabs((point[1]-point[0])*(point[2]-point[0])^(point[i]-point[0]))>eps) 113 { 114 Swap(point[3],point[i]); 115 allTheSamePlane=false; 116 break; 117 } 118 } 119 if(allTheSamePlane) 120 return -1; 121 这一段用来保证前四点不共面，若已保证，可略去************/ 122 plen=0;//计算过的面的个数 123 Plane add; 124 for(int i=0;i<4;i++) 125 { 126 add.a=(i+1)%4,add.b=(i+2)%4,add.c=(i+3)%4,add.in=true; 127 if(vol(point[i],add)>0) 128 swap(add.a,add.b); 129 edge[add.a][add.b]=edge[add.b][add.c]=edge[add.c][add.a]=plen;//记录与该边相邻的其中一个面，并且该顺序在该面内（从凸包外看）成右手系，因此，该面是唯一的 130 plane[plen++]=add; 131 } 132 for(int i=4;i<n;i++) 133 { 134 for(int j=0;j<plen;j++) 135 { 136 if(plane[j].in && vol(point[i],plane[j])>eps) 137 { 138 dfs(i,j); 139 break; 140 } 141 } 142 } 143 return plen; 144 } 145 double area(Plane a) 146 { 147 return vlen((point[a.b]-point[a.a])*(point[a.c]-point[a.a]))/2.0; 148 } 149 150 int main() 151 { 152 int n; 153 cin>>n; 154 for(int i=0;i<n;i++) 155 cin>>point[i].x>>point[i].y>>point[i].z; 156 int len=del(n); 157 if(len==-1) 158 printf("0.000\n"); 159 else 160 { 161 double ans=0.0; 162 for(int i=0;i<len;i++) 163 { 164 if(plane[i].in) 165 { 166 ans+=area(plane[i]); 167 } 168 } 169 printf("%.3lf\n",ans); 170 } 171 return 0; 172 }