UVA12305 Polishing a Extruded Polygon 多面体切割 [Rujia Liu's Presents, Dedicated to Geometry and CG Lovers]
不知道这算不算个神题了,AC的时候只有我和汝哥两人过这题……
想了一天,调了一天,竟然因为三角剖分时就差了一句叉积判断相邻边的凹凸性WA了好几天……
1、三角剖分:
所给多面体是非凸的,难以处理,剖分成一个个三棱柱就都是凸多面体了。最后才开始写陌生的三角剖分的,已经写+调了二百多行疲惫不堪的时候看到一计算几何书上好复杂的nlogn算法,竟然还要再搞平衡树,简直要崩溃。百度、谷歌都不给力了,期刊论文看着也晕乎,只好凭着自己的理解来暴力剖分了。
读入每个点建立双向链表(单向应该也没关系,灵活性差点),然后开始循环剖分,对每个相邻边,先判断凹凸性,然后枚举其他所有边判是否和如图所示虚线交叉,通过判断后可以确定这是个凸出的角,可以切掉,就得到了一个三角形。不停地对剩下的多边形进行这样的操作,每一时刻一定会有这样的凸出的角可以切,最后可以完成剖分。
2、模拟切割:
剖分的同时就建立了一个个三棱柱,对“外侧”的面打上标记,这样才能计算正确的表面积。我用链表来存多面体的各个面,用构成多面体的各个面来表示这个多面体。切割的时候就是切这些面,思考起来挺麻烦的,实现的时候思路还是比较清晰,切割凸多面体得到的任何面一定还是凸多边形,而切凸多边形只需要把切面一侧的点去掉,与切面相交的两点连起来,和剩下的点围成新的凸包。对于切掉的整个面,乃至整个凸多面体,就可以有效地利用链表轻松删除。
3、计算体积和面积:
面积就用多边形面积公式,和二维坐标计算的原理是一样的。由于是凸多面体,那么凸多面体内任选一点,就可以把多面体分割成一个个棱锥,棱锥的体积就是点到面的距离乘以面积除以3。这个点可以随意取,取多面体的一个顶点就非常方便了。
调试虽然花了不少时间,但是工作量并不算大,基本上调出样例就差不多了。
应该是没有三点共线的数据,我把处理这个的代码删了是可以AC的。
好长的代码……
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include<stdio.h> 2 #include<string.h> 3 #include<stdlib.h> 4 #include<math.h> 5 #include<vector> 6 #include<list> 7 #include<algorithm> 8 #include<iostream> 9 using namespace std; 10 const int maxn = 3011; 11 const double eps = 1e-8; 12 inline int dcmp(double x) 13 { 14 return x > eps ? 1 : (x < -eps ? -1 : 0); 15 } 16 inline double pz(double x) 17 { 18 return dcmp(x) ? x : 0; 19 } 20 /******************************************************************************/ 21 //定义三维点 22 struct Point3 23 { 24 double x, y, z; 25 Point3() 26 { 27 x = y = z = 0; 28 } 29 Point3(double a, double b, double c) 30 { 31 x = a, y = b, z = c; 32 } 33 Point3 cross(Point3 p) 34 { 35 return Point3(y * p.z - p.y * z, 36 z * p.x - x * p.z, 37 x * p.y - y * p.x); 38 } 39 double dot(Point3 p) 40 { 41 return x * p.x + y * p.y + z * p.z; 42 } 43 Point3 operator-(const Point3 &p)const 44 { 45 return Point3(x - p.x, y - p.y, z - p.z); 46 } 47 Point3 operator-()const 48 { 49 return Point3(-x, -y, -z); 50 } 51 Point3 operator+(const Point3 &p)const 52 { 53 return Point3(x + p.x, y + p.y, z + p.z); 54 } 55 Point3 operator*(const double b)const 56 { 57 return Point3(x * b, y * b, z * b); 58 } 59 bool operator==(const Point3 &p)const 60 { 61 return !dcmp(x - p.x) && !dcmp(y - p.y) && !dcmp(z - p.z); 62 } 63 bool operator!=(const Point3 &p)const 64 { 65 return !(*this == p); 66 } 67 bool operator<(const Point3 &p)const 68 { 69 if(!dcmp(x - p.x)) 70 { 71 if(!dcmp(y - p.y)) 72 return dcmp(z - p.z) < 0; 73 return dcmp(y - p.y) < 0; 74 } 75 return dcmp(x - p.x) < 0; 76 } 77 bool operator>(const Point3 &p)const 78 { 79 return p < *this; 80 } 81 bool operator>=(const Point3 &p)const 82 { 83 return !(*this < p); 84 } 85 bool operator<=(const Point3 &p)const 86 { 87 return !(*this > p); 88 } 89 Point3 fxl(Point3 b, Point3 c)//平面法向量 90 { 91 return (*this - b).cross(b - c); 92 } 93 double Dis(Point3 b) 94 { 95 return sqrt((*this - b).dot(*this - b)); 96 } 97 double vlen() 98 { 99 return sqrt(dot(*this)); 100 } 101 bool PinLine(Point3 b, Point3 c)//三点共线 102 { 103 return fxl(b, c).vlen() < eps; 104 } 105 bool PonPlane(Point3 b, Point3 c, Point3 d)//四点共面 106 { 107 return !dcmp(fxl(b, c).dot(d - *this)); 108 } 109 bool PonSeg(Point3 b, Point3 c)//点在线段上,包括端点 110 { 111 return !dcmp((*this - b).cross(*this - c).vlen()) && 112 (*this - b).dot(*this - c) <= 0; 113 } 114 bool PinSeg(Point3 b, Point3 c)//点在线段上,不包括端点 115 { 116 return PonSeg(b, c) && *this != b && *this != c; 117 } 118 double PtoLine(Point3 b, Point3 c)//点到直线距离 119 { 120 return (*this - b).cross(c - b).vlen() / b.Dis(c); 121 } 122 double PtoPlane(Point3 b, Point3 c, Point3 d)//点到平面距离 123 { 124 return fabs(b.fxl(c, d).dot(*this - b)) / b.fxl(c, d).vlen(); 125 } 126 }; 127 /******************************************************************************/ 128 //定义平面+空间平面凸包 129 struct Plane 130 { 131 double a, b, c, d; 132 bool outplane;//计入表面积的面 133 vector<Point3> p; 134 Plane() 135 { 136 a = b = c = d = 0, outplane = false; 137 p.clear(); 138 } 139 inline void init(double a_, double b_, double c_, double d_) 140 { 141 a = a_, b = b_, c = c_, d = d_; 142 p.clear(); 143 } 144 inline void init(Point3 pa, Point3 pb, Point3 pc) 145 { 146 Point3 t = (pa - pb).cross(pa - pc); 147 a = t.x, b = t.y, c = t.z; 148 d = -(pa.x * t.x + pa.y * t.y + pa.z * t.z); 149 p.clear(); 150 } 151 Plane(double a_, double b_, double c_, double d_) 152 { 153 init(a_, b_, c_, d_); 154 } 155 Plane(Point3 pa, Point3 pb, Point3 pc) 156 { 157 init(pa, pb, pc); 158 } 159 double PtoPlane(const Point3 &pa)const//点面距离 160 { 161 return fabs(Sub(pa)) / sqrt(a * a + b * b + c * c); 162 } 163 double Sub(const Point3 &p)const//点代入方程 164 { 165 return p.x * a + p.y * b + p.z * c + d; 166 } 167 Point3 PcrossPlane(Point3 a, Point3 b) 168 { 169 return a + (b - a) * (PtoPlane(a) / (PtoPlane(a) + PtoPlane(b))); 170 } 171 int Parallel(Plane pl)//面平行 172 { 173 if(!dcmp(a * pl.b - b * pl.a) && !dcmp(a * pl.c - c * pl.a) && !dcmp(b * pl.c - c * pl.b)) 174 return dcmp(pl.Sub(p[0])) > 0 ? 1 : -1; 175 return 0; 176 } 177 bool Cut(Plane &pl)//平面切割 178 { 179 switch(Parallel(pl)) 180 { 181 case -1: 182 return true; 183 case 1: 184 return false; 185 } 186 int i, j, k, n = p.size(); 187 bool flag1, flag2; 188 Point3 p1, p2; 189 vector<Point3> ret; 190 for(i = flag1 = flag2 = 0; i < n; ++ i) 191 { 192 flag1 |= pl.Sub(p[i]) < 0; 193 flag2 |= pl.Sub(p[i]) > 0; 194 } 195 if(flag1 != flag2) return flag1; 196 if(!flag1) return true; 197 for(i = 0; pl.Sub(p[i]) >= 0; ++ i); 198 for(; pl.Sub(p[i]) < 0; i = (i + 1) % n); 199 for(j = i; pl.Sub(p[j]) >= 0; j = (j + 1) % n); 200 for(k = j; k != i; k = (k + 1) % n) 201 ret.push_back(p[k]); 202 p1 = pl.PcrossPlane(p[i], p[(i + n - 1) % n]); 203 p2 = pl.PcrossPlane(p[j], p[(j + n - 1) % n]); 204 ret.push_back(p1), ret.push_back(p2); 205 p = ret; 206 pl.p.push_back(p1), pl.p.push_back(p2); 207 return true; 208 } 209 void Graham()//求空间平面凸包 210 { 211 int len, i, n, top = 1; 212 vector<Point3> res; 213 Point3 fx(a, b, c); 214 sort(p.begin(), p.end()); 215 vector<Point3>::iterator iter = unique(p.begin(), p.end()); 216 p.erase(iter, p.end()); 217 n = p.size(); 218 res.push_back(p[0]), res.push_back(p[1]); 219 for(i = 2; i < n; ++ i) 220 { 221 while(top && dcmp((res[top] - res[top - 1]).cross(p[i] - res[top]).dot(fx)) <= 0) 222 res.pop_back(), -- top; 223 res.push_back(p[i]), ++ top; 224 } 225 len = top; 226 res.push_back(p[n - 2]), ++ top; 227 for(i = n - 3; i >= 0; -- i) 228 { 229 while(top != len && dcmp((res[top] - res[top - 1]).cross(p[i] - res[top]).dot(fx)) <= 0) 230 res.pop_back(), -- top; 231 res.push_back(p[i]), ++ top; 232 } 233 res.pop_back(); 234 p = res; 235 } 236 double PolygonArea()//空间平面凸包面积 237 { 238 int n = p.size(); 239 double ret = 0; 240 for(int i = 2; i < n; ++ i) 241 ret += (p[i - 1] - p[0]).cross(p[i] - p[0]).vlen(); 242 return ret * 0.5; 243 } 244 }; 245 /******************************************************************************/ 246 //定义变量 247 struct Polyhedron//立方体面集合 248 { 249 list<Plane> pl; 250 }; 251 list<Polyhedron> pol; 252 int n, h, m;//点个数,高度,切割次数 253 struct PointChain 254 { 255 Point3 p; 256 int pre, nex; 257 bool outside;//标记发向nex的边为外部边,抬起的面计入总面积 258 } PC[maxn]; 259 /******************************************************************************/ 260 //判相交,辅助判断三角剖分合法性 261 inline double det(double x1, double y1, double x2, double y2) 262 { 263 return x1 * y2 - x2 * y1; 264 } 265 inline double cross2(Point3 a, Point3 b, Point3 c)//平面叉积 266 { 267 return det(b.x - a.x, b.y - a.y, c.x - a.x, c.y - a.y); 268 } 269 inline bool SegCross(Point3 u1, Point3 u2, Point3 v1, Point3 v2) 270 { 271 return ((dcmp(cross2(u1, u2, v1)) ^ dcmp(cross2(u1, u2, v2))) == -2 && 272 (dcmp(cross2(v1, v2, u1)) ^ dcmp(cross2(v1, v2, u2))) == -2) || 273 ((v1.PinSeg(u1, u2) || v2.PinSeg(u1, u2)) && 274 (u1.PonSeg(v1, v2) || u2.PonSeg(v1, v2))) || 275 ((v1.PonSeg(u1, u2) || v2.PonSeg(u1, u2)) && 276 (u1.PinSeg(v1, v2) || u2.PinSeg(v1, v2))); 277 } 278 /******************************************************************************/ 279 bool CanDo(int s, int e) 280 { 281 int tp = PC[s].nex; 282 while(tp != s) 283 { 284 if(SegCross(PC[s].p, PC[e].p, PC[tp].p, PC[PC[tp].pre].p)) 285 return false; 286 tp = PC[tp].nex; 287 } 288 return true; 289 } 290 void AddPolyhedron(int A, int B, int C) 291 { 292 Polyhedron t; 293 Plane pl; 294 Point3 a = PC[A].p, b = PC[B].p, c = PC[C].p; 295 Point3 a_ = Point3(a.x, a.y, h), b_ = Point3(b.x, b.y, h), c_ = Point3(c.x, c.y, h); 296 pl.init(a, b, c); 297 pl.p.push_back(a); 298 pl.p.push_back(b); 299 pl.p.push_back(c); 300 pl.outplane = true; 301 t.pl.push_front(pl); 302 pl.init(a_, b_, c_); 303 pl.p.push_back(a_); 304 pl.p.push_back(b_); 305 pl.p.push_back(c_); 306 t.pl.push_front(pl); 307 pl.init(a, b, b_); 308 pl.p.push_back(a); 309 pl.p.push_back(b); 310 pl.p.push_back(b_); 311 pl.p.push_back(a_); 312 pl.outplane = PC[A].outside, PC[A].outside = false; 313 t.pl.push_front(pl); 314 pl.init(b, c, c_); 315 pl.p.push_back(b); 316 pl.p.push_back(c); 317 pl.p.push_back(c_); 318 pl.p.push_back(b_); 319 pl.outplane = PC[B].outside, PC[B].outside = false; 320 t.pl.push_front(pl); 321 pl.init(c, a, a_); 322 pl.p.push_back(c); 323 pl.p.push_back(a); 324 pl.p.push_back(a_); 325 pl.p.push_back(c_); 326 if(PC[C].nex == A && PC[C].outside) 327 pl.outplane = true, PC[C].outside = false; 328 else pl.outplane = false; 329 t.pl.push_front(pl); 330 pol.push_front(t); 331 } 332 void Triangulation()//循环枚举相邻三点划分三角形 333 { 334 int tp = 0; 335 while(PC[PC[tp].nex].nex != tp) 336 { 337 while(dcmp(cross2(PC[tp].p, PC[PC[tp].nex].p, PC[PC[PC[tp].nex].nex].p)) <= 0 || 338 !CanDo(tp, PC[PC[tp].nex].nex)) tp = PC[tp].nex; 339 AddPolyhedron(tp, PC[tp].nex, PC[PC[tp].nex].nex); 340 PC[PC[PC[tp].nex].nex].pre = tp; 341 PC[tp].nex = PC[PC[tp].nex].nex; 342 } 343 } 344 void init() 345 { 346 double x, y; 347 int i, tp; 348 pol.clear(); 349 for(i = 0; i < n; ++ i) 350 { 351 scanf("%lf%lf", &x, &y); 352 PC[i].p = Point3(x, y, 0); 353 PC[i].pre = i - 1; 354 PC[i].nex = i + 1; 355 PC[i].outside = true; 356 } 357 PC[0].pre = n - 1, PC[n - 1].nex = 0; 358 Triangulation(); 359 } 360 /******************************************************************************/ 361 list<Polyhedron>::iterator iti; 362 list<Plane>::iterator itj; 363 double CalSumVol()//求总体积 364 { 365 double ans = 0; 366 Point3 tmp; 367 for(iti = pol.begin(); iti != pol.end(); ++ iti) 368 { 369 for(itj = (*iti).pl.begin(), tmp = (*itj).p[0]; itj != (*iti).pl.end(); ++ itj) 370 ans += (*itj).PolygonArea() * (*itj).PtoPlane(tmp); 371 } 372 return ans / 3; 373 } 374 double CalSumArea()//求总面积 375 { 376 double ans = 0; 377 for(iti = pol.begin(); iti != pol.end(); ++ iti) 378 for(itj = (*iti).pl.begin(); itj != (*iti).pl.end(); ++ itj) 379 if((*itj).outplane) ans += (*itj).PolygonArea(); 380 return ans; 381 } 382 383 void MakeAns() 384 { 385 double a, b, c, d; 386 bool flag; 387 while(m --) 388 { 389 scanf("%lf%lf%lf%lf", &a, &b, &c, &d); 390 for(iti = pol.begin(); iti != pol.end();) 391 { 392 Plane tmp(a, b, c, d); 393 tmp.outplane = true; 394 for(itj = (*iti).pl.begin(), flag = false; itj != (*iti).pl.end();) 395 { 396 if(!(*itj).Cut(tmp)) 397 { 398 list<Plane>::iterator tmpj = itj; 399 ++ itj; 400 (*iti).pl.erase(tmpj); 401 } 402 else flag = true, ++ itj; 403 } 404 if(!flag) 405 { 406 list<Polyhedron>::iterator tmpi = iti; 407 ++ iti; 408 pol.erase(tmpi); 409 } 410 else if(tmp.p.size()) 411 { 412 tmp.Graham(); 413 (*iti).pl.push_front(tmp); 414 ++ iti; 415 } 416 else ++ iti; 417 } 418 printf("%.3f %.3f\n", CalSumVol(), CalSumArea()); 419 } 420 } 421 int main() 422 { 423 while(scanf("%d%d%d", &n, &h, &m), n | h | m) 424 { 425 init(); 426 printf("%.3f %.3f\n", CalSumVol(), CalSumArea()); 427 MakeAns(); 428 } 429 return 0; 430 }