LA 3218 Find the Border PSLG *

题目地址：https://vjudge.net/problem/UVALive-3218

LRJ的算法看半天才懂，实际上是先遍历所有区域，然后其中一个面积为负的就是外轮廓

其中边的遍历方向类似上图所示，只有外轮廓面积为负数

#include <bits/stdc++.h>
using namespace std;
#define REP(i,a,b)  for(int i=a;i<(int)(b);++i)
#define REPD(i,a,b) for(int i=a;i>(int)(b);--i)
const double PI=acos(-1);
const double EPS=1e-8;   
int dcmp(double x){
	if(fabs(x)<EPS) return 0;
	return x > 0 ? 1 : -1;
}
struct Point{
	double x,y;
	Point(double x=0,double y=0):x(x),y(y){}
};
typedef Point Vector;
bool operator < (const Point& p1, const Point& p2){ return dcmp(p1.x-p2.x)<0 || (dcmp(p1.x-p2.x)==0 && dcmp(p1.y-p2.y)<0) ;}
bool operator == (const Point& a, const Point& b) {	return dcmp(a.x-b.x)==0 && dcmp(a.y-b.y)==0; }
Vector operator / (const Point& A, double x){ return Vector{A.x/x, A.y/x};}
Vector operator * (const Vector& A, double x){ return Vector{A.x*x, A.y*x};}
Vector operator - (const Vector& A, const Vector& B){ return Vector{A.x-B.x,A.y-B.y};}
Vector operator + (const Point& A, const Vector& v){ return Point{A.x+v.x,A.y+v.y};}
Vector Rotate(const Point& p,double ang){ return Vector{p.x*cos(ang)-p.y*sin(ang),p.x*sin(ang)+p.y*cos(ang)};}
double Cross(const Vector& A, const Vector& B){ return A.x*B.y-A.y*B.x;}
double Dot(Vector A, Vector B){ return A.x*B.x+A.y*B.y;}
double Length(Vector v){ return sqrt(fabs(Dot(v,v)));}
Vector Unit(Vector v){ return v/Length(v);}
double Angle(Vector v) { return atan2(v.y, v.x); }


const int maxn=100+5;
Point P[maxn],V[maxn*maxn*2];
int n,Vnum;
struct Edge{
	int from,to; double ang;
};
vector<Edge> edges;
vector<vector<int> > G(maxn*maxn);
vector<vector<double> > dist(maxn*maxn);
typedef	 vector<Point> Polygon;
vector<Polygon> faces;
int Prev[maxn*maxn*2];
double Area[maxn];
bool vis[maxn*maxn*2];

bool ProperInter(const Point& a1, const Point& a2, const Point& b1, const Point& b2){
	double c1=Cross(a2-a1, b1-a1), c2=Cross(a2-a1, b2-a1),
	c3=Cross(b2-b1, a1-b1), c4=Cross(b2-b1, a2-b1);
	return dcmp(c1)*dcmp(c2)<0 && dcmp(c3)*dcmp(c4)<0;
}
Point GetInterP(const Point& P, const Vector& v, const Point& Q, const Vector& w){
	Vector u = P-Q;
    double t = Cross(w, u) / Cross(v, w);
    return P+v*t;
}
void AddEdge(int from, int to){
	edges.push_back(Edge{from, to, Angle(V[to]-V[from])});    //每条边都有反向边，这样每条边的左边可以确定区域
	edges.push_back(Edge{to, from, Angle(V[from]-V[to])});
	int m=edges.size();
	G[from].push_back(m-2);
	G[to].push_back(m-1);
}
int ID(Point p){
	return lower_bound(V, V+Vnum, p)-V;
}
Polygon simplifly(const Polygon& poly){
	Polygon ans;
	int m=poly.size();
	REP(i,0,m){
		Point a=poly[(i-1+m)%m];
		Point b=poly[i];
		Point c=poly[(i+1)%m];
		if(dcmp(Cross(a-b, c-b))!=0) ans.push_back(b);
	}
	return ans;
}
double PolygonArae(const Polygon& poly){
	double x=0;
	REP(i,1,poly.size()-1) x+=Cross(poly[i]-poly[0], poly[i+1]-poly[0]);
	return x/2;
}
void PSLG(){

	REP(u,0,Vnum){
		int m=G[u].size();
		REP(i,0,m) REP(j,i+1,m) 
		if(edges[G[u][i]].ang>edges[G[u][j]].ang) swap(G[u][i], G[u][j]);
		REP(i,0,m) Prev[G[u][(i+1)%m]]=G[u][i];     //i顺时针转的下一条是prev[i];
	}

	Polygon poly; faces.clear();
	memset(vis,false,sizeof(vis));
	REP(u,0,Vnum) REP(i,0,G[u].size()) {
		int e=G[u][i];
		if(!vis[e]){
			poly.clear();
			for(;;){
				vis[e]=true; poly.push_back(V[edges[e].from]);
				e=Prev[e^1];         //反向边的顺时针转遇到的第一条边，所以找到的外轮廓是顺时针的，面积为负
				if(e==G[u][i]) break;
				assert(vis[e]==false);  //不可能走一半发现边走过了，那就不是连通图了
			} 
			faces.push_back(poly);
		}
	}

	REP(i,0,faces.size()) Area[i]=PolygonArae(faces[i]);
}
void Solve(){
	memcpy(V,P,sizeof(Point)*n);  //V中保存去重的所有点

	Vnum=n;
	P[n]=P[0];
	REP(i,0,n) dist[i].clear();   // dist[i][j]是第i条线段上的第j个点离起点（P[i]）的距离
	REP(i,0,n) REP(j,i+1,n) 
	if(ProperInter(P[i], P[i+1], P[j], P[j+1])){
		Point pi=GetInterP(P[i], P[i+1]-P[i], P[j], P[j+1]-P[j]);  
		V[Vnum++]=pi;          //求出所有线段的交点
		dist[i].push_back(Length(pi-P[i]));
		dist[j].push_back(Length(pi-P[j]));
	}

	
	sort(V, V+Vnum);
	Vnum=unique(V, V+Vnum)-V;   //去重
	
	edges.clear();
	REP(i,0,Vnum) G[i].clear();
	REP(i,0,n) {
		Vector v=P[i+1]-P[i];
		double len=Length(v);
		dist[i].push_back(0);    //两个端点放入
		dist[i].push_back(len);
		sort(dist[i].begin(), dist[i].end());
		REP(j,0,dist[i].size()-1){
			Point a=P[i]+v*(dist[i][j]/len);
			Point b=P[i]+v*(dist[i][j+1]/len);
			if(a==b) continue;
			AddEdge(ID(a), ID(b));            //把所有线段包括有交点分割新产生线段存入
		}
	}

	PSLG();   //遍历
	Polygon poly;
	REP(i,0,faces.size())
	if(Area[i]<0) {      // 对于连通图，惟一一个面积小于零的面是无限面
		poly=faces[i];  
		reverse(poly.begin(), poly.end());  // 对于内部区域来说，无限面多边形的各个顶点是顺时针的
		break;
	}

	poly=simplifly(poly);   //除去共线的点

	int ans=poly.size();
	printf("%d\n", ans);

	int start=0;
	REP(i,0,ans) if(poly[i] < poly[start]) start=i;  //找出左下的点

	REP(i,start,ans) printf("%.4lf %.4lf\n", poly[i].x, poly[i].y);
	REP(i,0,start) printf("%.4lf %.4lf\n", poly[i].x, poly[i].y);
}
int main(int argc, char const *argv[])
{
	// freopen("input.in","r",stdin); 
	while(scanf("%d",&n)==1){
		REP(i,0,n) scanf("%lf%lf",&P[i].x,&P[i].y);
		Solve();
	}
	return 0;
}

下面是自己写的卷包裹法：

tips:求右拐最厉害的边可以借助反向边的逆时针第一条边

#include <bits/stdc++.h>
using namespace std;
#define REP(i,a,b)  for(int i=a;i<(int)(b);++i)
#define REPD(i,a,b) for(int i=a;i>(int)(b);--i)
const double PI=acos(-1);
const double EPS=1e-8;   
int dcmp(double x){
	if(fabs(x)<EPS) return 0;
	return x > 0 ? 1 : -1;
}
struct Point{
	double x,y;
	Point(double x=0,double y=0):x(x),y(y){}
};
typedef Point Vector;
bool operator < (const Point& p1, const Point& p2){ return dcmp(p1.x-p2.x)<0 || (dcmp(p1.x-p2.x)==0 && dcmp(p1.y-p2.y)<0) ;}
bool operator == (const Point& a, const Point& b) {	return dcmp(a.x-b.x)==0 && dcmp(a.y-b.y)==0; }
Vector operator / (const Point& A, double x){ return Vector{A.x/x, A.y/x};}
Vector operator * (const Vector& A, double x){ return Vector{A.x*x, A.y*x};}
Vector operator - (const Vector& A, const Vector& B){ return Vector{A.x-B.x,A.y-B.y};}
Vector operator + (const Point& A, const Vector& v){ return Point{A.x+v.x,A.y+v.y};}
Vector Rotate(const Point& p,double ang){ return Vector{p.x*cos(ang)-p.y*sin(ang),p.x*sin(ang)+p.y*cos(ang)};}
double Cross(const Vector& A, const Vector& B){ return A.x*B.y-A.y*B.x;}
double Dot(Vector A, Vector B){ return A.x*B.x+A.y*B.y;}
double Length(Vector v){ return sqrt(fabs(Dot(v,v)));}
Vector Unit(Vector v){ return v/Length(v);}
double Angle(Vector v) { return atan2(v.y, v.x); }

const int maxn=100+5,maxe=10000+5;
int n,vCnt;
Point P[maxn],V[maxe];
struct Edge{
	int from,to; double ang;
};
vector<Edge> edges;
bool cmp(int i, int j){
	return edges[i].ang < edges[j].ang;
}
typedef vector<Point> Polygon;
vector<vector<int> > Seg(maxe);
vector<vector<double> > dist(maxn);
int Next[maxe*2];
bool SegIntSeg(const Point& a, const Point& b, const Point& c, const Point& d){
	double c1=Cross(b-a, c-a), c2=Cross(b-a, d-a),
	c3=Cross(d-c, a-c), c4=Cross(d-c, b-c);
	return dcmp(c1)*dcmp(c2)<0 && dcmp(c3)*dcmp(c4)<0;
}
Point IntPoint(const Point& P, const Vector& v, const Point& Q, const Vector& w){
	Vector u=P-Q;
	double t=Cross(u, v)/Cross(w, v);
	return Q+w*t;
}
void AddEdge(int from, int to){
	edges.push_back(Edge{from, to, Angle(V[to]-V[from])});
	edges.push_back(Edge{to, from, Angle(V[from]-V[to])});
	int m=edges.size();
	Seg[from].push_back(m-2);
	Seg[to].push_back(m-1);
}
int ID(const Point& p){
	return lower_bound(V, V+vCnt, p) - V;
}
Polygon Simplify(const Polygon& poly){
	Polygon ans;
	int m=poly.size();
	REP(i,0,m) {
		Point a=poly[i], b=poly[(i+1)%m], c=poly[(i+2)%m];
		if(dcmp(Cross(b-a, c-a))!=0) ans.push_back(a); 
	}
	return ans;
}
void GiftWrapping(){
	int st=0;
	REP(i,1,vCnt) if(V[i] < V[st]) st=i;
	

	Polygon poly;
	int e=Seg[st][0];
	do{
		poly.push_back(V[edges[e].from]);
		e=Next[e^1];   //反向边逆时针第一条
	}while(edges[e].from!=st);

	poly=Simplify(poly);

	int m=poly.size();
	printf("%d\n", m);
	REP(i,0,m) printf("%.4lf %.4lf\n", poly[i].x, poly[i].y);
}
void Solve(){
	REP(i,0,n) { dist[i].clear(); V[i]=P[i]; }
	vCnt=n;
	P[n]=P[0];
	REP(i,0,n) REP(j,i+1,n) if(SegIntSeg(P[i], P[i+1], P[j], P[j+1])) {
		Point p=IntPoint(P[i], P[i+1]-P[i], P[j], P[j+1]-P[j]);
		V[vCnt++]=p;
		dist[i].push_back(Length(P[i]-p));
		dist[j].push_back(Length(P[j]-p));
	}

	sort(V, V+vCnt);
	vCnt=unique(V, V+vCnt)-V;

	REP(i,0,vCnt) Seg[i].clear(); edges.clear();
	REP(i,0,n) {
		Vector v=P[i+1]-P[i];
		double len=Length(v);
		dist[i].push_back(0);
		dist[i].push_back(len);
		sort(dist[i].begin(), dist[i].end());
		REP(j,1,dist[i].size()) {
			Point a=P[i]+v*(dist[i][j-1]/len);
			Point b=P[i]+v*(dist[i][j]/len);
			if(a==b) continue;
			AddEdge(ID(a), ID(b));
		}
	}

	REP(u,0,vCnt){
		sort(Seg[u].begin(), Seg[u].end(), cmp);
		int m=Seg[u].size();
		REP(i,0,m) Next[Seg[u][i]]=Seg[u][(i+1)%m];
 	}

	GiftWrapping();
}

int main(int argc, char const *argv[])
{
	// freopen("input.in","r",stdin);
	while(scanf("%d",&n)==1&&n) {
		REP(i,0,n) scanf("%lf%lf",&P[i].x, &P[i].y);
		Solve();
	}
	return 0;
}
// 1.求出所有交点，每个点连接的线段按极角序排
// 2.从左下开始遍历，逆时针开始遍历，每次遇到第一个点就往反向边逆时针转的的第一条走去，直到回到起点
/*
那么问题就是如何保存每个点的边：
1.首先交点肯定都要求出来保存在V[]中，去重
在求1的同时也要把所有线段求出来，并且线段的起点也要知道(用vector，Seg[a][i]保存a点相连的边，注意按极角排序)，边的反向边也要容易得到(e^1方法得到)
要个数组保存边e逆时针转的第一条边

所以：
vector<double> dist[]; dist[u]保存与u成直线的点 
edges保存所有边
G[u][i]保存u出发的在edges中的编号，且已经排序
prev[e]保存e逆时针转的第一条边
逆时针卷包裹的时候不会遇到已经遍历的边，所以不用vis判重
vector<Point> ans保存逆时针遍历的点，也要去重共线的点
*/