UVA 1340-find the Border (roll wrap)

Topic Link: UVA 1340-find the Border

A volume wrapping algorithm that mimics the code written by others. is to walk along the perimeter of the polyline, each time to a new point need to turn to the right to turn the most.

#include <cstdio> #include <cstring> #include <cmath> #include <vector> #include <complex > #include <algorithm>using namespace std;typedef pair<int,int> pii;const Double pi = 4 * atan (1); Const Doub Le EPS = 1e-8;inline int dcmp (double x) {if (Fabs (x) < EPS) return 0; else return x < 0? -1:1;} Inline double getdistance (double x, double y) {return sqrt (x * x + y * y);} Inline double Torad (double deg) {return deg/180 * PI;} struct point {double x, y; Point (Double x = 0, double y = 0): x (x), Y (y) {}void read () {scanf ("%lf%lf", &x, &y);} void Write () {printf ("%lf%lf", x, y);} BOOL operator = = (CONST point& u) Const {return dcmp (x-u.x) = = 0 && dcmp (y-u.y) = = 0;} BOOL operator! = (const point& u) Const {return! ( *this = = u); }bool operator < (const point& u) Const {return dcmp (x-u.x) < 0 | | (DCMP (x-u.x) ==0 && dcmp (Y-U.Y) < 0); }bool operator > (const point& u) Const {return u < *this; }bool operator <= (const point& u) Const {return *this < u | | *this = u;} BOOL operator >= (const point& u) const {return *this > U | | *this = u;} Point operator + (const point& u) {return point (x + u.x, y + u.y);} Point operator-(const point& u) {return point (x-u.x, y-u.y);} Point operator * (const double u) {return point (x * u, y * u);} Point operator/(const double u) {return point (x/u, y/u);} Double operator * (const point& u) {return x*u.y-y*u.x;}}; typedef point VECTOR;TYPEDEF vector<point> polygon;struct line {Double A, B, C; Line (Double A = 0, double b = 0, double c = 0): A (a), B (b), C (c) {}};struct dirline {point P; Vector v;double Ang;dirline () {}dirline (point P, Vector v): P (P), V (v) {ang = atan2 (V.y, v.x);} BOOL operator < (const dirline& u) const {return Ang < U.ang;}}; struct Circle {point o;double R; Circle () {}circle (Point O, double r = 0): O (O), R (r) {}void read () {o.read (), scanf ("%lf", &r); }point point (Double rad) {return point (o.x + cos (rad) *r, O.y + sin (rad) *r);} Double Getarea (double rad) {return rad * R * R/2;}}; Namespace punctual {Double getdistance (point A, point B) {double x=a.x-b.x, y=a.y-b.y; return sqrt (x*x + y*y);}};  namespace Vectorial {/* dot product: The product of two vector lengths multiplied by the cosine of their angle, the angle is greater than 90 degrees the product is negative */double Getdot (vector A, vector b) {return a.x * b.x + A.Y * B.Y;  }/* Cross product: The cross product is equal to two vectors consisting of a triangle with a direction area of twice times, crosses (V, w) =-cross (W, v) */double Getcross (vector A, vector b) {return a.x * B.Y-A.Y * b.x; }double GetLength (Vector a) {return sqrt (Getdot (A, a));} Double Getplength (Vector a) {return Getdot (A, a);} Double Getangle (Vector u) {return atan2 (u.y, u.x);} Double Getangle (vector A, vector b) {return ACOs (Getdot (A, B)/GetLength (a)/GetLength (b));} Vector Rotate (vector A, double rad) {return vector (A.x*cos (RAD)-a.y*sin (RAD), A.x*sin (RAD) +a.y*cos (RAD));} /* unit normal */vector Getnormal (Vector a) {double L = getlength (a); return Vector (-a.y/l, a.x/l);}}; NamesPace Complexvector {typedef complex<double> POINT;TYPEDEF point vector;double getdot (vector A, vector b) {return re Al (Conj (a) *b); }double Getcross (vector A, vector b) {return imag (Conj (a) *b);} Vector Rotate (vector A, double rad) {return a*exp (point (0, RAD));}}; Namespace Linear {using namespace vectorial; Line GetLine (double x1, double y1, double x2, double y2) {return line (y2-y1, x1-x2, y1*x2-x1*y2);} Line GetLine (double A, double b, point u) {return line (A, B, u.y * b-u.x * a);} BOOL Getintersection (line p, line Q, point& o) {if (Fabs (P.A * q.b-q.a * p.b) < EPS) return false;o.x = (Q.C * p. B-P.C * q.b)/(P.A * q.b-q.a * p.b); o.y = (Q.C * p.a-p.c * q.a)/(P.B * q.a-q.b * p.a); return true;} /* Intersection of linear PV and line QW */bool getintersection (point P, Vector v, point Q, Vector W, point& o) {if (dcmp (Getcross (V, w)) = = 0 ) return false; Vector u = p-q;double k = Getcross (w, u)/Getcross (V, w); o = p + v * K;return true;} /* point P to line ab distance */double getdistancetoline (pOint p, point A, point B) {return fabs (Getcross (b-a, p-a)/GetLength (B-A));} Double getdistancetosegment (point P, point A, point B) {if (a = = b) return getlength (P-A); Vector v1 = b-a, v2 = p-a, V3 = p-b;if (dcmp (Getdot, V1) < 0) return v2 (getlength); else if (V2 (dcmp (Getdot, V 3)) > 0) return GetLength (V3), Else return fabs (Getcross (v1, v2)/GetLength (v1));} /* point P on the line ab projection */point getpointtoline (points P, a, point B) {Vector v = b-a; return a+v* (Getdot (V, p-a)/Getdot (V, v)); }/* to determine if a segment has an intersection */bool haveintersection (Point A1, Point A2, point B1, point B2) {double C1=getcross (a2-a1, b1-a1), c2=getc Ross (A2-A1, B2-A1), C3=getcross (B2-B1, A1-B1), C4=getcross (B2-B1,A2-B1), return dcmp (C1) *dcmp (C2) < 0 && dcmp (C3) *dcmp (C4) < 0;} /* Determine if the point is */bool onsegment on a line segment {return dcmp (Getcross (A-p, b-p)) = = 0 && dcmp (Getdot (A- P, b-p)) < 0; }bool Onleft (dirline L, point P) {return dcmp (L.V * (P-L.P)) > 0;}} Namespace Triangular {Using namespace Vectorial;double getangle (double A, double b, double c) {return ACOs ((a*a+b*b-c*c)/(2*a*b));} Double Getarea (double A, double b, double c) {double s = (a+b+c)/2; return sqrt (s* (s-a) * (s-b) * (S-C));} Double Getarea (Double A, double h) {return a * H/2;} Double Getarea (Point A, point B, point C) {return fabs (Getcross (b-a, c-a))/2;} Double Getdirarea (Point A, point B, point C) {return Getcross (b-a, c-a)/2;}};  Namespace polygonal {using namespace vectorial;using namespace linear;double getarea (point* p, int n) {double ret = 0;for (int i = 1; i < n-1; i++) ret + = Getcross (P[i]-p[0], p[i+1]-p[0]); return fabs (ret)/2;} /* Convex package */int getconvexhull (point* p, int n, point* ch) {sort (p, p + N); int m = 0;for (int i = 0; i < n; i++) {/* collinear */ while (M > 1 && dcmp (Getcross (ch[m-1]-ch[m-2), p[i]-ch[m-1]) < 0) M--;while (M > 1 && dcmp (get Cross (Ch[m-1]-ch[m-2], p[i]-ch[m-1])) <= 0) m--;ch[m++] = P[i];} int k = m;for (int i = n-2; I>= 0; i--) {/* can be collinear *///while (M > K && dcmp (Getcross (Ch[m-1]-ch[m-2], p[i]-ch[m-2]) < 0) M--;while (M > K &amp ;& dcmp (Getcross (Ch[m-1]-ch[m-2], p[i]-ch[m-2])) <= 0) m--;ch[m++] = P[i];} if (n > 1) m--;return m;} int Ispointinpolygon (Point O, point* p, int n) {int wn = 0;for (int i = 0; i < n; i++) {int J = (i + 1)% n;if (onsegm ENT (o, p[i], p[j])) return 0; the boundary int k = dcmp (Getcross (P[j]-p[i], o-p[i])); int D1 = DCMP (p[i].y-o.y); int d2 = DCMP (P[J].Y-O.Y); if (k > 0 &am p;& D1 <= 0 && d2 > 0) wn++;if (k < 0 && D2 <= 0 && d1 > 0) wn--;} Return WN? -1:1;}  /* Rotate Jam */void rotatingcalipers (point *p, int n, vector<pii>& sol) {sol.clear (); int j = 1; P[n] = p[0];for (int i = 0; I < n; i++) {while (Getcross (p[j+1]-p[i+1], p[i]-p[i+1]) > Getcross (p[j]-p[i+1], p[i]-p[i+1])) j = (j+1)% N;sol.push_back ( Make_pair (i, J)); Sol.push_back (Make_pair (i + 1, j + 1));}} /* Calculate half-plane intersection can be by increment method, O (n^2), initially set 4 infinite half-plane *//* Cut the polygon u with a straight a->b and return to the left.  May degenerate into a single point or segment */polygon Cutpolygon (Polygon u, point A, point B) {Polygon Ret;int n = u.size (); for (int i = 0; i < n; i++) {Point c = u[i], d = u[(i+1)%n];if (dcmp ((b-a) * (c-a)) >= 0) Ret.push_back (c); if (dcmp ((b-a) * (c-d))! = 0) {Point T;geti Ntersection (A, b-a, C, d-c, T); if (Onsegment (t, C, D)) Ret.push_back (t);}} return ret;} /* Half plane intersects */int halfplaneintersection (dirline* li, int n, point* poly) {sort (Li, Li + N); int first, last; point* p = new Point[n];D irline* q = new Dirline[n];q[first=last=0] = li[0];for (int i = 1; i < n; i++) {while (First & Lt Last &&!onleft (Li[i], p[last-1]) Last--;while (First < last &&!onleft (Li[i], P[first])) first++;q[+ +last] = li[i];if (dcmp (Q[LAST].V * q[last-1].v) = = 0) {last--;if (Onleft], q[last)) LI[I].P] = Q[last];} if (first < last) getintersection (Q[LAST-1].P, Q[LAST-1].V, Q[LAST].P, Q[LAST].V, p[last-1]);} while (first < last &&!onleft (Q[first], p[last-1]) last--;if (Last-first <= 1) Return 0;getintersection (Q[LAST].P, Q[LAST].V, Q[FIRST].P, Q[FIRST].P, p[last]); int m = 0;for (int i = first; I &L T;= last; i++) poly[m++] = P[i];return m;}; Namespace Circular {using namespace Linear;using namespace vectorial;using namespace triangular;/* straight and original intersection */int Getlinec  Ircleintersection (point P, point Q, Circle O, double& t1, double& T2, vector<point>& sol) {vector v = q -p;/* to be emptied before use Sol *///sol.clear ();d ouble a = v.x, B = p.x-o.o.x, c = v.y, d = p.y-o.o.y;double e = a*a+c*c, F = (a*b)  +c*d), g = b*b+d*d-o.r*o.r;double Delta = f*f-4*e*g;if (dcmp (Delta) < 0) return 0;if (dcmp (delta) = = 0) {T1 = t2 =-F /(2 * e); Sol.push_back (P + v * t1); return 1;} T1 = (-f-sqrt (delta))/(2 * e); Sol.push_back (p + v * t1); t2 = (-f + sqrt (delta))/(2 * e); Sol.push_back (p + v * T2); return 2;} /* Circle intersection */int getcirclecircleintersection (Circle O1, Circle O2, vector<point>& sol) {Double d = getlength (O1. O-O2.O); if (dcmp (d) = = 0) {if (DCMP (o1.r-o2.R) = = 0) Return-1;return 0;} if (dcmp (O1.R + o2.r-d) < 0) return 0;if (DCMP (fabs)-D) > 0) return O1.R-O2.R a = 0;double (Getangle) ;d ouble da = ACOs ((O1.R*O1.R + D*D-O2.R*O2.R)/(2*o1.r*d)); Point P1 = O1.point (a-da), p2 = O1.point (A+da), Sol.push_back (p1), if (P1 = = p2) return 1;sol.push_back (P2); return 2;} /* tangents */int gettangents (point P, Circle O, vector* v) {Vector U = o.o-p;double d = getlength (u); if (d < O.R) R Eturn 0;else if (dcmp (D-O.R) = = 0) {V[0] = rotate (u, PI/2); return 1;} else {double ang = asin (O.R/D); v[0] = rotate (u ,-ang); v[1] = rotate (U, ANG); return 2;}} /* A[i] and b[i] are the tangency points of the section I tangent on O1 and O2 respectively */int gettangents (Circle O1, Circle O2, point* A, point* b) {int cnt = 0;IF (O1.R < O2.R) {swap (O1, O2); swap (A, b);} Double D2 = GetLength (O1.O-O2.O); D2 = D2 * d2;double rdif = O1.R-O2.R, rsum = O1.R + o2.r;if (D2 < RDIF * Rdif) return 0;if (dcmp (d2) = = 0 && DCMP (O1.R-O2.R) = = 0) return-1;double base = Getangle (O2.O-O1.O)if (dcmp (D2-RDIF * rdif) = = 0) {a[cnt] = O1.point (base); b[cnt] = O2.point (base); Cnt++;return CNT;} Double ang = ACOs ((O1.R-O2.R)/sqrt (D2)); a[cnt] = O1.point (Base+ang); B[CNT] = O2.point (Base+ang); CNT++;A[CNT] = O1.point (Base-ang); B[CNT] = O2.point (Base-ang); Cnt++;if (dcmp (d2-rsum * rsum) = = 0) {a[cnt] = O1.point (base); b[cnt] = O2.point (base); cnt++;} else if (D2 > Rsum * rsum) {Double ang = ACOs ((O1.R + O2.R)/sqrt (D2)); a[cnt] = O1.point (Base+ang); b[cnt] = O2.point (Base+ang); cnt++;a[cnt ] = O1.point (Base-ang); B[CNT] = O2.point (Base-ang); cnt++;} return CNT;}  /* Three points to determine the tangent circle */circle circumscribedcircle (Point P1, Spot, point P3) {Double Bx = p2.x-p1.x, by = P2.Y-P1.Y;  Double Cx = p3.x-p1.x, Cy = p3.y-p1.y;  Double D = 2 * (Bx * cy-by * Cx);  Double CX = (Cy * (BX * bx + by *) – by * (CX * CX + cy * cy))/D + p1.x;  Double cy = (BX * (CX * CX + cy * CY)-CX * (BX * bx + by *))/D + p1.y;  Point P = Point (CX, CY); Return Circle (P, GetLength (P1- p));  }/* three to determine inscribed circle */circle inscribedcircle (Point P1, points P2, dot p3) {Double A = GetLength (P2-P3);  Double b = getlength (P3-P1);  Double c = getlength (P1-P2);  Point P = (P1 * A + P2 * b + p3 * c)/(A + B + c);  Return Circle (P, Getdistancetoline (P, p1, p2));  }/* Triangle vertex is center */double getpublicareatotriangle (Circle O, point A, point B) {if (dcmp ((A-O.O) * (b-o.o)) = = 0) return 0;int sig = 1;double da = getplength (o.o-a), db = Getplength (O.o-b), if (dcmp (DA-DB) > 0) {swap (DA, db), swap (A, b); sig =-1;} Double T1, t2;vector<point> sol;int n = getlinecircleintersection (A, B, O, T1, T2, Sol); if (dcmp (DA-O.R*O.R) <= 0 {if (dcmp (DB-O.R*O.R) <= 0) return Getdirarea (O.O, A, b) * Sig;int k = 0;if (Getplength (sol[0]-b) > Getplength (sol[ 1]-b)) K = 1;double ret = Getarea (o.o, A, sol[k]) + O.getarea (Getangle (SOL[K]-O.O, b-o.o));d ouble tmp = (a-o.o) * (B-O.O); re TURN ret * SIG * DCMP (TMP);} Double d = getdistancetosegment (o.o, A, b); if (dcmp (D-O.R) >= 0) {DOUBLE ret = O.getArea (Getangle (A-O.O, b-o.o));d ouble tmp = (a-o.o) * (B-O.O); return ret * sig * DCMP (TMP);} Double K1 = o.r/getlength (a-o.o), K2 = O.r/getlength (B-O.O); Point P = o.o + (a-o.o) * K1, q = o.o + (b-o.o) * k2;double Ret1 = O.getarea (Getangle (P-O.O, q-o.o));d ouble Ret2 = o.g Etarea (Getangle (SOL[0]-O.O, SOL[1]-O.O))-Getarea (O.O, sol[0], sol[1]);d ouble ret = (ret1-ret2), TMP = (A-O.O) * (B-O.O); return ret * SIG * DCMP (TMP);} Double Getpublicareatopolygon (Circle O, point* p, int n) {if (dcmp (O.R) = = 0) return 0;double area = 0;for (int i = 0; i < n; i++) {int u = (i + 1)% n;area + = Getpublicareatotriangle (O, P[i], p[u]);} return Fabs (area);}}; Using namespace vectorial;using namespace Linear;const int maxp = 1e5 + 10;struct PSLG {struct Edge {int u, v;double ang;e Dge (int u = 0, int v = 0, double ang = 0): U (U), V (v), Ang (ANG) {}};int n, M, cnt;double Xi[maxp], Yi[maxp];vector<edge > edges;vector<int> g[maxp];int VIS[MAXP * 2], LEFT[MAXP * 2], PREV[MAXP * 2];vector<polyGon> faces;double area[maxp];void init (int n) {this->n = n;for (int i = 0; i < n; i++) g[i].clear (); Edges.clear () ; Faces.clear ();} Double Getangle (int u, int v) {return atan2 (Yi[v]-yi[u], xi[v]-xi[u]);} Double Getarea (Polygon p) {int n = p.size ();d ouble area = 0;for (int i = 0; i < n-1; i++) area + = (p[i]-p[0]) * (P[i+1] -p[0]); return AREA/2;} void Addedge (int u, int v) {Edges.push_back (Edge (U, V, Getangle (u,v))), Edges.push_back (Edge (V, U, Getangle (v,u))); m = EDG Es.size (); G[u].push_back (m-2); G[v].push_back (m-1);} void Build () {for (int u = 0; u < n; u++) {int k = G[u].size (), for (int i = 0; i < K; i++) {for (int j = i + 1; J & Lt K J + +) {if (Edges[g[u][i]].ang > Edges[g[u][j]].ang) Swap (G[u][i], g[u][j]);}} for (int i = 0; i < K; i++) prev[g[u][(i+1)%k]] = g[u][i];} memset (Vis, 0, sizeof (VIS)), cnt = 0;for (int u = 0; u < n; u++) {for (int i = 0; i < g[u].size (); i++) {int e = G[u] [I];if (!vis[e]) {cnt++; Polygon Poly;while (True) {vis[e] = 1;left[e] = Cnt;int D =Edges[e].u;poly.push_back (Point (Xi[d], yi[d])); e = prev[e^1];if (E = = G[u][i]) break; Faces.push_back (poly);}}} for (int i = 0; i < cnt; i++) Area[i] = Getarea (Faces[i]);}} plane;const int maxn = 105;int N, C; Point P[MAXN], V[MAXN * maxn];int find (Point u) {return Lower_bound (V, v + C, u)-V; void Init () {for (int i = 0; i < N; i++) P[i].read (); P[n] = p[0]; C = N; Point tmp;vector<double> dist[maxn];for (int i = 0, i < n; i++) V[i] = p[i];for (int i = 0; i < N; i++) {for ( int j = i+1; J < N; J + +) {if (Haveintersection (P[i], p[i+1], p[j], p[j+1])) {getintersection (P[i], p[i+1]-p[i], p[j], p[j+1]-p[j], TMP); V[c++] = Tmp;dist[i].push_back (GetLength (Tmp-p[i]));d Ist[j].push_back (GetLength (Tmp-p[j]));}}} Sort (V, v + C); c = Unique (V, v + C)-V;plane.init (c), for (int i = 0; i < C; i++) {plane.xi[i] = v[i].x;plane.yi[i] = v[i].y;} for (int i = 0; i < N; i++) {Vector v = (p[i+1]-p[i])/GetLength (P[i+1]-p[i]);d ist[i].push_back (0);d Ist[i].push_back (GetLength (P[i+1]-p[i])); Sort (Dist[i].begin (), Dist[i].end ()), int sz = Dist[i].size (); for (int j = 1; j < Sz; J + +) {Point A = P[i] + V * dist[ I][J-1]; Point B = p[i] + V * dist[i][j];if (A = = b) Continue;plane.addedge (Find (a), find (b));}} Plane.build ();} Polygon Simplify (const polygon& poly) {Polygon ret;int n = poly.size (); for (int i = 0; i < n; i++) {Point A = Poly [i]; Point B = poly[(i+1)%n]; Point C = poly[(i+2)%n];if (dcmp ((b-a) * (c-b))! = 0) ret.push_back (b);} return ret;} int main () {while (scanf ("%d", &n) = = 1 && N) {init (); Polygon poly;for (int i = 0; i < plane.faces.size (); i++) {if (Plane.area[i] < 0) {poly = Plane.faces[i];reverse (pol Y.begin (), Poly.end ());p Oly = simplify (poly); int start = 0, M = poly.size ();p rintf ("%d\n", M), for (int i = 0; i < m; i++) {if (Poly[i] < Poly[start]) start = i;} for (int i = start; i < m; i++) printf ("%.4lf%.4lf\n", poly[i].x, Poly[i].y), for (int i = 0; i < start; i++) printf (" %.4lf%.4lf\n ", poly[i].x, POLY[I].Y);} return 0;}

Copyright NOTICE: This article for Bo Master original article, without Bo Master permission not reproduced.

UVA 1340-find the Border (roll wrap)