二维计算几何模板--圆

转载请注明出处 http://blog.csdn.net/xtulollipop

#include<cstdio>#include<cmath>#include<cstring>#include<iostream>#include<algorithm>#include<cstdlib>#include<queue>#include<vector>#include<map>#include<stack>#include<set>#include<complex>#define EPS 1e-6    //log(x)#define e exp(1.0); //2.718281828#define mod 1000000007#define INF 0x7fffffff#define inf 0x3f3f3f3ftypedef long long LL;using namespace std;//基本函数const double eps=1e-8;int cmp(double x) {    if(fabs(x)<eps) return 0;    if(x>0) return 1;    return -1;}const double pi=acos(-1.0);inline double sqr(double x) {    return x*x;}struct point {    double x,y;    point() {};    point(double a,double b):x(a),y(b) {};    void input() {        scanf("%lf %lf",&x,&y);    }    friend point operator + (const point &a,const point &b) {        return point(a.x+b.x,a.y+b.y);    }    friend point operator - (const point &a,const point &b) {        return point(a.x-b.x,a.y-b.y);    }    friend bool operator == (const point &a,const point &b) {        return cmp(a.x-b.x)==0&&cmp(a.y-b.y)==0;    }    friend point operator * (const point &a,const double &b) {        return point(a.x*b,a.y*b);    }    friend point operator * (const double &a,const point &b) {        return point(a*b.x,a*b.y);    }    friend point operator / (const point &a,const double &b) {        return point(a.x/b,a.y/b);    }    double norm() {        return sqrt(sqr(x)+sqr(y));    }};double det(const point &a,const point &b) {    return a.x*b.y-a.x*b.y;}double dot(const point &a,const point &b) {    return a.x*b.x+a.y*b.y;}double dis(const point &a,const point &b) {    return (a-b).norm();}point rotate_point(const point &p,double A) {    double tx=p.x,ty=p.y;    return point(tx*cos(A)-ty*sin(A),tx*sin(A)+ty*cos(A));}int dcmp(double k) {    return k<-eps?-1:k>eps?1:0;}double cross(const point &a,const point &b) {    return a.x*b.y-a.y*b.x;}double abs(const point &o) {    return sqrt(dot(o,o));}double mysqrt(double n){    return sqrt(max(0.0,n));}/**************************************/struct Circle {    point p;    double r;    Circle(){};    Circle(point p,double r):p(p),r(r){};    bool operator < (const Circle &o) const {        if(dcmp(r-o.r)!=0) return dcmp(r-o.r)==-1;        if(dcmp(p.x-o.p.x)!=0) {            return dcmp(p.x-o.p.x)==-1;        }        return dcmp(p.y-o.p.y)==-1;    }    bool operator == (const Circle &o) const {        return dcmp(r-o.r)==0&&dcmp(p.x-o.p.x)==0&&dcmp(p.y-o.p.y)==0;    }};//判断点是否在圆内  //误差范围内bool point_in_circle(point a,Circle p){    return dis(a,p.p)<p.r+eps;}//严格范围内bool point_in(const point &p,const Circle &a){    return cmp(dis(p,a.p)-a.r)<0;}/****************************************///圆与多边形的交面积int pon;            //多边形的点数point res[1000];  //逆/顺时针多边形的点double r;         //圆心在原点的圆的半径/***************************************///圆与线段/直线的交点//交点ret,个数numvoid circle_cross_line(point a,point b,Circle c,point ret[],int &num){    double x0=c.p.x,y0=c.p.y;    double x1=a.x,y1=a.y;    double x2=b.x,y2=b.y;    double dx=x2-x1,dy=y2-y1;    double A=dx*dx+dy*dy;    double B=2*dx*(x1-x0)+2*dy*(y1-y0);    double C=sqr(x1-x0)+sqr(y1-y0)-sqr(c.r);    double delta=B*B-4*A*C;    num=0;    if(dcmp(delta)>=0){        double t1=(-B-mysqrt(delta))/(2*A);        double t2=(-B+mysqrt(delta))/(2*A);        //圆与直线        //ret[num++]=point(x1+t1*dx,y1+t1*dy);        //ret[num++]=point(x1+t2*dx,y1+t2*dy);        //圆与线段        if(dcmp(t1-1)<=0&&dcmp(t1)>=0){ //            ret[num++]=point(x1+t1*dx,y1+t1*dy);        }        if(dcmp(t2-1)<=0&&dcmp(t2)>=0){            ret[num++]=point(x1+t2*dx,y1+t2*dy);        }    }    else if(dcmp(delta)==0){        double t1=(-B-mysqrt(delta))/(2*A);        //直线        //ret[num++]=point(x1+t1*dx,y1+t1*dy);        //线段        if(dcmp(t1-1)<=0&&dcmp(t1)>=0){ //            ret[num++]=point(x1+t1*dx,y1+t1*dy);        }    }}point crosspt(const point &a,const point &b,const point &p,const point &q){    double a1=cross(b-a,p-a);    double a2=cross(b-a,q-a);    return (p*a2-q*a1)/(a2-a1);}double sector_area(const point &a,const point &b){    double theta=atan2(a.y,a.x)-atan2(b.y,b.x);    while(theta<=0) theta+=2*pi;    while(theta>=2*pi) theta-=2*pi;    theta=min(theta,2*pi-theta);    return r*r*theta/2;}double calc(const point &a,const point &b){    point p[2];    int num=0;    int ina=dcmp(abs(a)-r)<0;    int inb=dcmp(abs(b)-r)<0;    if(ina){        if(inb) return fabs(cross(a,b))/2.0;        else{            circle_cross_line(a,b,Circle(point(0,0),r),p,num);            return sector_area(b,p[0])+fabs(cross(a,p[0]))/2.0;        }    }    else{        if(inb){            circle_cross_line(a,b,Circle(point(0,0),r),p,num);            return sector_area(p[0],a)+fabs(cross(p[0],b))/2.0;        }        else{            circle_cross_line(a,b,Circle(point(0,0),r),p,num);            if(num==2){                return sector_area(a,p[0])+sector_area(p[1],b)+fabs(cross(p[0],p[1]))/2.0;            }            else return sector_area(a,b);        }    }}//圆与多边形的交面积double circle_polygon_area(){    double ret=0;    for(int i=0;i<pon;i++){        int sgn=dcmp(cross(res[i],res[i+1]));        if(sgn!=0){            ret+=sgn*calc(res[i],res[i+1]);        }    }    return fabs(ret);}/****************************************///两个圆的交面积/几个圆的并面积/****************************************///两个圆的交面积double Circle_area(Circle x,Circle y){    double a=dis(x.p,y.p),b=x.r,c=y.r;    //相离、相切    if(a>=b+c) return 0.0;      //某个圆是点时    if(b<eps||c<eps) return 0.0;     //包含    if(b<c) swap(b,c);    if(a+c<=b) return pi*c*c;    //相交。    double cta1=acos((a*a+b*b-c*c)/2/(a*b)),           cta2=acos((a*a+c*c-b*b)/2/(a*c));    double s1=b*b*cta1-b*b*sin(cta1)*(a*a+b*b-c*c)/2/(a*b);    double s2=c*c*cta2-c*c*sin(cta2)*(a*a+c*c-b*b)/2/(a*c);    return fabs(s1+s2);}//圆的面积并。Circle tc[10]; //待求圆int cm;         //待求圆的个数//向量p旋转X角度后的向量，cost,sint为角度X的三角函数值point rotate(const point &p,double cost,double sint) {    double x=p.x,y=p.y;    return point(x*cost-y*sint,x*sint+y*cost);}//圆A,B的两个交点pair<point,point> crosspoint(point ap,double ar,point bp,double br) {    double d=(ap-bp).norm();    double cost=(ar*ar+d*d-br*br)/(2*ar*d);    double sint=sqrt(1.0-cost*cost);    point v=(bp-ap)/(bp-ap).norm()*ar;    return make_pair(ap+rotate(v,cost,-sint),ap+rotate(v,cost,sint));}inline pair<point ,point> crosspoint(const Circle &a,const Circle &b) {    return crosspoint(a.p,a.r,b.p,b.r);}struct Node {    point p;    double a;    int d;    Node(const point &p,double a,int d):p(p),a(a),d(d) {};    bool operator <(const Node &o)const {        return a<o.a;    }};double arg(point p) {    return arg(complex<double>(p.x,p.y));}double solve() {  //返回并面积    int n=0;    Circle c[10];    sort(tc,tc+cm);    cm=unique(tc,tc+cm)-tc;    for(int i=cm-1; i>=0; --i) {        bool ok=true;        for(int j=i+1; j<cm; j++) {            double d=(tc[i].p-tc[j].p).norm();            if(dcmp(d-abs(tc[i].r-tc[j].r))<=0) {                ok=false;                break;            }        }        if(ok) c[n++]=tc[i];    }    double ans=0;    for(int i=0; i<n; i++) {        vector<Node>event;        point boundary=c[i].p+point(-c[i].r,0);        event.push_back(Node(boundary,-pi,0));        event.push_back(Node(boundary,pi,0));        for(int j=0; j<n; ++j) {            if(i==j) continue;            double d=(c[i].p-c[j].p).norm();            if(dcmp(d-(c[i].r+c[j].r))<0) {                pair<point,point> ret=crosspoint(c[i],c[j]);                double x=arg(ret.first-c[i].p);                double y=arg(ret.second-c[i].p);                if(dcmp(x-y)>0) {                    event.push_back(Node(ret.first,x,1));                    event.push_back(Node(boundary,pi,-1));                    event.push_back(Node(boundary,-pi,1));                    event.push_back(Node(ret.second,y,-1));                } else {                    event.push_back(Node(ret.first,x,1));                    event.push_back(Node(ret.second,y,-1));                }            }        }        sort(event.begin(),event.end());        int sum=event[0].d;        for(int j=1; j<(int)event.size(); ++j) {            if(sum==0) {                ans+=cross(event[j-1].p,event[j].p)/2;                double x=event[j-1].a;                double y=event[j].a;                double area=c[i].r*c[i].r*(y-x)/2;                point v1=event[j-1].p-c[i].p;                point v2=event[j].p-c[i].p;                area-=cross(v1,v2)/2;                ans+=area;            }            sum+=event[j].d;        }    }    return ans;}/****************************************///圆的覆盖point p[533];  //初始点int pn;        //初始点的个数/****************************************///半径为r的圆最大能覆盖的点的个数 O(n^3)point center1,center2;//根据两个点确定的圆心//根据两个点,及半径确定圆心（有两个圆心）void get_center_point(point a,point b,double r){    double x0=(a.x+b.x)/2;    double y0=(a.y+b.y)/2;    double d=sqrt(r*r-((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y))/4);//圆心到中点的距离    if(fabs(a.y-b.y)<1e-6){        center1.x=x0;        center1.y=y0+d;        center2.x=x0;        center2.y=y0-d;    }    else{        double tmp=atan(-(a.x-b.x)/(a.y-b.y));        double dx=d*cos(tmp);        double dy=d*sin(tmp);        center1.x=x0+dx;        center1.y=y0+dy;        center2.x=x0-dx;        center2.y=y0-dy;    }}int max_point_cover() {    int ans=1;    //至少有1个点    for(int i=0;i<pn;i++){        for(int j=i+1;j<pn;j++){            if(dis(p[i],p[j])>2.0) continue;            get_center_point(p[i],p[j],1);            Circle temp1=Circle(center1,1);            Circle temp2=Circle(center2,1);            int ans1=0,ans2=0;            for(int k=0;k<pn;k++){                if(point_in_circle(p[k],temp1)) ans1++;                if(point_in_circle(p[k],temp2)) ans2++;            }            ans=max(ans,max(ans1,ans2));        }    }    return ans;}//给出n个点，求最小覆盖圆void circle_center(point p0,point p1,point p2,point &cp){    double a1=p1.x-p0.x,b1=p1.y-p0.y,c1=(a1*a1+b1*b1)/2;    double a2=p2.x-p0.x,b2=p2.y-p0.y,c2=(a2*a2+b2*b2)/2;    double d=a1*b2-a2*b1;    cp.x=p0.x+(c1*b2-c2*b1)/d;    cp.y=p0.y+(a1*c2-a2*c1)/d;}void circle_center(point p0,point p1,point &cp){    cp.x=(p0.x+p1.x)/2;    cp.y=(p0.y+p1.y)/2;}Circle min_circle_cover(){   // random_shuffle(p,p+pn);    Circle temp;    temp.p=p[0];    temp.r=0;    for(int i=1;i<pn;i++){        if(!point_in_circle(p[i],temp)){            temp.r=0;            temp.p=p[i];            for(int j=0;j<i;j++){                if(!point_in_circle(p[j],temp)){                    circle_center(p[i],p[j],temp.p);                    temp.r=dis(p[j],temp.p);                    for(int k=0;k<j;k++){                        if(!point_in_circle(p[k],temp)){                            circle_center(p[i],p[j],p[k],temp.p);                            temp.r=dis(p[k],temp.p);                        }                    }                }            }        }    }    return temp;}int main() {    return 0;}