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
=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
=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保存逆时针遍历的点，也要去重共线的点
*/