poj3528--Ultimate Weapon（三维凸包）
题目意思：求三维凸包的表面积，模板题
View Code
1 /*给出三维空间中的n个顶点,求解由这n个顶点构成的凸包表面的多边形个数. 2 增量法求解:首先任选4个点形成的一个四面体,然后每次新加一个点,分两种情况: 3 1> 在凸包内,则可以跳过 4 2> 在凸包外,找到从这个点可以"看见"的面,删除这些面, 5 然后对于一边没有面的线段,和新加的这个点新建一个面,至于这个点可以看见的面, 6 就是求出这个面的方程(可以直接求法向量). 7 下面是三维凸包的模板。。有了这个模板应该对付三维凸包的题就没问题了吧。。*/ 8 //Accepted 656K 0MS C++ 7720B 9 #include<iostream> 10 #include<cmath> 11 #include<cstring> 12 #include<cstdlib> 13 #include<algorithm> 14 using namespace std; 15 const int MAXN=505; 16 const double EPS=1e-8; 17 struct Point 18 { 19 double x,y,z; 20 Point(){} 21 Point(double xx,double yy,double zz):x(xx),y(yy),z(zz){} 22 23 Point operator -(const Point p1) //两向量之差 24 { 25 return Point(x-p1.x,y-p1.y,z-p1.z); 26 } 27 28 Point operator *(Point p) //叉乘 29 { 30 return Point(y*p.z-z*p.y,z*p.x-x*p.z,x*p.y-y*p.x); 31 } 32 33 double operator ^(Point p) //点乘 34 { 35 return (x*p.x+y*p.y+z*p.z); 36 } 37 }; 38 struct CH3D 39 { 40 struct face 41 { 42 int a,b,c; //表示凸包一个面上三个点的编号 43 bool ok; //表示该面是否属于最终凸包中的面 44 }; 45 46 int n; //初始顶点数 47 Point P[MAXN]; //初始顶点 48 49 int num; //凸包表面的三角形数 50 face F[8*MAXN]; 51 52 int g[MAXN][MAXN]; //凸包表面的三角形 53 54 double vlen(Point a) //向量长度 55 { 56 return sqrt(a.x*a.x+a.y*a.y+a.z*a.z); 57 } 58 59 Point cross(const Point &a, const Point &b, const Point &c) //叉乘 60 { 61 return Point((b.y-a.y)*(c.z-a.z)-(b.z-a.z)*(c.y-a.y),-((b.x-a.x)*(c.z-a.z) 62 -(b.z-a.z)*(c.x-a.x)),(b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x)); 63 } 64 double area(Point a,Point b,Point c) //三角形面积*2 65 { 66 return vlen((b-a)*(c-a)); 67 } 68 69 double volume(Point a,Point b,Point c,Point d) //四面体有向体积*6 70 { 71 return (b-a)*(c-a)^(d-a); 72 } 73 74 double dblcmp(Point &p,face &f) //正:点在面同向 75 { 76 Point m=P[f.b]-P[f.a]; 77 Point n=P[f.c]-P[f.a]; 78 Point t=p-P[f.a]; 79 return (m*n)^t; 80 } 81 82 void deal(int p,int a,int b) 83 { 84 int f=g[a][b]; 85 face add; 86 if(F[f].ok) 87 { 88 if(dblcmp(P[p],F[f])>EPS) 89 dfs(p,f); 90 else 91 { 92 add.a=b; 93 add.b=a; 94 add.c=p; 95 add.ok=1; 96 g[p][b]=g[a][p]=g[b][a]=num; 97 F[num++]=add; 98 } 99 } 100 } 101 102 void dfs(int p,int now) 103 { 104 F[now].ok=0; 105 deal(p,F[now].b,F[now].a); 106 deal(p,F[now].c,F[now].b); 107 deal(p,F[now].a,F[now].c); 108 } 109 110 bool same(int s,int t) 111 { 112 Point &a=P[F[s].a]; 113 Point &b=P[F[s].b]; 114 Point &c=P[F[s].c]; 115 return fabs(volume(a,b,c,P[F[t].a]))<EPS && fabs(volume(a,b,c,P[F[t].b]))<EPS 116 && fabs(volume(a,b,c,P[F[t].c]))<EPS; 117 } 118 119 void solve() //构建三维凸包 120 { 121 int i,j,tmp; 122 face add; 123 bool flag=true; 124 num=0; 125 if(n<4) 126 return; 127 for(i=1;i<n;i++) //此段是为了保证前四个点不共面,若以保证,则可去掉 128 { 129 if(vlen(P[0]-P[i])>EPS) 130 { 131 swap(P[1],P[i]); 132 flag=false; 133 break; 134 } 135 } 136 if(flag) 137 return; 138 flag=true; 139 for(i=2;i<n;i++) //使前三点不共线 140 { 141 if(vlen((P[0]-P[1])*(P[1]-P[i]))>EPS) 142 { 143 swap(P[2],P[i]); 144 flag=false; 145 break; 146 } 147 } 148 if(flag) 149 return; 150 flag=true; 151 for(i=3;i<n;i++) //使前四点不共面 152 { 153 if(fabs((P[0]-P[1])*(P[1]-P[2])^(P[0]-P[i]))>EPS) 154 { 155 swap(P[3],P[i]); 156 flag=false; 157 break; 158 } 159 } 160 if(flag) 161 return; 162 for(i=0;i<4;i++) 163 { 164 add.a=(i+1)%4; 165 add.b=(i+2)%4; 166 add.c=(i+3)%4; 167 add.ok=true; 168 if(dblcmp(P[i],add)>0) 169 swap(add.b,add.c); 170 g[add.a][add.b]=g[add.b][add.c]=g[add.c][add.a]=num; 171 F[num++]=add; 172 } 173 for(i=4;i<n;i++) 174 { 175 for(j=0;j<num;j++) 176 { 177 if(F[j].ok && dblcmp(P[i],F[j])>EPS) 178 { 179 dfs(i,j); 180 break; 181 } 182 } 183 } 184 tmp=num; 185 for(i=num=0;i<tmp;i++) 186 if(F[i].ok) 187 { 188 F[num++]=F[i]; 189 } 190 } 191 192 double area() //表面积 193 { 194 double res=0.0; 195 if(n==3) 196 { 197 Point p=cross(P[0],P[1],P[2]); 198 res=vlen(p)/2.0; 199 return res; 200 } 201 for(int i=0;i<num;i++) 202 res+=area(P[F[i].a],P[F[i].b],P[F[i].c]); 203 return res/2.0; 204 } 205 206 double volume() //体积 207 { 208 double res=0.0; 209 Point tmp(0,0,0); 210 for(int i=0;i<num;i++) 211 res+=volume(tmp,P[F[i].a],P[F[i].b],P[F[i].c]); 212 return fabs(res/6.0); 213 } 214 215 int triangle() //表面三角形个数 216 { 217 return num; 218 } 219 220 int polygon() //表面多边形个数 221 { 222 int i,j,res,flag; 223 for(i=res=0;i<num;i++) 224 { 225 flag=1; 226 for(j=0;j<i;j++) 227 if(same(i,j)) 228 { 229 flag=0; 230 break; 231 } 232 res+=flag; 233 } 234 return res; 235 } 236 }; 237 CH3D hull; 238 int main() 239 { 240 int i; 241 double res; 242 while(scanf("%d",&hull.n)!=EOF) 243 { 244 for(i=0;i<hull.n;i++) 245 scanf("%lf%lf%lf",&hull.P[i].x,&hull.P[i].y,&hull.P[i].z); 246 hull.solve(); 247 res=hull.area(); 248 printf("%.3lf\n",res); 249 } 250 return 0; 251 }
