「SOL」Bridge Across Islands(POJ)

对踵点是旋转卡壳的核心

---

# 题面

> Link POJ 3608

给两个凸包（凸多边形）(A,B)，保证不相交。求两凸包的最小距离。

数据规模：凸包的大小不超过 (10^4)。

---

# 解析

求一个凸包的直径可以用旋转卡壳，那两个凸包呢？

我们尝试将旋转的两条有向直线分别放在两个凸包上，并重新定义对踵点：

「定义」

对于凸包 $A$ 上的有向直线 $l=A_iA_{i+1}$，定义其方向为凸包的逆时针方向。

则在 $l$ 左侧的点到 $l$ 的有向距离为正，右侧有向距离为负，有向距离的绝对值为该点到 $l$ 的欧几里得距离。

$A_iA_{i+1}$ 的对踵点定义为 $B$ 上到 $l$ 的有向距离最大的点。

举个例子，下图点 (a) 是直线 (u) 的对踵点，点 (b) 是直线 (v) 的对踵点：

由于凸包之间的距离一定可以表示为其中一个凸包的点到另一个凸包的线段的距离（不一定是 (A) 上的点，(B) 上的线段；所以两种情况都要计算）。于是这样写旋转卡壳就可以求得凸包之间的距离。

---

# 源代码

/*Lucky_Glass*/
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;

namespace GEO{
	#define cdou const double &
	const double EPS=1e-16,INF=1e9;
	inline int sgn(cdou key){return fabs(key)<EPS? 0:(key<0? -1:1);}
	struct Vector{
		double x,y;
		Vector(cdou _x=0,cdou _y=0):x(_x),y(_y){}
		inline Vector operator *(cdou k)const{return Vector(x*k,y*k);}
		inline double len()const{return sqrt(x*x+y*y);}
		inline friend double cross(const Vector &u,const Vector &v){
			return u.x*v.y-u.y*v.x;
		}
		inline friend double dot(const Vector &u,const Vector &v){
			return u.x*v.x+u.y*v.y;
		}
		inline Vector operator -()const{return Vector(-x,-y);}
	};
	struct Point{
		double x,y;
		Point(cdou _x=0,cdou _y=0):x(_x),y(_y){}
		friend double distPoint(const Point &u,const Point &v){
			return sqrt((u.x-v.x)*(u.x-v.x)+(u.y-v.y)*(u.y-v.y));
		}
		inline Point operator +(const Vector &v)const{return Point(x+v.x,y+v.y);}
		inline Point operator -(const Vector &v)const{return Point(x-v.x,y-v.y);}
		inline Vector operator -(const Point &v)const{return Vector(x-v.x,y-v.y);}
		inline bool operator ==(const Point &v)const{return !sgn(x-v.x) && !sgn(y-v.y);}
		inline bool operator !=(const Point &v)const{return sgn(x-v.x) || sgn(y-v.y);}
	};
	struct Line{
		Point p;Vector d;
		double ang;
		Line(){}
		Line(const Point &_p,const Vector &_d):p(_p),d(_d){}
		Line(const Point &s,const Point &t):p(s),d(t-s){}
		void getAngle(){ang=atan2(d.y,d.x);}
	};
	//notice that if the head/tail of a segment is on the other segment, it will return 'true'
	bool ifSegInter(const Point &u1,const Point &u2,const Point &v1,const Point &v2){
		int re1=sgn(cross(v1-u1,u2-u1)*cross(v2-u1,u2-u1)),
			re2=sgn(cross(u1-v1,v2-v1)*cross(u2-v1,v2-v1));
		return re1<=0 && re2<=0;
	}
	//1=left / 0=on / -1=right
	int fixSide(const Point &s,const Point &t,const Point &p){return sgn(cross(t-s,p-s));}
	int fixSide(const Line &l,const Point &p){return sgn(cross(l.d,p-l.p));}
	Point getLineInter(const Line &u,const Line &v){
		double h=cross(u.p-v.p,v.d),d=cross(u.d,v.d);
		if(!sgn(d)) return u.p;
		return u.p-u.d*(h/d);
	}
	bool cmpPointToX(const Point &u,const Point &v){return sgn(u.x-v.x)? sgn(u.x-v.x)<0:sgn(u.y-v.y)<0;}
	bool cmpLineToAngle(const Line &u,const Line &v){return sgn(u.ang-v.ang)<0;}
	/*
	* notice:
	* 'org' is 0-indexed;
	* points on the edge of the convex is not counted
	*/
	bool buildConvex(Point *org,int n,Point *res,int &m){
		if(n<=2) return false;
		m=0;
		sort(org,org+n,cmpPointToX);
		for(int i=0;i<n;i++){
			while(m>1 && fixSide(res[m-2],res[m-1],org[i])<=0) m--;
			res[m++]=org[i];
		}
		int m0=m;
		for(int i=n-2;~i;i--){
			while(m>m0 && fixSide(res[m-2],res[m-1],org[i])<=0) m--;
			res[m++]=org[i];
		}
		m--;
		return m>2;
	}
	//notice: 'lin','rep' and 'rel' are 0-indexed
	bool getSI(Line *lin,int n,Point *rep,Line *rel,int &m){
		for(int i=0;i<n;i++) lin[i].getAngle();
		sort(lin,lin+n,cmpLineToAngle);
		int indl=0,indr=0;
		rel[indr++]=lin[0];
		for(int i=1;i<n;i++){
			while(indr-indl>1 && fixSide(lin[i],rep[indr-2])<=0) indr--;
			while(indr-indl>1 && fixSide(lin[i],rep[indl])<=0) indl++;
			if(!sgn(cross(lin[i].d,rel[indr-1].d))){
				if(fixSide(rel[indr-1],lin[i].p)>0){
					rel[indr-1]=lin[i];
					rep[indr-2]=getLineInter(rel[indr-1],rel[indr-2]);
				}
				continue;
			}
			rel[indr++]=lin[i];
			rep[indr-2]=getLineInter(rel[indr-1],rel[indr-2]);
		}
		while(indr-indl>1 && fixSide(rel[indl],rep[indr-2])<=0) indr--;
		rep[indr-1]=getLineInter(rel[indr-1],rel[indl]);
		m=indr-indl;
		if(m<=2) return false;
		for(int i=indl;i<indr;i++) rel[i-indl]=rel[i],rep[i-indl]=rep[i];
		return true;
	}
	double distLinePoint(const Line &l,const Point &p){
		return fabs(cross(l.d,p-l.p))/l.d.len();
	}
	double distSegPoint(const Point &s1,const Point &s2,const Point &p){
		if(sgn(dot(p-s1,s2-s1))>=0 && sgn(dot(p-s2,s1-s2))>=0) return distLinePoint(Line(s1,s2),p);
		return min(distPoint(p,s1),distPoint(p,s2));
	}
}
using namespace GEO;

const int N=1e4+10;

int n,m;
Point p1[N],p2[N];

#define debug(it) it.x,it.y

double directDist(const Point &s,const Point &t,const Point &p){
	return cross(t-s,p-s)/distPoint(s,t);
}
int main(){
	// freopen("input.in","r",stdin);
	while(~scanf("%d%d",&n,&m) && n && m){
		for(int i=0;i<n;i++) scanf("%lf%lf",&p1[i].x,&p1[i].y);
		for(int i=0;i<m;i++) scanf("%lf%lf",&p2[i].x,&p2[i].y);
		if(fixSide(p1[0],p1[1],p1[2])<0) reverse(p1,p1+n);
		if(fixSide(p2[0],p2[1],p2[2])<0) reverse(p2,p2+m);
		p1[n]=p1[0],p2[m]=p2[0];
		double ans=1e9;
		for(int i=0,j=0;i<n;i++){
			while(sgn(directDist(p1[i],p1[i+1],p2[j])-directDist(p1[i],p1[i+1],p2[j+1]))<0) j=(j+1)%m;
			ans=min(ans,distSegPoint(p1[i],p1[i+1],p2[j]));
			ans=min(ans,distSegPoint(p1[i],p1[i+1],p2[j+1]));
		}
		for(int i=0,j=0;i<m;i++){
			while(sgn(directDist(p2[i],p2[i+1],p1[j])-directDist(p2[i],p2[i+1],p1[j+1]))<0) j=(j+1)%n;
			ans=min(ans,distSegPoint(p2[i],p2[i+1],p1[j]));
			ans=min(ans,distSegPoint(p2[i],p2[i+1],p1[j+1]));
		}
		printf("%.6f
",ans);
	}
	return 0;
}

---

THE END

Thanks for reading!

[egin{split} “ &你是春山 你是岁酒\ &你是自由 你是误谬\ &你是颠沛流离之后 我绮丽的愁 ”\ ——& ext{《你是我遥不可及的梦》 By 苍穹} end{split} ]

> Link 你是我遥不可及的梦-Bilibili