计算几何模板
2017-07-31 本文已影响0人
1QzUPm_09F
/*///////////////////////////////////
***************Content***************
范数（模的平方）
向量的模(求线段距离/两点距离)
点乘
叉乘
判断向量垂直
判断a1-a2与b1-b2垂直
判断线段垂直
判断向量平行
判断a1-a2与b1-b2平行
判断线段平行
求点在线段上的垂足
求点关于直线的对称点
判断P0,P1,P2三点位置关系,Vector p0-p2 在p0-p1的相对位置
判断p1-p2与p3-p4是否相交
判断线段是否相交
a,b两点间的距离
点到直线的距离
点到线段的距离
两线段间的距离(相交为0)
求两线段的交点
求直线交点
中垂线
向量a,b的夹角，范围[0,180]
向量（点）极角，范围[-180,180]
角度排序(从x正半轴起逆时针一圈)范围为[0,180)
判断点在多边形内部
凸包(CCW/CW)
求向量A,向量B构成三角形的面积
(旋转卡壳)求平面中任意两点的最大距离
三点所成的外接圆
三点所成的内切圆
过点p求过该点的两条切线与圆的两个切点
极角(直线与x轴的角度)
点与圆的位置关系
已知点与切线求圆心
已知两直线和半径,求夹在两直线间的圆
求与给定两圆相切的圆
多边形(存储方式:点->线)
半平面交
求最近点对的距离(分治)
*////////////////////////////////////
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
#include<vector>
using namespace std;
#define EPS (1e-10)
#define equals(a, b) (fabs(a-b)<EPS)
const double PI = acos(-1);
const double INF = 1e20;
static const int CCW=1;//逆时针
static const int CW=-1;//顺时针
static const int BACK=-2;//后面
static const int FRONT=2;//前面
static const int ON=0;//线段上
struct Point
{
double x, y;
Point(double x=0, double y=0):x(x), y(y) {}
};//点
typedef Point Vector;//向量
struct Segment
{
Point p1, p2;
Segment(Point p1=Point(), Point p2=Point()):p1(p1), p2(p2) {}
double angle;
};//线段
typedef Segment Line;//直线
typedef vector<Point> Polygon;//多边形(存点)
typedef vector<Segment> Pol;//多边形（存边）
class Circle
{
public:
Point c;
double r;
Circle(Point c=Point(), double r=0.0):c(c), r(r) {}
Point point(double a)
{
return Point(c.x + cos(a)*r, c.y + sin(a)*r);
}
};//圆
Point operator + (Point a, Point b)
{
return Point(a.x+b.x, a.y+b.y);
}
Point operator - (Point a, Point b)
{
return Point(a.x-b.x, a.y-b.y);
}
Point operator * (Point a, double p)
{
return Point(a.x*p, a.y*p);
}
Point operator / (Point a, double p)
{
return Point(a.x/p, a.y/p);
}
//排序左下到右上
bool operator < (const Point &a,const Point &b)
{
return a.x<b.x||(a.x==b.x&&a.y<b.y);
}
bool operator == (Point a, Point b)
{
return fabs(a.x-b.x)<EPS && fabs(a.y-b.y)<EPS;
}
//范数（模的平方）
double norm(Vector a)
{
return a.x*a.x+a.y*a.y;
}
//向量的模(求线段距离/两点距离)
double abs(Vector a)
{
return sqrt(norm(a));
}
//点乘
double dot(Vector a, Vector b)
{
return a.x*b.x+a.y*b.y;
}
//叉乘
double cross(Vector a, Vector b)
{
return a.x*b.y-a.y*b.x;
}
//判断向量垂直
bool isOrthgonal(Vector a, Vector b)
{
return equals(dot(a, b), 0.0);
}
//判断a1-a2与b1-b2垂直
bool isOrthgonal(Point a1, Point a2, Point b1, Point b2)
{
return isOrthgonal(a1-a2, b1-b2);
}
//判断线段垂直
bool isOrthgonal(Segment s1, Segment s2)
{
return equals(dot(s1.p2-s1.p1, s2.p2-s2.p1), 0.0);
}
//判断向量平行
bool isParallel(Vector a, Vector b)
{
return equals(cross(a, b), 0.0);
}
//判断a1-a2与b1-b2平行
bool isParallel(Point a1, Point a2, Point b1, Point b2)
{
return isParallel(a1-a2, b1-b2);
}
//判断线段平行
bool isParallel(Segment s1, Segment s2)
{
return equals(cross(s1.p2-s1.p1, s2.p2-s2.p1), 0.0);
}
//求点在线段上的垂足
Point project(Segment s, Point p)
{
Vector base=s.p2-s.p1;
double r=dot(p-s.p1, base)/norm(base);
return s.p1+base*r;
}
//求点关于直线的对称点
Point reflect(Segment s, Point p)
{
return p+(project(s, p)-p)*2.0;
}
//判断P0,P1,P2三点位置关系,Vector p0-p2 在p0-p1的相对位置
int ccw(Point p0, Point p1, Point p2)
{
Vector a=p1-p0;
Vector b=p2-p0;
if( cross(a, b)>EPS ) return CCW;
if( cross(a, b)<-EPS ) return CW;
if( dot(a, b)<-EPS ) return BACK;
if( norm(a)<norm(b) ) return FRONT;
return ON;
}
//判断p1-p2与p3-p4是否相交
bool intersect(Point p1, Point p2, Point p3, Point p4)
{
return (ccw(p1, p2, p3)*ccw(p1, p2, p4)<=0 &&ccw(p3, p4, p1)*ccw(p3, p4, p2)<=0);
}
//判断线段是否相交
bool intersect(Segment s1, Segment s2)
{
return intersect(s1.p1, s1.p2, s2.p1, s2.p2);
}
//a,b两点间的距离
double getDistance(Point a, Point b)
{
return abs(a-b);
}
//点到直线的距离
double getDistanceLP(Line l, Point p)
{
return abs(cross(l.p2-l.p1, p-l.p1)/abs(l.p2-l.p1));
}
//点到线段的距离
double getDistanceSP(Segment s, Point p)
{
if(dot(s.p2-s.p1, p-s.p1)<0.0) return abs(p-s.p1);
if(dot(s.p1-s.p2, p-s.p2)<0.0) return abs(p-s.p2);
return getDistanceLP(s, p);
}
//两线段间的距离(相交为0)
double getDistanceSS(Segment s1, Segment s2)
{
if(intersect(s1, s2)) return 0.0;
return min( min(getDistanceSP(s1,s2.p1), getDistanceSP(s1,s2.p2)),
min(getDistanceSP(s2,s1.p1), getDistanceSP(s2,s1.p2)) );
}
//求两线段的交点
Point getCrossPoint(Segment s1, Segment s2)
{
Vector base=s2.p2-s2.p1;
double d1=abs(cross(base, s1.p1-s2.p1));
double d2=abs(cross(base, s1.p2-s2.p1));
double t=d1/(d1+d2);
return s1.p1+(s1.p2-s1.p1)*t;
}
//求直线交点
Point intersectL(Segment a, Segment b)
{
double x1=a.p1.x,y1=a.p1.y,x2=a.p2.x,y2=a.p2.y;
double x3=b.p1.x,y3=b.p1.y,x4=b.p2.x,y4=b.p2.y;
double k1=(x4-x3)*(y2-y1),k2=(x2-x1)*(y4-y3);
double ans_x=(k1*x1-k2*x3+(y3-y1)*(x2-x1)*(x4-x3))/(k1-k2);
double ans_y=(k2*y1-k1*y3+(x3-x1)*(y2-y1)*(y4-y3))/(k2-k1);
return Point(ans_x,ans_y);
}
//中垂线
Line mid_vert(Line l)
{
double x1=l.p1.x,y1=l.p1.y;
double x2=l.p2.x,y2=l.p2.y;
double xm=(x1+x2)/2,ym=(y1+y2)/2;
Line s;
s.p1.x=xm+ym-y1;
s.p1.y=ym-xm+x1;
s.p2.x=xm-ym+y1;
s.p2.y=ym+xm-x1;
return s;
}
//向量a,b的夹角，范围[0,180]
double Angle(Vector a, Vector b)
{
return acos(dot(a, b)/(abs(a)*abs(b)));
}
//向量（点）极角，范围[-180,180]
double angle(Vector v)
{
return atan2(v.y, v.x);
}
//角度排序(从x正半轴起逆时针一圈)范围为[0,180)
double SortAngle(Vector a)
{
Point p0(0.0, 0.0);
Point p1(a.x, a.y);
Point p2(-1.0, 0.0);
Vector b=p2;
if(ccw(p0, p1, p2)==CW) return acos(dot(a, b)/(abs(a)*abs(b)));
if(ccw(p0, p1, p2)==CCW) return 2*PI-acos(dot(a, b)/(abs(a)*abs(b)));
if(ccw(p0, p1, p2)==BACK) return PI;
else return 0;
}
/*
判断点在多边形内部
IN 2
ON 1
OUT 0
*/
int contains(Polygon g, Point p)
{
int n=g.size();
bool x=false;
for(int i=0; i<n; i++)
{
Point a=g[i]-p;
Point b=g[(i+1)%n]-p;
if(abs(cross(a, b))<EPS && dot(a, b)<EPS) return 1;
if(a.y>b.y) swap(a,b);
if(a.y<EPS && EPS<b.y && cross(a,b)>EPS) x=!x;
}
return (x? 2 : 0);
}
//凸包(CCW/CW)
Polygon andrewScan(Polygon s)
{
Polygon u, l;
if(s.size()<3) return s;
sort(s.begin(), s.end());
u.push_back(s[0]);
u.push_back(s[1]);
l.push_back(s[s.size()-1]);
l.push_back(s[s.size()-2]);
for(int i=2; i<s.size(); i++)
{
for(int n=u.size(); n>=2 && ccw(u[n-2], u[n-1], s[i])!=CW; n--)
u.pop_back();
u.push_back(s[i]);
}
for(int i=s.size()-3; i>=0; i--)
{
for(int n=l.size(); n>=2 && ccw(l[n-2], l[n-1], s[i])!=CW; n--)
l.pop_back();
l.push_back(s[i]);
}
reverse(l.begin(), l.end());
for(int i=u.size()-2; i>=1; i--) l.push_back(u[i]);
return l;
}
//求向量A,向量B构成三角形的面积
double TriArea(Vector a, Vector b)
{
return 0.5*abs(cross(a,b));
}
//求平面中任意两点的最大距离(旋转卡壳)
double RotatingCalipers(const Polygon& s)
{
Polygon l;
double dis, maxn=0.0;
int len, i, k;
l=andrewScan(s);
len=l.size();
if(len>=3)
{
for(i=0, k=2; i<len; i++)
{
while(cross(l[(k+1)%len]-l[i], l[(k+1)%len]-l[(i+1)%len])>=cross(l[k%len]-l[i], l[k%len]-l[(i+1)%len]))
k++;
dis=max(norm(l[k%len]-l[i]), norm(l[k%len]-l[(i+1)%len]));
if(dis>maxn) maxn=dis;
}
}
else maxn=norm(l[1]-l[0]);
return maxn;
}
//三点所成的外接圆
Circle CircumscribedCircle(Point a, Point b, Point c)
{
double x=0.5*(norm(b)*c.y+norm(c)*a.y+norm(a)*b.y-norm(b)*a.y-norm(c)*b.y-norm(a)*c.y)
/(b.x*c.y+c.x*a.y+a.x*b.y-b.x*a.y-c.x*b.y-a.x*c.y);
double y=0.5*(norm(b)*a.x+norm(c)*b.x+norm(a)*c.x-norm(b)*c.x-norm(c)*a.x-norm(a)*b.x)
/(b.x*c.y+c.x*a.y+a.x*b.y-b.x*a.y-c.x*b.y-a.x*c.y);
Point O(x, y);
double r=abs(O-a);
Circle m(O, r);
return m;
}
//三点所成的内切圆
Circle InscribedCircle(Point a, Point b, Point c)
{
double A=abs(b-c), B=abs(a-c), C=abs(a-b);
double x=(A*a.x+B*b.x+C*c.x)/(A+B+C);
double y=(A*a.y+B*b.y+C*c.y)/(A+B+C);
Point O(x, y);
Line l(a, b);
double r=getDistanceLP(l, O);
Circle m(O, r);
return m;
}
//过点p求过该点的两条切线与圆的两个切点
Segment TangentLineThroughPoint(Circle m, Point p)
{
Point c=m.c;
double l=abs(c-p);
double r=m.r;
double k=(2*r*r-l*l+norm(p)-norm(c)-2*p.y*c.y+2*c.y*c.y)/(2*(p.y-c.y));
double A=1+(p.x-c.x)*(p.x-c.x)/((p.y-c.y)*(p.y-c.y));
double B=-(2*k*(p.x-c.x)/(p.y-c.y)+2*c.x);
double C=c.x*c.x+k*k-r*r;
double x1, x2, y1, y2;
x1=(-B-sqrt(B*B-4*A*C))/(2*A);
x2=(-B+sqrt(B*B-4*A*C))/(2*A);
y1=(2*r*r-l*l+norm(p)-norm(c)-2*(p.x-c.x)*x1)/(2*(p.y-c.y));
y2=(2*r*r-l*l+norm(p)-norm(c)-2*(p.x-c.x)*x2)/(2*(p.y-c.y));
Point p1(x1, y1), p2(x2, y2);
Segment L(p1, p2);
return L;
}
//极角(直线与x轴的角度)
double PolarAngle(Vector a)
{
Point p0(0.0, 0.0);
Point p1(1.0, 0.0);
Point p2(a.x, a.y);
Vector b=p1;
double ans=0;
if(ccw(p0, p1, p2)==CW) ans=180-acos(dot(a, b)/(abs(a)*abs(b)))*180/acos(-1);
else if(ccw(p0, p1, p2)==CCW) ans=acos(dot(a, b)/(abs(a)*abs(b)))*180/acos(-1);
else ans=0;
if(ans>=180) ans-=180;
if(ans<0) ans+=180;
return ans;
}
//点与圆的位置关系
int CircleContain(Circle m, Point p)
{
double r=m.r;
double l=abs(p-m.c);
if(r>l) return 2;
if(r==l) return 1;
if(r<l) return 0;
}
//已知点与切线求圆心
void CircleThroughAPointAndTangentToALineWithRadius(Point p, Line l, double r)
{
Point m=project(l, p);
if(abs(p-m)>2*r)
{
printf("[]\n");
}
else if(abs(p-m)==2*r)
{
Circle c((p+m)/2, r);
printf("[(%.6f,%.6f)]\n", c.c.x, c.c.y);
}
else if(abs(p-m)<EPS)
{
Point m0(m.x+10, m.y);
if(abs(m0-project(l, m0))<EPS) m0.y+=20;
Point m1=project(l, m0);
Circle c1(m-(m0-m1)/abs(m0-m1)*r, r);
Circle c2(m+(m0-m1)/abs(m0-m1)*r, r);
if(c1.c.x>c2.c.x) swap(c1, c2);
else if(c1.c.x==c2.c.x && c1.c.y>c2.c.y) swap(c1, c2);
printf("[(%.6f,%.6f),(%.6f,%.6f)]\n", c1.c.x, c1.c.y, c2.c.x, c2.c.y);
}
else if(abs(p-m)<2*r)
{
double s=abs(p-m);
double d=sqrt(r*r-(r-s)*(r-s));
Point m1, m2;
m1=(m+(l.p1-l.p2)/abs(l.p1-l.p2)*d);
m2=(m-(l.p1-l.p2)/abs(l.p1-l.p2)*d);
Circle c1(m1+(p-m)/abs(p-m)*r, r);
Circle c2(m2+(p-m)/abs(p-m)*r, r);
if(c1.c.x>c2.c.x) swap(c1, c2);
else if(c1.c.x==c2.c.x && c1.c.y>c2.c.y) swap(c1, c2);
printf("[(%.6f,%.6f),(%.6f,%.6f)]\n", c1.c.x, c1.c.y, c2.c.x, c2.c.y);
}
return ;
}
bool cmp_CircleTangentToTwoLinesWithRadius(Circle x, Circle y)
{
if(x.c.x==y.c.x) return x.c.y<y.c.y;
else return x.c.x<y.c.x;
}
//已知两直线和半径,求夹在两直线间的圆
void CircleTangentToTwoLinesWithRadius(Line l1, Line l2, double r)
{
Point p=intersectL(l1, l2);
Vector a, b;
l1.p2=p+p-l1.p1;
l2.p2=p+p-l2.p1;
a=(l1.p1-p)/abs(l1.p1-p), b=(l2.p2-p)/abs(l2.p2-p);
Circle c1(p+(a+b)/abs(a+b)*(r/sin(Angle(a,a+b))),r);
a=(l1.p1-p)/abs(l1.p1-p), b=(l2.p1-p)/abs(l2.p1-p);
Circle c2(p+(a+b)/abs(a+b)*r/sin(Angle(a,a+b)),r);
a=(l1.p2-p)/abs(l1.p2-p), b=(l2.p1-p)/abs(l2.p1-p);
Circle c3(p+(a+b)/abs(a+b)*(r/sin(Angle(a,a+b))),r);
a=(l1.p2-p)/abs(l1.p2-p), b=(l2.p2-p)/abs(l2.p2-p);
Circle c4(p+(a+b)/abs(a+b)*(r/sin(Angle(a,a+b))),r);
vector<Circle> T;
T.push_back(c1);
T.push_back(c2);
T.push_back(c3);
T.push_back(c4);
sort(T.begin(), T.end(), cmp_CircleTangentToTwoLinesWithRadius);
printf("[(%.6f,%.6f),(%.6f,%.6f),(%.6f,%.6f),(%.6f,%.6f)]\n", T[0].c.x, T[0].c.y, T[1].c.x, T[1].c.y, T[2].c.x, T[2].c.y, T[3].c.x, T[3].c.y);
}
//求与给定两圆相切的圆
void CircleTangentToTwoDisjointCirclesWithRadius(Circle C1, Circle C2)
{
double d = abs(C1.c-C2.c);
if(d<EPS) printf("[]\n");
else if(fabs(C1.r+C2.r)<d) printf("[]\n");
else if(fabs(C1.r-C2.r)>d) printf("[]\n");
else
{
double sita = angle(C2.c - C1.c);
double da = acos((C1.r*C1.r + d*d - C2.r*C2.r) / (2 * C1.r*d));
Point p1 = C1.point(sita - da), p2 = C1.point(sita + da);
if(p1.x>p2.x) swap(p1, p2);
else if(p1.x==p2.x && p1.y>p2.y) swap(p1, p2);
if (p1 == p2) printf("[(%.6f,%.6f)]\n", p1.x, p1.y);
else printf("[(%.6f,%.6f),(%.6f,%.6f)]\n", p1.x, p1.y, p2.x, p2.y);
}
}
//多边形(存储方式:点->线)
Pol Polygon_to_Pol(Polygon L)
{
Pol S;
Segment l;
for(int i=0; i+1<L.size(); i++)
{
l.p1=L[i];
l.p2=L[i+1];
swap(l.p1, l.p2);//注意顺序
l.angle=angle(l.p2-l.p1);
S.push_back(l);
}
l.p1=L[L.size()-1];
l.p2=L[0];
swap(l.p1, l.p2);//注意顺序
l.angle=angle(l.p2-l.p1);
S.push_back(l);
return S;
}
//半平面交
bool SortAnglecmp(Segment a, Segment b)
{
if(fabs(a.angle-b.angle)>EPS)
return a.angle>b.angle;
return ccw(b.p1, b.p2, a.p1)!=CW;
}
int intersection_of_half_planes(Pol s)
{
Segment deq[1505];
Segment l[1505];
Point p[1505];
memset(deq, 0, sizeof(deq));
memset(l, 0, sizeof(l));
memset(p, 0, sizeof(p));
sort(s.begin(), s.end(), SortAnglecmp);
int cnt=0;
for(int i=0; i<s.size(); i++)
if(fabs(s[i].angle-l[cnt].angle)>EPS)
l[++cnt]=s[i];
int le=1,ri=1;
for(int i=1; i<=cnt; i++)
{
while(ri>le+1 && ccw(l[i].p1, l[i].p2, intersectL(deq[ri-1],deq[ri-2]))==CW) ri--;
while(ri>le+1 && ccw(l[i].p1, l[i].p2, intersectL(deq[le],deq[le+1]))==CW) le++;
deq[ri++]=l[i];
}
while(ri>le+2 && ccw(deq[le].p1, deq[le].p2, intersectL(deq[ri-1],deq[ri-2]))==CW) ri--;
while(ri>le+2 && ccw(deq[ri-1].p1, deq[ri-1].p2, intersectL(deq[le],deq[le+1]))==CW) le++;
//***************getArea******************
/*
if(ri<=le+2)
{
printf("0.00\n");
return 0;
}
deq[ri]=deq[le];
cnt=0;
for(int i=le; i<ri; i++)
p[++cnt]=intersectL(deq[i],deq[i+1]);
double ans=0.0;
for(int i=2; i<cnt; i++)
ans+=fabs(TriArea(p[i]-p[1],p[i+1]-p[1]));
printf("%.2f\n", ans);
return 0;
*/
//****************************************
if(ri>le+2) return 1;
return 0;
}
//求最近点对的距离(分治)
//************************************************************************
bool cmpxy(Point a, Point b)
{
if(a.x != b.x) return a.x < b.x;
else return a.y < b.y;
}
bool cmpy(Point a, Point b)
{
return a.y < b.y;
}
int n;
Point closest_p[100010];//将点都存入closest_p中！
Point closest_tmpt[100010];
double Closest_Pair(int left,int right)
{
double d = INF;
if(left == right) return d;
if(left + 1 == right)
return getDistance(closest_p[left],closest_p[right]);
int mid = (left+right)/2;
double d1 = Closest_Pair(left,mid);
double d2 = Closest_Pair(mid+1,right);
d = min(d1,d2);
int k = 0;
for(int i = left; i <= right; i++)
if(fabs(closest_p[mid].x - closest_p[i].x) <= d)
closest_tmpt[k++] = closest_p[i];
sort(closest_tmpt,closest_tmpt+k,cmpy);
for(int i = 0; i <k; i++)
for(int j = i+1; j < k && closest_tmpt[j].y - closest_tmpt[i].y < d; j++)
d = min(d,getDistance(closest_tmpt[i],closest_tmpt[j]));
return d;
}
double Closest()//直接调用此函数即可
{
sort(closest_p,closest_p+n,cmpxy);
return Closest_Pair(0,n-1)/2;
}
//************************************************************************
int main()
{
while(scanf("%d",&n)==1 && n)
{
for(int i = 0; i < n; i++)
scanf("%lf%lf",&closest_p[i].x,&closest_p[i].y);
printf("%.2lf\n",Closest());
}
return 0;
}
