新模板
1 /* 2 HDU 4273 Rescue 3 给一个三维凸包，求重心到表面的最短距离 4 模板题：三维凸包+多边形重心+点面距离 5 */ 6 7 #include<stdio.h> 8 #include<algorithm> 9 #include<string.h> 10 #include<math.h> 11 #include<stdlib.h> 12 using namespace std; 13 const int MAXN=550; 14 const double eps=1e-8; 15 16 struct Point 17 { 18 double x,y,z; 19 Point(){} 20 21 Point(double xx,double yy,double zz):x(xx),y(yy),z(zz){} 22 23 //两向量之差 24 Point operator -(const Point p1) 25 { 26 return Point(x-p1.x,y-p1.y,z-p1.z); 27 } 28 29 //两向量之和 30 Point operator +(const Point p1) 31 { 32 return Point(x+p1.x,y+p1.y,z+p1.z); 33 } 34 35 //叉乘 36 Point operator *(const Point p) 37 { 38 return Point(y*p.z-z*p.y,z*p.x-x*p.z,x*p.y-y*p.x); 39 } 40 41 Point operator *(double d) 42 { 43 return Point(x*d,y*d,z*d); 44 } 45 46 Point operator / (double d) 47 { 48 return Point(x/d,y/d,z/d); 49 } 50 51 //点乘 52 double operator ^(Point p) 53 { 54 return (x*p.x+y*p.y+z*p.z); 55 } 56 }; 57 58 struct CH3D 59 { 60 struct face 61 { 62 //表示凸包一个面上的三个点的编号 63 int a,b,c; 64 //表示该面是否属于最终凸包上的面 65 bool ok; 66 }; 67 //初始顶点数 68 int n; 69 //初始顶点 70 Point P[MAXN]; 71 //凸包表面的三角形数 72 int num; 73 //凸包表面的三角形 74 face F[8*MAXN]; 75 //凸包表面的三角形 76 int g[MAXN][MAXN]; 77 //向量长度 78 double vlen(Point a) 79 { 80 return sqrt(a.x*a.x+a.y*a.y+a.z*a.z); 81 } 82 //叉乘 83 Point cross(const Point &a,const Point &b,const Point &c) 84 { 85 return Point((b.y-a.y)*(c.z-a.z)-(b.z-a.z)*(c.y-a.y), 86 (b.z-a.z)*(c.x-a.x)-(b.x-a.x)*(c.z-a.z), 87 (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x) 88 ); 89 } 90 //三角形面积*2 91 double area(Point a,Point b,Point c) 92 { 93 return vlen((b-a)*(c-a)); 94 } 95 //四面体有向体积*6 96 double volume(Point a,Point b,Point c,Point d) 97 { 98 return (b-a)*(c-a)^(d-a); 99 } 100 //正：点在面同向 101 double dblcmp(Point &p,face &f) 102 { 103 Point m=P[f.b]-P[f.a]; 104 Point n=P[f.c]-P[f.a]; 105 Point t=p-P[f.a]; 106 return (m*n)^t; 107 } 108 void deal(int p,int a,int b) 109 { 110 int f=g[a][b];//搜索与该边相邻的另一个平面 111 face add; 112 if(F[f].ok) 113 { 114 if(dblcmp(P[p],F[f])>eps) 115 dfs(p,f); 116 else 117 { 118 add.a=b; 119 add.b=a; 120 add.c=p;//这里注意顺序，要成右手系 121 add.ok=true; 122 g[p][b]=g[a][p]=g[b][a]=num; 123 F[num++]=add; 124 } 125 } 126 } 127 void dfs(int p,int now)//递归搜索所有应该从凸包内删除的面 128 { 129 F[now].ok=0; 130 deal(p,F[now].b,F[now].a); 131 deal(p,F[now].c,F[now].b); 132 deal(p,F[now].a,F[now].c); 133 } 134 bool same(int s,int t) 135 { 136 Point &a=P[F[s].a]; 137 Point &b=P[F[s].b]; 138 Point &c=P[F[s].c]; 139 return fabs(volume(a,b,c,P[F[t].a]))<eps && 140 fabs(volume(a,b,c,P[F[t].b]))<eps && 141 fabs(volume(a,b,c,P[F[t].c]))<eps; 142 } 143 //构建三维凸包 144 void create() 145 { 146 int i,j,tmp; 147 face add; 148 149 num=0; 150 if(n<4)return; 151 //********************************************** 152 //此段是为了保证前四个点不共面 153 bool flag=true; 154 for(i=1;i<n;i++) 155 { 156 if(vlen(P[0]-P[i])>eps) 157 { 158 swap(P[1],P[i]); 159 flag=false; 160 break; 161 } 162 } 163 if(flag)return; 164 flag=true; 165 //使前三个点不共线 166 for(i=2;i<n;i++) 167 { 168 if(vlen((P[0]-P[1])*(P[1]-P[i]))>eps) 169 { 170 swap(P[2],P[i]); 171 flag=false; 172 break; 173 } 174 } 175 if(flag)return; 176 flag=true; 177 //使前四个点不共面 178 for(int i=3;i<n;i++) 179 { 180 if(fabs((P[0]-P[1])*(P[1]-P[2])^(P[0]-P[i]))>eps) 181 { 182 swap(P[3],P[i]); 183 flag=false; 184 break; 185 } 186 } 187 if(flag)return; 188 //***************************************** 189 for(i=0;i<4;i++) 190 { 191 add.a=(i+1)%4; 192 add.b=(i+2)%4; 193 add.c=(i+3)%4; 194 add.ok=true; 195 if(dblcmp(P[i],add)>0)swap(add.b,add.c); 196 g[add.a][add.b]=g[add.b][add.c]=g[add.c][add.a]=num; 197 F[num++]=add; 198 } 199 for(i=4;i<n;i++) 200 { 201 for(j=0;j<num;j++) 202 { 203 if(F[j].ok&&dblcmp(P[i],F[j])>eps) 204 { 205 dfs(i,j); 206 break; 207 } 208 } 209 } 210 tmp=num; 211 for(i=num=0;i<tmp;i++) 212 if(F[i].ok) 213 F[num++]=F[i]; 214 215 } 216 //表面积 217 double area() 218 { 219 double res=0; 220 if(n==3) 221 { 222 Point p=cross(P[0],P[1],P[2]); 223 res=vlen(p)/2.0; 224 return res; 225 } 226 for(int i=0;i<num;i++) 227 res+=area(P[F[i].a],P[F[i].b],P[F[i].c]); 228 return res/2.0; 229 } 230 double volume() 231 { 232 double res=0; 233 Point tmp(0,0,0); 234 for(int i=0;i<num;i++) 235 res+=volume(tmp,P[F[i].a],P[F[i].b],P[F[i].c]); 236 return fabs(res/6.0); 237 } 238 //表面三角形个数 239 int triangle() 240 { 241 return num; 242 } 243 //表面多边形个数 244 int polygon() 245 { 246 int i,j,res,flag; 247 for(i=res=0;i<num;i++) 248 { 249 flag=1; 250 for(j=0;j<i;j++) 251 if(same(i,j)) 252 { 253 flag=0; 254 break; 255 } 256 res+=flag; 257 } 258 return res; 259 } 260 //三维凸包重心 261 Point barycenter() 262 { 263 Point ans(0,0,0),o(0,0,0); 264 double all=0; 265 for(int i=0;i<num;i++) 266 { 267 double vol=volume(o,P[F[i].a],P[F[i].b],P[F[i].c]); 268 ans=ans+(o+P[F[i].a]+P[F[i].b]+P[F[i].c])/4.0*vol; 269 all+=vol; 270 } 271 ans=ans/all; 272 return ans; 273 } 274 //点到面的距离 275 double ptoface(Point p,int i) 276 { 277 return fabs(volume(P[F[i].a],P[F[i].b],P[F[i].c],p)/vlen((P[F[i].b]-P[F[i].a])*(P[F[i].c]-P[F[i].a]))); 278 } 279 }; 280 CH3D hull; 281 int main() 282 { 283 // freopen("in.txt","r",stdin); 284 // freopen("out.txt","w",stdout); 285 while(scanf("%d",&hull.n)==1) 286 { 287 for(int i=0;i<hull.n;i++) 288 { 289 scanf("%lf%lf%lf",&hull.P[i].x,&hull.P[i].y,&hull.P[i].z); 290 } 291 hull.create(); 292 Point p=hull.barycenter(); 293 double minn=1e20; 294 for(int i=0;i<hull.num;i++) 295 { 296 minn=min(minn,hull.ptoface(p,i)); 297 } 298 printf("%.3lf ",minn); 299 } 300 return 0; 301 }