Find the Border
PSGL
1 #include <bits/stdc++.h> 2 using namespace std; 3 #define pb push_back 4 const double eps = 1e-8; 5 const int inf = 0x3f3f3f3f; 6 7 int dcmp(double x) { 8 if(fabs(x) < eps) return 0; 9 return x < 0 ? -1 : 1; 10 } 11 12 struct Point { 13 double x,y; 14 Point (double x=0, double y=0): x(x), y(y) {} 15 }; 16 typedef Point Vector; 17 18 Vector operator + (Vector a, Vector b) { 19 return Vector(a.x+b.x, a.y+b.y); 20 } 21 Vector operator - (Point a, Point b) { 22 return Vector(a.x-b.x, a.y-b.y); 23 } 24 Vector operator * (Vector a, double p) { 25 return Vector(a.x*p, a.y*p); 26 } 27 Vector operator / (Vector a, double p) { 28 return Vector(a.x/p, a.y/p); 29 } 30 // 理论上这个“小于”运算符是错的，因为可能有三个点a, b, c, a和b很接近（即a<b好b<a都不成立），b和c很接近，但a和c不接近 31 // 所以使用这种“小于”运算符的前提是能排除上述情况 32 bool operator < (Point a, Point b) { 33 return dcmp(a.x-b.x)<0 || (dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)<0); 34 } 35 bool operator == (Point a, Point b) { 36 return dcmp(a.x-b.x)==0 && dcmp(a.y-b.y)==0; 37 } 38 double Dot(Vector a, Vector b) { 39 return a.x*b.x + a.y*b.y; 40 } 41 double Cross(Vector a, Vector b) { 42 return a.x*b.y - a.y*b.x; 43 } 44 double Length(Vector a) { 45 return sqrt(Dot(a,a)); 46 } 47 48 typedef vector<Point> Poly; 49 50 bool SegProIn(Point a1, Point b1, Point a2, Point b2) { 51 double c1 = Cross(b1-a1, b2-a1), c2 = Cross(b1-a1,a2-a1); 52 double c3 = Cross(b2-a2, a1-a2), c4 = Cross(b2-a2, b1-a2); 53 return dcmp(c1)*dcmp(c2)<0 && dcmp(c3)*dcmp(c4)<0; 54 } 55 Point GetLineIn(Point p, Vector v, Point q, Vector w) { 56 Vector u = p-q; 57 double t = Cross(w,u) / Cross(v,w); 58 return p+v*t; 59 } 60 bool OnSeg(Point p, Point a, Point b) { 61 return dcmp(Cross(a-p, b-p))==0 && dcmp(Dot(a-p, b-p))<0; 62 } 63 // 多边形的有向面积 64 double GetArea(Poly p) { 65 double area = 0; 66 int n = p.size(); 67 for(int i = 1; i < n-1; i++) area += Cross(p[i]-p[0], p[(i+1)%n]-p[0]); 68 return area/2; 69 } 70 71 struct Edge{ 72 int u, v; 73 double ang; 74 Edge(int u=0, int v=0, double ang=0): u(u), v(v), ang(ang){} 75 }; 76 const int maxe = 10000+10; // 最大边数 77 78 // 平面直线图（PSGL）实现 79 struct PSGL{ 80 int n, m, face_cnt; 81 double x[maxe], y[maxe]; 82 vector<Edge> edges; 83 vector<int> G[maxe]; 84 int vis[maxe<<1]; // 每条边是否已经访问过 85 int left[maxe<<1]; // 左面的编号 86 int pre[maxe<<1]; // 相同起点的上一条边（即顺时针旋转碰到的下一条边）的编号 87 88 vector<Poly> faces; 89 double area[maxe]; 90 91 void init(int n) { 92 this->n = n; 93 for(int i = 0; i < n; i++) G[i].clear(); 94 edges.clear(); 95 faces.clear(); 96 } 97 double getAngle(int u,int v) { 98 return atan2(y[v]-y[u], x[v]-x[u]); 99 } 100 void add(int u, int v) { 101 edges.pb(Edge(u,v,getAngle(u,v))); 102 edges.pb(Edge(v,u,getAngle(v,u))); 103 m = edges.size(); 104 G[u].pb(m-2); 105 G[v].pb(m-1); 106 } 107 //找出faces并计算面积 108 void build(){ 109 for(int u = 0; u < n; u++) { 110 int d = G[u].size(); 111 for(int i = 0; i < d; i++) { 112 for(int j = i+1; j < d; j++) { 113 if(edges[G[u][i]].ang > edges[G[u][j]].ang) swap(G[u][i],G[u][j]); 114 } 115 } 116 for(int i = 0; i < d; i++) pre[G[u][(i+1)%d]] = G[u][i]; 117 } 118 memset(vis, 0, sizeof(vis)); 119 face_cnt = 0; 120 for(int u = 0; u < n; u++) { 121 for(int i = 0; i < G[u].size(); i++){ 122 int e = G[u][i]; 123 if(!vis[e]){ //逆时针找圈 124 face_cnt++; 125 Poly p; 126 while(1){ 127 vis[e] = 1; 128 left[e] = face_cnt; 129 int from = edges[e].u; 130 p.pb(Point(x[from], y[from])); 131 e = pre[e^1]; 132 if(e == G[u][i]) break; 133 } 134 faces.pb(p); 135 } 136 } 137 } 138 for(int i = 0; i < faces.size(); i++) { 139 area[i] = GetArea(faces[i]); 140 141 } 142 } 143 144 }; 145 146 PSGL g; 147 148 const int maxv = 110; 149 int n, c; 150 151 Point p[maxv]; 152 Point V[maxv*maxv/2+maxv]; 153 154 int ID(Point p) { 155 return lower_bound(V, V+c, p) - V; 156 } 157 //假定poly没有相邻点重合的情况，只需要删除三点共线的情况 158 Poly simplify(const Poly& p) { 159 Poly temp; 160 int n = p.size(); 161 for(int i = 0; i < n; i++) { 162 Point a = p[i]; 163 Point b = p[(i+1)%n]; 164 Point c = p[(i+2)%n]; 165 if(dcmp(Cross(a-b, c-b)) != 0) temp.pb(b); 166 } 167 return temp; 168 } 169 170 void build_graph(){ 171 c = n; 172 for(int i = 0; i < n; i++) V[i] = p[i]; 173 vector<double> dis[maxv]; // dist[i][j]是第i条线段上的第j个点离起点（P[i]）的距离 174 for(int i = 0; i < n; i++) { 175 for(int j = i+1; j < n; j++) { 176 if(SegProIn(p[i], p[(i+1)%n], p[j], p[(j+1)%n])) { 177 Point xy = GetLineIn(p[i], p[(i+1)%n]-p[i], p[j], p[(j+1)%n]-p[j]); 178 V[c++] = xy; 179 dis[i].pb(Length(xy - p[i])); 180 dis[j].pb(Length(xy - p[j])); 181 } 182 } 183 } 184 // 为了保证“很接近的点”被看作同一个，这里使用了sort+unique的方法 185 // 必须使用前面提到的“理论上是错误”的小于运算符，否则不能保证“很接近的点”在排序后连续排列 186 // 另一个常见的处理方式是把坐标扩大很多倍（比如100000倍），然后四舍五入变成整点（计算完毕后再还原），用少许的精度损失换来鲁棒性和速度。 187 sort(V, V+c); 188 c = unique(V, V+c) - V; 189 190 g.init(c); //c是平面图的点数 191 for(int i = 0; i < c; i++) { 192 g.x[i] = V[i].x; 193 g.y[i] = V[i].y; 194 } 195 for(int i = 0; i < n; i++) { 196 Vector v = p[(i+1)%n] - p[i]; 197 double len = Length(v); 198 dis[i].pb(0); 199 dis[i].pb(len); 200 sort(dis[i].begin(), dis[i].end()); 201 int sz = dis[i].size(); 202 for(int j = 1; j < sz; j++) { 203 Point a = p[i] + v*(dis[i][j-1]/len); 204 Point b = p[i] + v*(dis[i][j]/len); 205 if(a == b) continue; 206 g.add(ID(a), ID(b)); 207 } 208 } 209 g.build(); 210 211 Poly po; 212 for(int i = 0; i < g.faces.size(); i++) { 213 if(g.area[i] < 0) { // 对于连通图，惟一一个面积小于零的面是无限面 214 po = g.faces[i]; 215 reverse(po.begin(),po.end()); // 对于内部区域来说，无限面多边形的各个顶点是顺时针的 216 po = simplify(po); // 无限面多边形上可能会有相邻共线点 217 break; 218 } 219 } 220 int m = po.size(); 221 printf("%d ", m); 222 int s = 0; 223 for(int i = 0; i < m; i++) if(po[i] < po[s]) s = i; 224 for(int i = s; i < m; i++) printf("%.4lf %.4lf ", po[i].x, po[i].y); 225 for(int i = 0; i < s; i++) printf("%.4lf %.4lf ", po[i].x, po[i].y); 226 } 227 228 229 int main(){ 230 while(scanf("%d", &n)!=EOF) { 231 int x,y; 232 for(int i = 0; i < n; i++) { 233 scanf("%d %d", &x, &y); 234 p[i] = Point(x,y); 235 } 236 build_graph(); 237 } 238 return 0; 239 }