卡了一个月啊，都快想放弃了。今天看了一下LRJ的PSLG版的代码，才想起来我的模拟捆包裹法存在的漏洞。那就是对点的坐标进行排序的时候，应该用的sgn()函数来进行比较，而不是我一开始的直接对浮点数进行比较。如果直接对浮点数进行比较，本来距离差不多的点就会被拉伸开来，导致应该是相邻点的点被分隔了。
做法比较简单，我的捆包裹法是先将所有的交点离散出来，然后处理出所有的相邻的点。最后，就是模拟向顺时针转动最多的方向的相邻点移动的所有过程。
代码如下：
1 #include <cstdio> 2 #include <cstring> 3 #include <cmath> 4 #include <set> 5 #include <vector> 6 #include <iostream> 7 #include <algorithm> 8 9 using namespace std; 10 11 // Point class 12 struct Point { 13 double x, y; 14 Point() {} 15 Point(double x, double y) : x(x), y(y) {} 16 } ; 17 template<class T> T sqr(T x) { return x * x;} 18 inline double ptDis(Point a, Point b) { return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y));} 19 20 // basic calculations 21 typedef Point Vec; 22 Vec operator + (Vec a, Vec b) { return Vec(a.x + b.x, a.y + b.y);} 23 Vec operator - (Vec a, Vec b) { return Vec(a.x - b.x, a.y - b.y);} 24 Vec operator * (Vec a, double p) { return Vec(a.x * p, a.y * p);} 25 Vec operator / (Vec a, double p) { return Vec(a.x / p, a.y / p);} 26 27 const double EPS = 1e-6; 28 const double PI = acos(-1.0); 29 inline int sgn(double x) { return fabs(x) < EPS ? 0 : (x < 0 ? -1 : 1);} 30 bool operator < (Point a, Point b) { return sgn(a.x - b.x) < 0 || (sgn(a.x - b.x) == 0 && a.y < b.y);} 31 bool operator == (Point a, Point b) { return sgn(a.x - b.x) == 0 && sgn(a.y - b.y) == 0;} 32 33 inline double dotDet(Vec a, Vec b) { return a.x * b.x + a.y * b.y;} 34 inline double crossDet(Vec a, Vec b) { return a.x * b.y - a.y * b.x;} 35 inline double crossDet(Point o, Point a, Point b) { return crossDet(a - o, b - o);} 36 inline double vecLen(Vec x) { return sqrt(sqr(x.x) + sqr(x.y));} 37 inline double toRad(double deg) { return deg / 180.0 * PI;} 38 inline double angle(Vec v) { return atan2(v.y, v.x);} 39 inline double angle(Vec a, Vec b) { return acos(dotDet(a, b) / vecLen(a) / vecLen(b));} 40 inline double triArea(Point a, Point b, Point c) { return fabs(crossDet(b - a, c - a));} 41 inline Vec vecUnit(Vec x) { return x / vecLen(x);} 42 inline Vec rotate(Vec x, double rad) { return Vec(x.x * cos(rad) - x.y * sin(rad), x.x * sin(rad) + x.y * cos(rad));} 43 Vec normal(Vec x) { 44 double len = vecLen(x); 45 return Vec(- x.y / len, x.x / len); 46 } 47 48 // Line class 49 struct Line { 50 Point s, t; 51 Line() {} 52 Line(Point s, Point t) : s(s), t(t) {} 53 Point point(double x) { 54 return s + (t - s) * x; 55 } 56 Line move(double x) { // while x > 0 move to (s->t)'s left 57 Vec nor = normal(t - s); 58 nor = nor * x; 59 return Line(s + nor, t + nor); 60 } 61 Vec vec() { return t - s;} 62 } ; 63 typedef Line Seg; 64 65 inline bool onSeg(Point x, Point a, Point b) { return sgn(crossDet(a - x, b - x)) == 0 && sgn(dotDet(a - x, b - x)) < 0;} 66 inline bool onSeg(Point x, Seg s) { return onSeg(x, s.s, s.t);} 67 68 // 0 : not intersect 69 // 1 : proper intersect 70 // 2 : improper intersect 71 int segIntersect(Point a, Point c, Point b, Point d) { 72 Vec v1 = b - a, v2 = c - b, v3 = d - c, v4 = a - d; 73 int a_bc = sgn(crossDet(v1, v2)); 74 int b_cd = sgn(crossDet(v2, v3)); 75 int c_da = sgn(crossDet(v3, v4)); 76 int d_ab = sgn(crossDet(v4, v1)); 77 if (a_bc * c_da > 0 && b_cd * d_ab > 0) return 1; 78 if (onSeg(b, a, c) && c_da) return 2; 79 if (onSeg(c, b, d) && d_ab) return 2; 80 if (onSeg(d, c, a) && a_bc) return 2; 81 if (onSeg(a, d, b) && b_cd) return 2; 82 return 0; 83 } 84 inline int segIntersect(Seg a, Seg b) { return segIntersect(a.s, a.t, b.s, b.t);} 85 86 // point of the intersection of 2 lines 87 Point lineIntersect(Point P, Vec v, Point Q, Vec w) { 88 Vec u = P - Q; 89 double t = crossDet(w, u) / crossDet(v, w); 90 return P + v * t; 91 } 92 inline Point lineIntersect(Line a, Line b) { return lineIntersect(a.s, a.t - a.s, b.s, b.t - b.s);} 93 94 // Warning: This is a DIRECTED Distance!!! 95 double pt2Line(Point x, Point a, Point b) { 96 Vec v1 = b - a, v2 = x - a; 97 return crossDet(v1, v2) / vecLen(v1); 98 } 99 inline double pt2Line(Point x, Line L) { return pt2Line(x, L.s, L.t);} 100 101 double pt2Seg(Point x, Point a, Point b) { 102 if (a == b) return vecLen(x - a); 103 Vec v1 = b - a, v2 = x - a, v3 = x - b; 104 if (sgn(dotDet(v1, v2)) < 0) return vecLen(v2); 105 if (sgn(dotDet(v1, v3)) > 0) return vecLen(v3); 106 return fabs(crossDet(v1, v2)) / vecLen(v1); 107 } 108 inline double pt2Seg(Point x, Seg s) { return pt2Seg(x, s.s, s.t);} 109 110 // Polygon class 111 struct Poly { 112 vector<Point> pt; 113 Poly() { pt.clear();} 114 Poly(vector<Point> pt) : pt(pt) {} 115 Point operator [] (int x) const { return pt[x];} 116 int size() { return pt.size();} 117 double area() { 118 double ret = 0.0; 119 int sz = pt.size(); 120 for (int i = 1; i < sz; i++) { 121 ret += crossDet(pt[i], pt[i - 1]); 122 } 123 return fabs(ret / 2.0); 124 } 125 } ; 126 127 // Circle class 128 struct Circle { 129 Point c; 130 double r; 131 Circle() {} 132 Circle(Point c, double r) : c(c), r(r) {} 133 Point point(double a) { 134 return Point(c.x + cos(a) * r, c.y + sin(a) * r); 135 } 136 } ; 137 138 // Cirlce operations 139 int lineCircleIntersect(Line L, Circle C, double &t1, double &t2, vector<Point> &sol) { 140 double a = L.s.x, b = L.t.x - C.c.x, c = L.s.y, d = L.t.y - C.c.y; 141 double e = sqr(a) + sqr(c), f = 2 * (a * b + c * d), g = sqr(b) + sqr(d) - sqr(C.r); 142 double delta = sqr(f) - 4.0 * e * g; 143 if (sgn(delta) < 0) return 0; 144 if (sgn(delta) == 0) { 145 t1 = t2 = -f / (2.0 * e); 146 sol.push_back(L.point(t1)); 147 return 1; 148 } 149 t1 = (-f - sqrt(delta)) / (2.0 * e); 150 sol.push_back(L.point(t1)); 151 t2 = (-f + sqrt(delta)) / (2.0 * e); 152 sol.push_back(L.point(t2)); 153 return 2; 154 } 155 156 int lineCircleIntersect(Line L, Circle C, vector<Point> &sol) { 157 Vec dir = L.t - L.s, nor = normal(dir); 158 Point mid = lineIntersect(C.c, nor, L.s, dir); 159 double len = sqr(C.r) - sqr(ptDis(C.c, mid)); 160 if (sgn(len) < 0) return 0; 161 if (sgn(len) == 0) { 162 sol.push_back(mid); 163 return 1; 164 } 165 Vec dis = vecUnit(dir); 166 len = sqrt(len); 167 sol.push_back(mid + dis * len); 168 sol.push_back(mid - dis * len); 169 return 2; 170 } 171 172 // -1 : coincide 173 int circleCircleIntersect(Circle C1, Circle C2, vector<Point> &sol) { 174 double d = vecLen(C1.c - C2.c); 175 if (sgn(d) == 0) { 176 if (sgn(C1.r - C2.r) == 0) { 177 return -1; 178 } 179 return 0; 180 } 181 if (sgn(C1.r + C2.r - d) < 0) return 0; 182 if (sgn(fabs(C1.r - C2.r) - d) > 0) return 0; 183 double a = angle(C2.c - C1.c); 184 double da = acos((sqr(C1.r) + sqr(d) - sqr(C2.r)) / (2.0 * C1.r * d)); 185 Point p1 = C1.point(a - da), p2 = C1.point(a + da); 186 sol.push_back(p1); 187 if (p1 == p2) return 1; 188 sol.push_back(p2); 189 return 2; 190 } 191 192 void circleCircleIntersect(Circle C1, Circle C2, vector<double> &sol) { 193 double d = vecLen(C1.c - C2.c); 194 if (sgn(d) == 0) return ; 195 if (sgn(C1.r + C2.r - d) < 0) return ; 196 if (sgn(fabs(C1.r - C2.r) - d) > 0) return ; 197 double a = angle(C2.c - C1.c); 198 double da = acos((sqr(C1.r) + sqr(d) - sqr(C2.r)) / (2.0 * C1.r * d)); 199 sol.push_back(a - da); 200 sol.push_back(a + da); 201 } 202 203 int tangent(Point p, Circle C, vector<Vec> &sol) { 204 Vec u = C.c - p; 205 double dist = vecLen(u); 206 if (dist < C.r) return 0; 207 if (sgn(dist - C.r) == 0) { 208 sol.push_back(rotate(u, PI / 2.0)); 209 return 1; 210 } 211 double ang = asin(C.r / dist); 212 sol.push_back(rotate(u, -ang)); 213 sol.push_back(rotate(u, ang)); 214 return 2; 215 } 216 217 // ptA : points of tangency on circle A 218 // ptB : points of tangency on circle B 219 int tangent(Circle A, Circle B, vector<Point> &ptA, vector<Point> &ptB) { 220 if (A.r < B.r) { 221 swap(A, B); 222 swap(ptA, ptB); 223 } 224 double d2 = sqr(A.c.x - B.c.x) + sqr(A.c.y - B.c.y); 225 double rdiff = A.r - B.r, rsum = A.r + B.r; 226 if (d2 < sqr(rdiff)) return 0; 227 double base = atan2(B.c.y - A.c.y, B.c.x - A.c.x); 228 if (d2 == 0 && A.r == B.r) return -1; 229 if (d2 == sqr(rdiff)) { 230 ptA.push_back(A.point(base)); 231 ptB.push_back(B.point(base)); 232 return 1; 233 } 234 double ang = acos((A.r - B.r) / sqrt(d2)); 235 ptA.push_back(A.point(base + ang)); 236 ptB.push_back(B.point(base + ang)); 237 ptA.push_back(A.point(base - ang)); 238 ptB.push_back(B.point(base - ang)); 239 if (d2 == sqr(rsum)) { 240 ptA.push_back(A.point(base)); 241 ptB.push_back(B.point(PI + base)); 242 } else if (d2 > sqr(rsum)) { 243 ang = acos((A.r + B.r) / sqrt(d2)); 244 ptA.push_back(A.point(base + ang)); 245 ptB.push_back(B.point(PI + base + ang)); 246 ptA.push_back(A.point(base - ang)); 247 ptB.push_back(B.point(PI + base - ang)); 248 } 249 return (int) ptA.size(); 250 } 251 252 // -1 : onside 253 // 0 : outside 254 // 1 : inside 255 int ptInPoly(Point p, Poly &poly) { 256 int wn = 0, sz = poly.size(); 257 for (int i = 0; i < sz; i++) { 258 if (onSeg(p, poly[i], poly[(i + 1) % sz])) return -1; 259 int k = sgn(crossDet(poly[(i + 1) % sz] - poly[i], p - poly[i])); 260 int d1 = sgn(poly[i].y - p.y); 261 int d2 = sgn(poly[(i + 1) % sz].y - p.y); 262 if (k > 0 && d1 <= 0 && d2 > 0) wn++; 263 if (k < 0 && d2 <= 0 && d1 > 0) wn--; 264 } 265 if (wn != 0) return 1; 266 return 0; 267 } 268 269 // if DO NOT need a high precision 270 /* 271 int ptInPoly(Point p, Poly poly) { 272 int sz = poly.size(); 273 double ang = 0.0, tmp; 274 for (int i = 0; i < sz; i++) { 275 if (onSeg(p, poly[i], poly[(i + 1) % sz])) return -1; 276 tmp = angle(poly[i] - p) - angle(poly[(i + 1) % sz] - p) + PI; 277 ang += tmp - floor(tmp / (2.0 * PI)) * 2.0 * PI - PI; 278 } 279 if (sgn(ang - PI) == 0) return -1; 280 if (sgn(ang) == 0) return 0; 281 return 1; 282 } 283 */ 284 285 286 // Convex Hull algorithms 287 // return the number of points in convex hull 288 289 // andwer's algorithm 290 // if DO NOT want the points on the side of convex hull, change all "<" into "<=" 291 int andrew(Point *pt, int n, Point *ch) { 292 sort(pt, pt + n); 293 int m = 0; 294 for (int i = 0; i < n; i++) { 295 while (m > 1 && crossDet(ch[m - 1] - ch[m - 2], pt[i] - ch[m - 2]) <= 0) m--; 296 ch[m++] = pt[i]; 297 } 298 int k = m; 299 for (int i = n - 2; i >= 0; i--) { 300 while (m > k && crossDet(ch[m - 1] - ch[m - 2], pt[i] - ch[m - 2]) <= 0) m--; 301 ch[m++] = pt[i]; 302 } 303 if (n > 1) m--; 304 return m; 305 } 306 307 // graham's algorithm 308 // if DO NOT want the points on the side of convex hull, change all "<=" into "<" 309 Point origin; 310 inline bool cmpAng(Point p1, Point p2) { return crossDet(origin, p1, p2) > 0;} 311 inline bool cmpDis(Point p1, Point p2) { return ptDis(p1, origin) > ptDis(p2, origin);} 312 313 void removePt(Point *pt, int &n) { 314 int idx = 1; 315 for (int i = 2; i < n; i++) { 316 if (sgn(crossDet(origin, pt[i], pt[idx]))) pt[++idx] = pt[i]; 317 else if (cmpDis(pt[i], pt[idx])) pt[idx] = pt[i]; 318 } 319 n = idx + 1; 320 } 321 322 int graham(Point *pt, int n, Point *ch) { 323 int top = -1; 324 for (int i = 1; i < n; i++) { 325 if (pt[i].y < pt[0].y || (pt[i].y == pt[0].y && pt[i].x < pt[0].x)) swap(pt[i], pt[0]); 326 } 327 origin = pt[0]; 328 sort(pt + 1, pt + n, cmpAng); 329 removePt(pt, n); 330 for (int i = 0; i < n; i++) { 331 if (i >= 2) { 332 while (!(crossDet(ch[top - 1], pt[i], ch[top]) <= 0)) top--; 333 } 334 ch[++top] = pt[i]; 335 } 336 return top + 1; 337 } 338 339 340 // Half Plane 341 // The intersected area is always a convex polygon. 342 // Directed Line class 343 struct DLine { 344 Point p; 345 Vec v; 346 double ang; 347 DLine() {} 348 DLine(Point p, Vec v) : p(p), v(v) { ang = atan2(v.y, v.x);} 349 bool operator < (const DLine &L) const { return ang < L.ang;} 350 DLine move(double x) { // while x > 0 move to v's left 351 Vec nor = normal(v); 352 nor = nor * x; 353 return DLine(p + nor, v); 354 } 355 356 } ; 357 358 Poly cutPoly(Poly &poly, Point a, Point b) { 359 Poly ret = Poly(); 360 int n = poly.size(); 361 for (int i = 0; i < n; i++) { 362 Point c = poly[i], d = poly[(i + 1) % n]; 363 if (sgn(crossDet(b - a, c - a)) >= 0) ret.pt.push_back(c); 364 if (sgn(crossDet(b - a, c - d)) != 0) { 365 Point ip = lineIntersect(a, b - a, c, d - c); 366 if (onSeg(ip, c, d)) ret.pt.push_back(ip); 367 } 368 } 369 return ret; 370 } 371 inline Poly cutPoly(Poly &poly, DLine L) { return cutPoly(poly, L.p, L.p + L.v);} 372 373 inline bool onLeft(DLine L, Point p) { return crossDet(L.v, p - L.p) > 0;} 374 Point dLineIntersect(DLine a, DLine b) { 375 Vec u = a.p - b.p; 376 double t = crossDet(b.v, u) / crossDet(a.v, b.v); 377 return a.p + a.v * t; 378 } 379 380 Poly halfPlane(DLine *L, int n) { 381 Poly ret = Poly(); 382 sort(L, L + n); 383 int fi, la; 384 Point *p = new Point[n]; 385 DLine *q = new DLine[n]; 386 q[fi = la = 0] = L[0]; 387 for (int i = 1; i < n; i++) { 388 while (fi < la && !onLeft(L[i], p[la - 1])) la--; 389 while (fi < la && !onLeft(L[i], p[fi])) fi++; 390 q[++la] = L[i]; 391 if (fabs(crossDet(q[la].v, q[la - 1].v)) < EPS) { 392 la--; 393 if (onLeft(q[la], L[i].p)) q[la] = L[i]; 394 } 395 if (fi < la) p[la - 1] = dLineIntersect(q[la - 1], q[la]); 396 } 397 while (fi < la && !onLeft(q[fi], p[la - 1])) la--; 398 if (la - fi <= 1) return ret; 399 p[la] = dLineIntersect(q[la], q[fi]); 400 for (int i = fi; i <= la; i++) ret.pt.push_back(p[i]); 401 return ret; 402 } 403 404 405 // 3D Geometry 406 void getCoor(double R, double lat, double lng, double &x, double &y, double &z) { 407 lat = toRad(lat); 408 lng = toRad(lng); 409 x = R * cos(lat) * cos(lng); 410 y = R * cos(lat) * sin(lng); 411 z = R * sin(lat); 412 } 413 414 415 /****************** template above *******************/ 416 417 const int N = 111; 418 Point ori[N]; 419 vector<Point> ex; 420 //vector<int> nx[5555], on[N]; 421 int nx[5555][N << 1], on[N][N << 1]; 422 423 int main() { 424 // freopen("in", "r", stdin); 425 // freopen("out", "w", stdout); 426 int n; 427 while (cin >> n) { 428 ex.clear(); 429 for (int i = 0; i < 5555; i++) 430 // nx[i].clear(); 431 nx[i][0] = 0; 432 for (int i = 0; i < n; i++) { 433 // on[i].clear(); 434 on[i][0] = 0; 435 cin >> ori[i].x >> ori[i].y; 436 ex.push_back(ori[i]); 437 } 438 ori[n] = ori[0]; 439 for (int i = 0; i < n; i++) { 440 for (int j = 0; j < n; j++) { 441 if (segIntersect(ori[i], ori[i + 1], ori[j], ori[j + 1]) == 1) { 442 ex.push_back(lineIntersect(Line(ori[i], ori[i + 1]), Line(ori[j], ori[j + 1]))); 443 } 444 } 445 } 446 sort(ex.begin(), ex.end()); 447 int t = unique(ex.begin(), ex.end()) - ex.begin(); 448 while (ex.size() > t) ex.pop_back(); 449 int sz = ex.size(); 450 for (int i = 0; i < sz; i++) { 451 for (int j = 0; j < n; j++) { 452 if (ex[i] == ori[j] || ex[i] == ori[j + 1] || onSeg(ex[i], Seg(ori[j], ori[j + 1]))) 453 // on[j].push_back(i); 454 on[j][++on[j][0]] = i; 455 } 456 } 457 // cout << ex.size() << " ~~~" << endl; 458 // for (int i = 0; i < ex.size(); i++) { 459 // cout << i << " : " << ex[i].x << ' ' << ex[i].y << endl; 460 // } 461 // for (int i = 0; i < n; i++) { 462 // cout << "on Line " << i << endl; 463 // for (int j = 1; j <= on[i][0]; j++) { 464 // cout << ex[on[i][j]].x << ' ' << ex[on[i][j]].y << endl; 465 // } 466 // } 467 // cout << "Debug!!~~~~~~~~~~~~" << endl; 468 for (int i = 0; i < n; i++) { 469 // int csz = on[i].size() - 1; 470 int csz = on[i][0] - 1; 471 for (int j = 0; j < csz; j++) { 472 // nx[on[i][j]].push_back(on[i][j + 1]); 473 // nx[on[i][j + 1]].push_back(on[i][j]); 474 nx[on[i][j + 1]][++nx[on[i][j + 1]][0]] = on[i][j + 2]; 475 nx[on[i][j + 2]][++nx[on[i][j + 2]][0]] = on[i][j + 1]; 476 } 477 } 478 vector<int> path; 479 path.clear(); 480 path.push_back(0); 481 // int sz0 = nx[0].size(); 482 int sz0 = nx[0][0]; 483 // int cur = nx[0][0]; 484 int cur = nx[0][1]; 485 for (int i = 1; i < sz0; i++) { 486 // int id = nx[0][i]; 487 int id = nx[0][i + 1]; 488 if (sgn(crossDet(ex[cur] - ex[0], ex[id] - ex[0])) < 0) cur = id; 489 } 490 // cout << cur << ' ' << ex[cur].x << ' ' << ex[cur].y << endl; 491 int cnt = 500000; 492 while (cnt-- && cur) { 493 path.push_back(cur); 494 int l1 = path[path.size() - 1], l0 = path[path.size() - 2]; 495 // int t = (nx[cur][0] == l0) ? nx[cur][1] : nx[cur][0], csz = nx[cur].size(); 496 int t = (nx[cur][1] == l0) ? nx[cur][2] : nx[cur][1], csz = nx[cur][0]; 497 Vec dr = ex[l1] - ex[l0]; 498 double ang = angle(dr); 499 for (int i = 1; i < csz; i++) { 500 // int id = nx[cur][i]; 501 int id = nx[cur][i + 1]; 502 if (id == l0) continue; 503 Vec t1 = ex[t] - ex[l1], t2 = ex[id] - ex[l1]; 504 t1 = rotate(t1, -ang), t2 = rotate(t2, -ang); 505 if (angle(t1) > angle(t2)) t = id; 506 } 507 cur = t; 508 // cout << ex[cur].x << ' ' << ex[cur].y << endl; 509 } 510 // cout << path.size() << endl; 511 int psz = path.size(); 512 vector<int> out; 513 out.clear(); 514 path.push_back(path[0]); 515 for (int i = 0; i < psz; i++) { 516 if (i > 0 && sgn(crossDet(ex[path[i]] - ex[path[i - 1]], ex[path[i + 1]] - ex[path[i]])) == 0) continue; 517 // cout << ex[path[i]].x << ' ' << ex[path[i]].y << endl; 518 out.push_back(path[i]); 519 } 520 cout << out.size() << endl; 521 psz = out.size(); 522 for (int i = 0; i < psz; i++) { 523 // cout << "!! " << out[i] << endl; 524 cout << ex[out[i]].x << ' ' << ex[out[i]].y << endl; 525 // printf("%.5f %.5f\n", ex[out[i]].x, ex[out[i]].y); 526 } 527 } 528 return 0; 529 }
PSLG法在LRJ的代码里面有，不放上来了。
——written by Lyon
621
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
