二维计算几何模板--圆

转载请注明出处 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 0x3f3f3f3f

typedef 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,个数num
void 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;
}