计算几何模板

#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<vector>
#include<algorithm>
#define hhh puts("hhh")
#define see(x) (cerr<<(#x)<<'='<<(x)<<endl)
using namespace std;
typedef long long ll;
const int maxn = 1e5+10;
const int mod = 1e9+7;
const double eps = 1e-9;

//求正负符号
int sgn(double d) {
    if(fabs(d)<eps) return 0;
    if(d>0) return 1;
    return -1;
}

//求x,y大小关系,y默认为0
int dcmp(double x, double y=0) {
    if(fabs(x-y)<eps) return 0;
    if(x>y) return 1;
    return -1;
}

struct Point{
    double x, y;
    Point(double x=0, double y=0):x(x), y(y) {}
};

typedef Point Vector;

Vector operator + (Vector A, Vector B) {return Vector(A.x+B.x,A.y+B.y);}

Vector operator - (Point A, Point B) {return Vector(A.x-B.x,A.y-B.y);}

Vector operator * (Vector A, double p) {return Vector(A.x*p,A.y*p);}

Vector operator / (Vector A, double p) {return Vector(A.x/p,A.y/p);}


/*
typedef Point Vector;
Vector operator + (Vector A, Vector B) {
    Vector ret;
    ret.x = A.x + B.x, ret.y = A.y + B.y;
    return ret;
}

Vector operator - (Point A, Point B) {
    Vector ret;
    ret.x = A.x - B.x, ret.y = A.y - B.y;
    return ret;
}

Vector operator * (Vector A, double p) {
    Vector ret;
    ret.x = A.x * p, ret.y = A.y * p;
    return ret;
}


Vector operator / (Vector A, double p) {
    Vector ret;
    ret.x = A.x / p, ret.y = A.y / p;
    return ret;
}
*/

bool operator < (const Point &a, const Point &b) {
    if(a.x==b.x) return a.y<b.y;
    return a.x<b.x;
}

bool operator == (const Point &a, const Point &b) {
    if(sgn(a.x-b.x)==0 && sgn(a.y-b.y)==0) return true;
    return false;
}

//向量点乘,内积
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;
}

//向量长度
double Length(Vector A) {
    return sqrt(Dot(A,A));
}

//二维向量夹角,rad值
double Angle(Vector A, Vector B) {
    return acos(Dot(A,B)/Length(A)/Length(B));
}

//三角形面积的二倍
double Area2(Point A, Point B, Point C) {
    return Cross(B-A,C-A);
}

//向量逆时针旋转rad角度
Vector Rotate(Vector A, double rad) {
    return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad));
}

//向量A逆时针转90°后的单位向量
Vector Normal(Vector A) {
    double L=Length(A);
    return Vector(-A.y/L,A.x/L);
}

//判断bc是否在ab的逆时针方向,构造凸包时常用
bool ToLeftTest(Point a, Point b, Point c) {
    return Cross(b-a,c-b)>0;
}

struct Line{
    Point v, p;
    Line(){}
    Line(Point v, Point p):v(v), p(p) {}
    Point point(double t) {
        return v+(p-v)*t;
    }
};

//求两直线交点,但需保证Cross(v,w)!=0,即两直线不能夹直角
Point GetLineIntersection(Point P, Vector v, Point Q, Vector w) {
    Vector u=P-Q;
    double t=Cross(w,u)/Cross(v,w);
    return P+v*t;
}

//点P到直线AB的距离公式
double DistanceToLine(Point P, Point A, Point B) {
    Vector v1=B-A, v2=P-A;
    return fabs(Cross(v1,v2)/Length(v1));
}

//点P到线段AB的距离公式
double DistanceToSegment(Point P, Point A, Point B) {
    if(A==B) return Length(P-A);
    Vector v1=B-A, v2=P-A, v3=P-B;
    if(dcmp(Dot(v1,v2))<0) return Length(v2);
    if(dcmp(Dot(v1,v3))>0) return Length(v3);
    return DistanceToLine(P,A,B);
}

//点P在直线AB上的投影点
Point GetLineProjection(Point P, Point A, Point B) {
    Vector v=B-A;
    return A+v*(Dot(v,P-A)/Dot(v,v));
}

//判断点P是否在线段AB上,包含端点
bool OnSegment(Point p, Point a1, Point a2) {
    return dcmp(Cross(a1-p,a2-p))==0 && dcmp(Dot(a1-p,a2-p))<=0;
}

//判断两条线段是否相交
bool SegmentProperIntersection(Point a1, Point a2, Point b1, Point b2) {
    double c1=Cross(a2-a1,b1-a1), c2=Cross(a2-a1,b2-a1);
    double c3=Cross(b2-b1,a1-b1), c4=Cross(b2-b1,a2-b1);
    //若允许线段在端点处相交,则加入如下if判断,否则不需要
    if(!sgn(c1)||!sgn(c2)||!sgn(c3)||!sgn(c4)) {
        bool f1=OnSegment(b1,a1,a2);
        bool f2=OnSegment(b2,a1,a2);
        bool f3=OnSegment(a1,b1,b2);
        bool f4=OnSegment(a2,b1,b2);
        bool f=(f1|f2|f3|f4);
        return f;
    }
    return sgn(c1)*sgn(c2)<0 && sgn(c3)*sgn(c4)<0;
}

//求多边形面积,p为顶点集合,n为顶点个数
//方法为从第一个顶点出发把凸多边形分成n-2个三角形(另外一个常用的方法就自己写吧)
//顶点应按照逆时针顺序排列!!!关键点
double PolygonArea(Point *p, int n) {
    double s=0;
    for(int i=1; i<n-1; ++i) s+=Cross(p[i]-p[0],p[i+1]-p[0]);
    return s/2;
}

//判断点是否在多边形内(不保证凸多边形)
//若点在多边形内返回1,在多边形外部返回0,在多边形上返回-1
//若要判断在凸多边形内,则只需判断是否此点在所有边的左边
int isPointInPolygon(Point p, vector<Point> poly) {
    int wn=0;
    int n=poly.size();
    for(int i=0; i<n; ++i) {
        if(OnSegment(p,poly[i],poly[(i+1)%n])) return -1;
        int k=sgn(Cross(poly[(i+1)%n]-poly[i],p-poly[i]));
        int d1=sgn(poly[i].y-p.y);
        int d2=sgn(poly[(i+1)%n].y-p.y);
        if(k>0&&d1<=0&&d2>0) wn++;
        if(k<0&&d2<=0&&d1>0) wn--;
    }
    if(wn!=0) return 1;
    return 0;
}

struct Circle{
    Point c;
    double r;
    Circle(Point c, double r):c(c), r(r) {}
    Point point(double a) { //通过圆心角求坐标
        return Point(c.x+cos(a)*r,c.y+sin(a)*r);
    }
};

//求圆与直线的交点,交点存放在sol里面
int getLineCircleIntersection(Line L, Circle C, double &t1, double &t2, vector<Point> &sol) {
    double a=L.v.x, b=L.p.x-C.c.x, c=L.v.y, d=L.p.y-C.c.y;
    double e=a*a+c*c, f=2*(a*b+c*d), g=b*b+d*d-C.r*C.r;
    double delta=f*f-4*e*g;
    if(sgn(delta)<0) return 0; //相离
    if(sgn(delta)==0) { //相切
        t1=-f/(2*e);
        t2=-f/(2*e);
        sol.push_back(L.point(t1));
        return 1;
    }
    //相交
    t1=(-f-sqrt(delta))/(2*e);
    sol.push_back(L.point(t1));
    t2=(-f+sqrt(delta))/(2*e);
    sol.push_back(L.point(t2));
    return 2;
}

//两个圆相交的面积
double AreaOfOverlap(Point c1, double r1, Point c2, double r2) {
    double d=Length(c1-c2);
    if(r1+r2<d+eps) return 0.0;
    if(d<fabs(r1-r2)+eps) {
        double r=min(r1,r2);
        return acos(-1)*r*r;
    }
    double x=(d*d+r1*r1-r2*r2)/(2.0*d);
    double p=(r1+r2+d)/2.0;
    double t1=acos(x/r1);
    double t2=acos((d-x)/r2);
    double s1=r1*r1*t1;
    double s2=r2*r2*t2;
    double s3=2*sqrt(p*(p-r1)*(p-r2)*(p-d));
    return s1+s2-s3;
}

//求解凸包,GrahamScan算法,复杂度nlogn
//顶点编号从0到n-1
//凸包顶点(编号为排序后的编号)保存在stk数组中,从0到top-1
Point lst[maxn];
int stk[maxn], top;

bool cmp(Point p1, Point p2) {
    int tmp=sgn(Cross(p1-lst[0],p2-lst[0]));
    if(tmp>0) return true;
    if(tmp==0&&Length(lst[0]-p1)<Length(lst[0]-p2)) return true;
    return false;
}

void Graham(int n) {
    int k=0;
    Point p0;
    p0.x=lst[0].x;
    p0.y=lst[0].y;
    for(int i=1; i<n; ++i) {
        if(p0.y>lst[i].y || (p0.y==lst[i].y&&p0.x>lst[i].x)) {
            p0.x=lst[i].x;
            p0.y=lst[i].y;
            k=i;
        }
    }
    lst[k]=lst[0];
    lst[0]=p0;
    sort(lst+1,lst+n,cmp);
    if(n==1) { top=1; stk[0]=0; return; }
    if(n==2) { top=2; stk[0]=0; stk[1]=1; return; }
    stk[0]=0, stk[1]=1, top=2;
    for(int i=2; i<n; ++i) {
        while(top>1&&Cross(lst[stk[top-1]]-lst[stk[top-2]],lst[i]-lst[stk[top-2]])<=0) --top;
        stk[top]=i, ++top;
    }
}

//三角形
double sqr(double x){
	return x*x;
}
double dis(Point a,Point b){ //求ab的长度 
	return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y));
}
//三角形内心,角平分线的交点,到三角形三边的距离相同
Point Incenter(Point a,Point b,Point c){
     double A=dis(b,c);
     double B=dis(a,c);
     double C=dis(a,b);
     double S=A+B+C;
     double x=(A*a.x+B*b.x+C*c.x)/S;
     double y=(A*a.y+B*b.y+C*c.y)/S;
	 return Point(x,y);	
} 
//三条中线的交点,到三角形三顶点距离的平方和最小的点,三角形内到三边距离之积最大的点,重心和3个顶点组成的3个三角形面积相等。
Point gravity(Point a,Point b,Point c){ 
	double x=(a.x+b.x+c.x)/3;
	double y=(a.y+b.y+c.y)/3;
	return Point(x,y);
}
//三角形外心,三边中垂线交点,到三角形三个顶点距离相同
Point Circum(Point a,Point b,Point c){ 
	double x1=a.x,y1=a.y;
	double x2=b.x,y2=b.y;
	double x3=c.x,y3=c.y;
	
	double a1=2*(x2-x1);
	double b1=2*(y2-y1);
	double c1=x2*x2+y2*y2-x1*x1-y1*y1;
	
	double a2=2*(x3-x2);
	double b2=2*(y3-y2);
	double c2=x3*x3+y3*y3-x2*x2-y2*y2;
	
	double x=(c1*b2-c2*b1)/(a1*b2-a2*b1);
	double y=(a1*c2-a2*c1)/(a1*b2-a2*b1);
	
	return Point(x,y);
} 
//三角形垂心,三条高线的交点 
Point ortho(Point a,Point b,Point c){
	double A1=b.x-c.x;
	double B1=b.y-c.y;
	double C1=A1*a.y-B1*a.x;
	
	double A2=a.x-c.x;
	double B2=a.y-c.y;
	double C2=A2*b.y-B2*b.x;
	
	double x=(A1*C2-A2*C1)/(A2*B1-A1*B2);
	double y=(B1*C2-B2*C1)/(A2*B1-A1*B2);
	
	return Point(x,y);
}
//两条线段的叉积 
double crossProduct(double x1, double y1, double x2, double y2) {
    return x1 * y2 - x2 * y1;
}
//判断是否点是否在三角形内部 
bool isInside(double x, double y, double x1, double y1, double x2, double y2, double x3, double y3) {
    //特判是否形成三角形
    if (crossProduct(x3 - x1, y3 - y1, x2 - x1, y2 - y1) == 0) {
       return false;
    }
    //注意输入点的顺序不一定是逆时针,需要判断一下
    if(crossProduct(x3 - x1, y3 - y1, x2 - x1, y2 - y1) >= 0)
    {
        double tmpx = x2;
        double tmpy = y2;
        x2 = x3;
        y2 = y3;
        x3 = tmpx;
        y3 = tmpy;
    }
    if(crossProduct(x2 - x1, y2 - y1, x - x1, y - y1) < 0) return false;
    if(crossProduct(x3 - x2, y3 - y2, x - x2, y - y2) < 0) return false;
	if(crossProduct(x1 - x3, y1 - y3, x - x3, y - y3) < 0) return false;
    return true;
}

//最近点对模板
#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
const int MAXN=1e5+5;
struct Point{
	double x,y;	
}p[MAXN];
int temp[MAXN];
bool cmpx(const Point &a,const Point &b){ 
	return a.x<b.x;
}
bool cmpy(const int &a,const int &b){
	return p[a].y<p[b].y;
}
double dis(int a,int b){//两点之间的距离 
	return sqrt((p[a].x-p[b].x)*(p[a].x-p[b].x)+(p[a].y-p[b].y)*(p[a].y-p[b].y));
}
double min(double a,double b){//重载min函数 
	return a<b?a:b;
}
double closestPair(int l,int r)
{
	if(r==l)	return (1<<30);//只有一个点返回无穷大 
	if(r-l==1)	return dis(l,r); 
	int mid=(r+l)>>1;
	double d1=closestPair(l,mid);//递归左边 
	double d2=closestPair(mid+1,r);
	double d=min(d1,d2);
	int cnt=0;
	for(int i=l;i<=r;i++)	
		if(fabs(p[i].x-p[mid].x)<d)	temp[cnt++]=i;
	sort(temp,temp+cnt,cmpy);
	for(int i=0;i<cnt;i++)
		for(int j=i+1;j<=i+6 && j<cnt;j++){
			if(p[temp[j]].y-p[temp[i]].y>=d)	break;
			d=min(d,dis(temp[i],temp[j]));
		}
	return d;
}
int main(){
	int n;
	while(scanf("%d",&n),n){
		for(int i=0;i<n;i++)	scanf("%lf%lf",&p[i].x,&p[i].y);
		sort(p,p+n,cmpx);//对点进行x轴排序 
		printf("%.2f\n",closestPair(0,n-1)/2);
	}
	return 0;
}
//求圆心为(x0, y0), AB两个点的弧长公式
double arc(double x0, double y0, double x1, double y1, double x2, double y2) {
    double r,ab,k,b,h,so,C,l,O,s;
    r=sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0));
    double r2=(x1-x0)*(x1-x0)+(y1-y0)*(y1-y0);
    ab=sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
    double ab2=(x1-x2)*(x1-x2)+(y1-y2)*(y1-y2);
    double co=(2*r2-ab2)/(2*r2);
    O=acos(co);
    if(fabs(ab-2*r)<eps)
    {
        l=PI*r;
        return l;
    }
    else 			
    {
        if(O>PI)
        O=2*PI-O;
        l=O*r;
        return l;
    }		
}

详细介绍请看:(20条消息) 计算几何基础【用图来助你理解几何算法】_星空皓月的博客-CSDN博客

 (23条消息) 计算几何基础【用图来助你理解几何算法】_星空皓月的博客-CSDN博客

posted @   Landnig_on_Mars  阅读(7)  评论(0编辑  收藏  举报  
相关博文:
阅读排行:
· 地球OL攻略 —— 某应届生求职总结
· 周边上新:园子的第一款马克杯温暖上架
· Open-Sora 2.0 重磅开源!
· 提示词工程——AI应用必不可少的技术
· .NET周刊【3月第1期 2025-03-02】
点击右上角即可分享
微信分享提示