陈家泽 计算几何模板

7 计 算 几 何
7.1 几 何 基 础
7.1.1 求 两 向 量 的 叉积

int cross(int x1,int y1,int x2,int y2) //int型
{
    return x1*y2-x2*y1;
}

double cross(double x1,double y1,double x2,double y2) //double型
{
    return x1*y2-x2*y1;
}

double cross(Point a,Point b) // 向量叉积
{
    return a.x*b.y-b.x*a.y;
}
View Code

7.1.2 求 平 面 两 点 欧 氏 距离

double distance(Point a,Point b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
View Code

7.1.5 判 断double 型 变 量 的 符号(ps:重点,替换kuangbin sgn函数)

const double eps=1e-8;
int dlcmp(double x)
{
    return x<-eps?-1:x>eps;
}
View Code

7.1.7 求 两 线 段 间 最 短 距离

double dis_segments(Segment seg1,Segment seg2)
{
    double m1=dis_point_segment(seg1.s,seg2);

    double m2=dis_point_segment(seg1.e,seg2);
    double m3=dis_point_segment(seg2.s,seg1);
    double m4=dis_point_segment(seg2.e,seg1);
    return min(min(m1,m2),min(m3,m4));
}
View Code

7.1.8 判 断 两 直 线 是 否 相 交 ( 共 线 及 平 行 , 若 相 交 则 求 出 交 点 )  (请认真考虑一下直线与线段,以及与kuangbin模板里的&操作符)

#include <math.h>
#define eps 1e-8
struct Point
{
    double x;
    double y;
};

struct Line
{
    Point s;
    Point e;
    Point v;
};
Line l1,l2;
double x,y交点坐标;//
int line_intersect(Line l1,Line l2)
{
    Point vec;
    double r;
    vec.x=l1.s.x-l2.s.x;
    vec.y=l1.s.y-l2.s.y;
    if (fabs(l1.v.x*l2.v.y-l2.v.x*l1.v.y)<eps)
    {
        if (fabs(l1.v.x*vec.y-vec.x*l1.v.y)<eps)
            return 2; //共线
        else
            return 0; //平行
    }
    else
    {
        r=((l1.s.x-l2.s.x)*l2.v.y-(l1.s.y-l2.s.y)*l2.v.x)
        /(l1.v.x*l2.v.y-l1.v.y*l2.v.x);
        x=-r*l1.v.x+l1.s.x;
        y=-r*l1.v.y+l1.s.y;
        return 1; //相交
    }
}
View Code

7.2 多 边 形
7.2.1 判 断 线 段 是 否 与 矩 形 相 交 ( 包 括 线 段 在 矩 形 内 部 )

struct Rectangle//这好像是矩形
{
    Point lt;//lefttop
    Point rb;//rightbottom
};

int segement_rectangle_intersect(Segment l,Rectangle r)
{
    Segment d1,d2;//retangle ’s diagonal

    d1.s=r.lt;
    d1.e=r.rb;
    d2.s.x=d1.e.x;
    d2.s.y=d1.s.y;
    d2.e.x=d1.s.x;
    d2.e.y=d1.e.y;

    if (l.s.x>=r.lt.x&&l.s.x<=r.rb.x&&
            l.s.y<=r.lt.y&&l.s.y>=r.rb.y||
            l.e.x>=r.lt.x&&l.e.x<=r.rb.x&&
            l.e.y<=r.lt.y&&l.e.y>=r.rb.y)

        return 1;

    if (segment_intersect(l,d1)|| segment_intersect(l,d2))
        //segment_intersect(endpoint inclusive)
        return 1;

    return 0;
}
View Code

7.2.3 求 简 单 多 边 形 重 心

Point get_center(Point pt[],int n)
{
    double sum,area;
    Point res(0,0),o(0,0);
    int i;
    sum=0;
    for (i=0; i<n; i++)
    {
        area=cross(pt[i]-o,pt[(i+1)%n]-o);
        res=res+(pt[i]+pt[(i+1)%n])/3*area;
        sum+=area;
    }
    res=res/sum;
    return res;
}
View Code

7.2.7 判 断 两 凸 多 边 形 是 否 相 交 (graham scan 求 凸 包+ 枚 举 边 、 点)

#include <iostream >
#include <cstdio >
#include <algorithm >
#define N 110
using namespace std;

struct Point
{
    int x,y;
};

struct Polygon
{
    Point p[N];
    int n;
};

Point pt[N];
int stack[N];

int cross(Point a,Point b,Point s)
{
    return (a.x-s.x)*(b.y-s.y)-(b.x-s.x)*(a.y-s.y);
}

int dist(Point a,Point b)
{
    return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);
}

int cmp(Point a,Point b)
{
    if (cross(a,b,pt[0]) >0||cross(a,b,pt[0])==0&&dist(a,pt[0])<dist(b,pt[0]))
        return 1;
    else
        return 0;
}

int on_segment(Point s,Point e,Point o)
{

    if(cross(s,e,o)==0&&
       o.x>=min(s.x,e.x)&&
       o.x<=max(s.x,e.x)&&
       o.y>=min(s.y,e.y)&&
       o.y<=max(s.y,e.y))

        return 1;
    else
        return 0;
}

int graham_scan(int n)
{
    int i,top,t;
    if (n<=1)
    {
        stack[0]=0;
        return n;
    }
    for (t=0,i=1; i<n; i++)
        if (pt[i].y<pt[t].y||pt[i].y==pt[t].y&&pt[i].x<pt[t].x)
            t=i;
    swap(pt[0],pt[t]);
    sort(pt+1,pt+n,cmp);
    top=2;
    for (i=0; i<2; i++)
        stack[i]=i;
    for (i=2; i<n; i++)
    {
        while (top >1&&cross(pt[stack[top -1]],pt[i],pt[stack[top -2]]) <=0)
            top--;
        stack[top++]=i;
    }

    return top;
}

int segment_intersect(Point s1,Point e1,Point s2,Point e2)
{
    if (max(s1.x,e1.x)>=min(s2.x,e2.x)&&
        max(s2.x,e2.x)>=min(s1.x,e1.x)&&
        max(s1.y,e1.y)>=min(s2.y,e2.y)&&
        max(s2.y,e2.y)>=min(s1.y,e2.y)&&
        (double)cross(s2,e1,s1)*(double)cross(e2,e1,s1)<=0&&
        (double)cross(s1,e2,s2)*(double)cross(e1,e2,s2)<=0)

        return 1;
    else
        return 0;
}

int point_inside(Point o,Polygon pln)
{
    int i,a1=0,a2=0,n=pln.n;

    pln.p[n]=pln.p[0];

    for (i=0; i<n; i++)
        a1+=abs(cross(pln.p[i],pln.p[i+1],o));

    for (i=1; i<n; i++)
        a2+=abs(cross(pln.p[i],pln.p[i+1],pln.p[0]));

    if (a1==a2)
        return 1;
    else
        return 0;
}

int convex_polygon_intersect(Polygon pln1,Polygon pln2)
{
    int i,j;
    pln1.p[pln1.n]=pln1.p[0];
    pln2.p[pln2.n]=pln2.p[0];

    for (i=0; i<pln1.n; i++)
        for (j=0; j<pln2.n; j++)
            if (segment_intersect(pln1.p[i],pln1.p[i+1],pln2.p[j],pln2.p[j+1]))
                return 1;

    if (point_inside(pln1.p[0],pln2)||point_inside(pln2.p[0],pln1))
        return 1;

    return 0;
}

int main()
{
    int n,m,i,vertexnum,ans;
    Polygon pln1,pln2;

    while (cin>>n>>m,n||m)
    {
        for (i=0; i<n; i++)
            cin>>pt[i].x>>pt[i].y;
        vertexnum=graham_scan(n);
        pln1.n=vertexnum;

        for (i=0; i<vertexnum; i++)
            pln1.p[i]=pt[stack[i]];

        for (i=0; i<m; i++)
            cin>>pt[i].x>>pt[i].y;

        vertexnum=graham_scan(m);
        pln2.n=vertexnum;

        for (i=0; i<vertexnum; i++)
            pln2.p[i]=pt[stack[i]];

        if (pln1.n==1&&pln2.n==1)
            ans=1;
        else if (pln1.n==1&&pln2.n==2)
        {
            if (on_segment(pln2.p[0],pln2.p[1],pln1.p[0]))
                ans=0;
            else
                ans=1;
        }
        else if (pln1.n==2&&pln2.n==1)
        {
            if (on_segment(pln1.p[0],pln1.p[1],pln2.p[0]))
                ans=0;
            else
                ans=1;
        }
        else if (pln1.n==2&&pln2.n==2)
        {
            if (segment_intersect(pln1.p[0],pln1.p[1],pln2.p[0],pln2.p[1]))
                ans=0;
            else
                ans=1;
        }
        else if (pln1.n==1)
        {
            if (point_inside(pln1.p[0],pln2))
                ans=0;
            else
                ans=1;
        }
        else if (pln2.n==1)
        {
            if (point_inside(pln2.p[0],pln1))
                ans=0;
            else
                ans=1;
        }
        else
        {
            if (convex_polygon_intersect(pln1,pln2)==0)
                ans=1;
            else
                ans=0;
        }

        if (ans==0)
            cout <<"NO"<<endl;
        else
            cout <<"YES"<<endl;
    }

    return 0;
}
View Code

7.3.3 两 个 简 单 多 边 形 求 面 积 并 、交

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#define eps 1e-8
#define N 550
using namespace std;

struct Point
{
    double x,y;

    Point () {}
    Point (double a,double b):x(a),y(b) {}

    Point operator - (const Point a) const
    {
        return Point(x-a.x,y-a.y);
    }
};

Point zero(0,0);
int dlcmp(double x)
{
    return x<-eps?-1:x>eps;
}
double cross(Point v1,Point v2)
{
    return v1.x*v2.y-v2.x*v1.y;
}

struct Polygon
{
    Point p[N];
    int n;

    Polygon():n(0) {}
    void clear()
    {
        n=0;
    }
    void add(Point a)
    {
        p[n++]=a;
    }

    double area()
    {
        double res=0;
        for (int i=1; i<n-1; i++)
            res+=cross(p[i]-p[0],p[i+1]-p[0]);
        return fabs(res/2);
    }
};

Polygon A,B,rec;


Point line_intersect(Point s1,Point e1,Point s2,Point e2) //两直线交点
{
    Point v1,v2,res;
    double r;

    v1=s1-e1;
    v2=s2-e2;
    r=((s1.x-s2.x)*v2.y-(s1.y-s2.y)*v2.x)/(v1.x*v2.y-v1.y*v2.x);
    res.x=-r*v1.x+s1.x;
    res.y=-r*v1.y+s1.y;

    return res;
}

void cut(Point s,Point e) //半平面交
{
    int i,j,d1,d2;
    Polygon ker;

    for (i=0; i<rec.n; i++)
    {
        j=(i+1)%rec.n;
        d1=dlcmp(cross(e-s,rec.p[i]-s));
        d2=dlcmp(cross(e-s,rec.p[j]-s));

        if (d1>=0)
            ker.add(rec.p[i]);
        if (d1*d2<0)
            ker.add(line_intersect(s,e,rec.p[i],rec.p[j]));
    }

    rec=ker;
}

double calc(Point p1,Point p2,Point q1,Point q2)
{
    int dp=dlcmp(cross(p1,p2)),dq=dlcmp(cross(q1,q2));
    int sgn=dp*dq;

    if (sgn==0)
        return 0;

    rec.clear();
    rec.add(zero);
    rec.add(p1);
    rec.add(p2);
    if (dp<0)
        swap(rec.p[1],rec.p[2]);
    if (dq>0)
    {
        cut(zero,q1);
        cut(q1,q2);
        cut(q2,zero);
    }
    else
    {
        cut(zero,q2);
        cut(q2,q1);
        cut(q1,zero);
    }

    return sgn*rec.area();
}

double solve()
{
    double res=A.area()+B.area();
    double sum=0;
    int i,j;

//对两个多边形三角剖分,分别求两个三角形的面积交
    for (i=0; i<A.n; i++)
        for (j=0; j<B.n; j++)
            sum+=calc(A.p[i],A.p[(i+1)%A.n],B.p[j],B.p[(j+1)%B.n]);
    res-=fabs(sum); //fabs(sum)为两个多边形的面积交

    return res; //面积并
}
//hdu3060(题目数据有误)
int main()
{
    int n,m,i;
    Point pt;

    double ans;

    while (scanf("%d%d",&n,&m)!=EOF)
    {
        A.clear();
        B.clear();
        for (i=0; i<n; i++)
        {
            scanf("%lf%lf",&pt.x,&pt.y);
            A.add(pt);
        }
        for (i=0; i<m; i++)
        {
            scanf("%lf%lf",&pt.x,&pt.y);
            B.add(pt);
        }

        ans=solve();

        printf("%.2f\n",ans+eps);
    }

    return 0;
}
View Code

 

 

7.4 圆
7.4.1 点类

struct Point
{
    double x,y;

    Point() {}
    Point(double a,double b):x(a),y(b) {}
    Point operator + (const Point a) const
    {
        return Point(x+a.x,y+a.y);
    }
    Point operator - (const Point a) const
    {
        return Point(x-a.x,y-a.y);
    }
    Point operator * (const double a) const
    {
        return Point(x*a,y*a);
    }
    Point operator / (const double a) const
    {
        return Point(x/a,y/a);
    }

    bool operator < (const Point a) const
    {
        if (dlcmp(x-a.x)==0)
            return dlcmp(x-a.y)<0;
        else
            return dlcmp(x-a.x)<0;
    }
    bool operator == (const Point a) const
    {
        return !dlcmp(x-a.x)&&!dlcmp(y-a.y);
    }

//向量长度定为d
    Point trunc(double d)
    {
        double dis(Point,Point);
        double len=dis(*this,Point(0,0));
        return Point(x*d/len,y*d/len);
    }

//坐标逆时针旋转a度
    Point rotate(double a)
    {
        return Point(x*cos(a)-y*sin(a),y*cos(a)+x*sin(a));
    }
};

double dis(Point a,Point b)
{
    return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y));
}

double cross(Point a,Point b,Point s)
{
    double x1=a.x-s.x,y1=a.y-s.y;
    double x2=b.x-s.x,y2=b.y-s.y;

    return x1*y2-x2*y1;
}

double cross(Point a,Point b)
{
    return a.x*b.y-b.x*a.y;
}

double dot(Point a,Point b,Point s)
{
    double x1=a.x-s.x,y1=a.y-s.y;
    double x2=b.x-s.x,y2=b.y-s.y;

    return x1*x2+y1*y2;
}

double dot(Point a,Point b)
{
    return a.x*b.x+a.y*b.y;
}
View Code

 7.4.2 圆 类

struct Circle
{
    Point o;
    double r;
    Circle() {}
    Circle(Point a,double l):o(a),r(l) {}

    double area()
    {
        return sqr(r)*PI;
    }
};

//判断圆a是否含于圆b
int inner_circle(Circle a,Circle b)
{
    if (dlcmp(a.r-b.r)>0)
        return 0;
    return dlcmp(dis(a.o,b.o)+a.r-b.r)<=0;
}

//以base点为基点,极角排序,排序前base需赋初值
Point base;
int cmp(const Point a,const Point b)
{
    return atan2(a.y-base.y,a.x-base.x)<atan2(b.y-base.y,b.x-base.x);
}

//向量a,b的夹角
double vec_angle(Point a,Point b)
{
    double tmp=dot(a,b)/(dis(a,Point(0,0))*dis(b,Point(0,0)));
    if (dlcmp(tmp -1) >=0) tmp=1;
    if (dlcmp(tmp+1) <=0) tmp=-1;

    return acos(tmp);
}

//计算由a到b逆时针方向的弓形面积
double arc_area(Point a,Point b,Circle c)
{
    double theta=vec_angle(a-c.o,b-c.o);
    double sf=sqr(c.r)*theta/2.0;
    double st=sqr(c.r)*sin(theta)/2.0;

    if (dlcmp(cross(a,b,c.o))>0)
        return sf-st;
    else
        return c.area()-sf+st;
}

double arc_area(double th,double r)
{
    return 0.5*sqr(r)*(th-sin(th));
}
View Code

7.4.3 圆 面 积 交 、并

//求两圆交点,排除相切的情况,不考虑内含
int inter_circle_or(Circle c1,Circle c2,Point &p1,Point &p2)
{
    double len=dis(c1.o,c2.o);

    if (dlcmp(len-c1.r-c2.r)>=0)
        return 0;

    double s=(sqr(c1.r)-sqr(c2.r)+sqr(len))/len/2;
    double h=sqrt(sqr(c1.r)-sqr(s));
    Point vec=c2.o-c1.o;
    Point p0=c1.o+vec.trunc(s);
    //圆类里的函数,可能和kuangbin相同
    p1=p0+vec.rotate(PI/2).trunc(h);
    p2=p0-vec.rotate(PI/2).trunc(h);
    return 1;
}
View Code

 7.4.4 简 单 多 边 形 与 圆 求 面 积 交

 #include <cstdio>
 #include <cstring>
 #include <algorithm>
 #include <cmath>
 #define N 200
 #define eps 1e-8
 using namespace std;

 const double PI=acos(-1.0);

 struct Point
 {
 double x,y;
 };

 Point pt[N];
 int n;

 int dlcmp(double x)
 {
 return x<-eps?-1:x>eps;
 }

 double sqr(double x)
 {
 return x*x;
 }

 double dis(Point a,Point b)
 {
 return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y));
 }

 double outer(Point a,Point b,Point c)
 {
 return (a.x-c.x)*(b.y-c.y)-(a.y-c.y)*(b.x-c.x);
 }

 double inner(Point a,Point b,Point c)
 {
 return (a.x-c.x)*(b.x-c.x)+(a.y-c.y)*(b.y-c.y);
 }

 double calc_area(Point a,Point b,Point c,double r)
 {
 double A,B,C,x,y,tS;

 A=dis(b,c);
 B=dis(a,c);
 C=dis(b,a);

 if (A<r&&B<r)
 return outer(a,b,c)/2;

 else if (A<r&&B>=r)
 {
 x=(inner(a,c,b)+sqrt(sqr(r)*sqr(C)-sqr(outer(a,c,b))))/C;
 tS=outer(a,b,c)/2;

 return asin(tS*(1-x/C)*2/r/B)*sqr(r)/2+tS*x/C;
 }
 else if (A>=r&&B<r)
 {
 y=(inner(b,c,a)+sqrt(sqr(r)*sqr(C)-sqr(outer(b,c,a))))/C;
 tS=outer(a,b,c)/2;

 return asin(tS*(1-y/C)*2/r/A)*sqr(r)/2+tS*y/C;
 }
 else if (fabs(outer(a,b,c))>=r*C||inner(b,c,a) <=0||inner(a,c,b)<=0)
 {
 if (inner(a,b,c)<0)
 {
 if (outer(a,b,c)<0)
 return (-PI-asin(outer(a,b,c)/A/B))*sqr(r)/2;
 else
 return (PI-asin(outer(a,b,c)/A/B))*sqr(r)/2;
 }
 else
 return asin(outer(a,b,c)/A/B)*sqr(r)/2;
 }
 else
 {
 x=(inner(a,c,b)+sqrt(sqr(r)*sqr(C)-sqr(outer(a,c,b))))/C;
 y=(inner(b,c,a)+sqrt(sqr(r)*sqr(C)-sqr(outer(b,c,a))))/C;
 tS=outer(a,b,c)/2;

 return (asin(tS*(1-x/C)*2/r/B)+asin(tS*(1-y/C)*2/r/A))*sqr(r)/2+tS*((y+x)/C-1);
 }
 }

 //计算一般多边形与圆的交面积(将多边形划分为三角形,然后有向三角形与圆求有向面积交)
 double solve(Point o,double r)
 {
 int i,j;
double res,sum;
 Point tri[3];

 res=0;
 for (i=1;i<n-1;i++)
 {
 tri[0]=pt[0];
 tri[1]=pt[i];
 tri[2]=pt[i+1];
 sum=0;

 for (j=0;j<3;j++)
 sum+=calc_area(tri[j],tri[(j+1)%3],o,r);

 //sum为三角形与圆交的有向面积
 res+=sum;
 }

 return fabs(res);
 }

 //poj3675
 int main()
 {
 double x0,y0,v,vx,vy,g,r,t,theta ,ans;
 Point o;
 int i;

 o.x=o.y=0;

 while (scanf("%lf",&r)!=EOF)
 {
 scanf("%d",&n);
 for (i=0;i<n;i++)
 scanf("%lf%lf",&pt[i].x,&pt[i].y);

 ans=solve(o,r);

 printf("%.2f\n",ans);
 }

 return 0;
 }
View Code

7.4.5 求 线 段 与 圆 的 交 点
若求直线与圆的交点类似,无需讨论k1、k2的取值范围

int inter_circle_segment(Circle c,Point a,Point b,Point &p1,Point &p2)
{
    Point vec=b-a;
    double A=sqr(vec.x)+sqr(vec.y);
    double B=2*(vec.x*(a.x-c.o.x)+vec.y*(a.y-c.o.y));
    double C=sqr(a.x-c.o.x)+sqr(a.y-c.o.y)-sqr(c.r);
    double delta=sqr(B)-4*A*C;

    if (dlcmp(delta)<0)
        return 0;

    double k1=(-B-sqrt(fabs(delta)))/(2*A);
    double k2=(-B+sqrt(fabs(delta)))/(2*A);
    int res=0;

    if (dlcmp(k1) >=0&&dlcmp(k1-1) <=0)
    {
        res++;
        p1=a+vec*k1;
    }
    if (dlcmp(k2) >=0&&dlcmp(k2-1) <=0)
    {
        res++;
        if (res==1)
            p1=a+vec*k2;
        else
            p2=a+vec*k2;
    }
    return res;
}
View Code

7.4.6 求两圆公切线

//求两相离的圆的两条内共切线
void get_InCommonTangent(Circle c1,Circle c2,Point &s1,Point &e1,Point &s2,Point &e2)
{
    double l=dis(c1.o,c2.o);
    double d=l*c1.r/(c1.r+c2.r);
    double tmp=c1.r/d;
    tmp=fix(tmp);
    double theta=acos(tmp);
    Point vec=c2.o-c1.o;

    vec=vec.trunc(c1.r);
    s1=c1.o+vec.rotate(theta);
    s2=c1.o+vec.rotate(-theta);

    vec=c1.o-c2.o;
    vec=vec.trunc(c2.r);
    e1=c2.o+vec.rotate(theta);
    e2=c2.o+vec.rotate(-theta);
}

//求两相离的圆的两条外公切线
void get_OutCommonTangent(Circle c1,Circle c2,Point &s1,Point &e1,Point &s2,Point &e2)
{
    double l=dis(c1.o,c2.o);
    double d=fabs(c1.r-c2.r);
    double theta=acos(d/l);

    if (dlcmp(c1.r-c2.r)>0)
        swap(c1,c2);

    Point vec=c1.o-c2.o;
    vec=vec.trunc(c1.r);
    s1=c1.o+vec.rotate(theta);
    s2=c1.o+vec.rotate(-theta);
    vec=vec.trunc(c2.r);
    e1=c2.o+vec.rotate(theta);
    e2=c2.o+vec.rotate(-theta);
}
View Code

7.4.7 最 小 圆 覆盖(平面上n个点,求一个半径最小的圆,能够覆盖所有的点。)

#include <Bits/stdc++.h>
using namespace std;
#define eps 1e-8
#define MAX_P 2000
struct Point
{
    double x,y;

    Point operator -(Point &a)
    {
        Point t;
        t.x=x-a.x;
        t.y=y-a.y;
        return t;
    }
};
struct Circle
{
    double r;
    Point center;
};

struct Triangle
{
    Point t[3];
};

Point pt[MAX_P]; //点集
Circle c; //最小圆

double distance(Point a,Point b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

double cross(Point a,Point b)
{
    return a.x*b.y-b.x*a.y;
}

double triangle_area(Triangle tri) //三角形距离
{
    Point v1=tri.t[1]-tri.t[0];
    Point v2=tri.t[2]-tri.t[0];

    return fabs(cross(v1,v2))/2;
}

Circle circumcircle_triangle(Triangle tri) //三角形外接圆
{
    Circle res;

    double a,b,c,c1,c2;
    double xA,yA,xB,yB,xC,yC;

    a=distance(tri.t[0],tri.t[1]);
    b=distance(tri.t[1],tri.t[2]);
    c=distance(tri.t[2],tri.t[0]);

    res.r=a*b*c/triangle_area(tri)/4;

    xA=tri.t[0].x;
    yA=tri.t[0].y;
    xB=tri.t[1].x;
    yB=tri.t[1].y;
    xC=tri.t[2].x;
    yC=tri.t[2].y;

    c1=(xA*xA+yA*yA-xB*xB-yB*yB)/2;
    c2 = (xA*xA+yA*yA-xC*xC-yC*yC)/2;

    res.center.x=(c1*(yA-yC)-c2*(yA-yB))/ ((xA-xB)*(yA-yC)-(xA-xC)*(yA-yB));
    res.center.y = (c1*(xA-xC)-c2*(xA-xB))/ ((yA-yB)*(xA-xC)-(yA-yC)*(xA-xB));

    return res;
}

Circle mincircle_triangle(int trinum,Triangle tri)
{
    Circle res;

    if (trinum==0)
        res.r=-2;
    else if (trinum==1)
    {
        res.center=tri.t[0];
        res.r=0;
    }
    else if (trinum==2)
    {
        res.center.x=(tri.t[0].x+tri.t[1].x)/2;
        res.center.y=(tri.t[0].y+tri.t[1].y)/2;
        res.r=distance(tri.t[0],tri.t[1])/2;
    }
    else if (trinum==3)
        res=circumcircle_triangle(tri);
    return res;
}

void mincircle_pointset(int m,int trinum,Triangle tri)  //求点集的最小覆盖圆
{
    int i,j;
    Point tmp;

    c=mincircle_triangle(trinum,tri);

    if (trinum==3)
        return;

    for (i=0; i<m; i++)
        if (distance(pt[i],c.center)>c.r)
        {
            tri.t[trinum]=pt[i];
            mincircle_pointset(i,trinum+1,tri);
            tmp=pt[i];

            for (j=i; j>=1; j--)
                pt[j]=pt[j-1];

            pt[0]=tmp;
        }
}

int main()
{
    int n,i,f1,f2;
    Triangle tri;
    while (scanf("%d%d%d",&f1,&f2,&n)!=EOF)
    {
        for (i=0; i<n; i++)
            scanf("%lf%lf",&pt[i].x,&pt[i].y);

        mincircle_pointset(n,0,tri);
        printf("%lf %lf %lf\n",c.center.x,c.center.y,c.r);
    }
    return 0;
}
View Code

 7.4.8 单 位 圆 覆盖(一个单位圆最多能覆盖平面上多少点)

#include <math.h>
#include <bits/stdc++.h>
using namespace std;
#define eps 1e-8
#define MAX_P 505
const double r=1.0;//单位圆半径

struct Point
{
    double x,y;
};

Point pt[MAX_P];

double distance(Point a,Point b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
//sqrt函数速度较慢,应尽量避免出现,此处可优化为距离的平方和的形式
}

Point find_center(Point a,Point b)
{
    Point v,mid,center;
    double d,s,ang;

    v.x=a.x-b.x;
    v.y=a.y-b.y;

    mid.x=(a.x+b.x)/2;
    mid.y=(a.y+b.y)/2;

    d=distance(a,mid);
    s=sqrt(r*r-d*d); //优化为s=sqrt(r*r-d);

    if (fabs(v.y)<eps)
    {
        center.x=mid.x;
        center.y=mid.y+s;
    }
    else
    {
        ang=atan(-v.x/v.y);
        center.x=mid.x+s*cos(ang);
        center.y=mid.y+s*sin(ang);
    }
    return center;
}

int main()
{

    int n,i,j,k,ans,cnt;
    double tmp;
    Point center;

    while (~scanf("%d",&n),n)
    {
        for (i=0; i<n; i++)
            scanf("%lf%lf",&pt[i].x,&pt[i].y);
        ans=1;
        for (i=0; i<n; i++)
            for (j=i+1; j<n; j++)
            {
                if (distance(pt[i],pt[j])>2*r) //优化为distance(pt[i],pt[j])>2*2*r*r
                    continue;
                cnt=0;
                center=find_center(pt[i],pt[j]);

                for (k=0; k<n; k++)
                    if (distance(pt[k],center)<=r+eps)
                        cnt++;
                if (ans<cnt)
                    ans=cnt;
            }
        printf("%d\n",ans);
    }

    return 0;
}
View Code

 

7.5 模 拟 退 火
7.5.1 求 多 边 形 费 马点

#include <iostream >
 #include <cstdio >
 #include <cmath >
 #define eps 1e-6
 #define N 105
 using namespace std;

 struct Point
 {
 double x,y;
 };


 double point_dis(Point a,Point b)
 {
 return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));

 }

 double sum_dis(Point pt[],int n,Point o)
 {
 double res=0;
 int i;

 for (i=0;i<n;i++)
 res+=point_dis(pt[i],o);

 return res;
 }

 double polygon_Fermatpoint(Point pln[],int n)
{
 Point cp,np,tmp;
 double min,step,d;
 int flag;

 cp=pln[0]; //cp保存当前更新后最优的费马点
 min=sum_dis(pln,n,cp);
 step=10000; //选取坐标范围的最大值

 while (step>eps)
 {
 flag=1;
 while (flag)
 {
 flag=0;
 np=cp;

 tmp=cp,tmp.x+=step;
 d=sum_dis(pln,n,tmp);

 if (min>d)
 min=d, np=tmp, flag=1;

 tmp=cp,tmp.x-=step;
 d=sum_dis(pln,n,tmp);

 if (min>d)

 min=d, np=tmp,flag=1;

 tmp=cp,tmp.y+=step;
 d=sum_dis(pln,n,tmp);

 if (min>d)
 min=d, np=tmp,flag=1;

 tmp=cp,tmp.y-=step;
 d=sum_dis(pln,n,tmp);

 if (min>d)
 min=d, np=tmp,flag=1;

 cp=np;
 }

 step*=0.98; //系数根据精度要求修改
 }

 return min;
 }

int main()
 {
 int n,i;
 double min;
 Point pt[N];

 cin>>n;

 for (i=0;i<n;i++)
 cin>>pt[i].x>>pt[i].y;

 min=polygon_Fermatpoint(pt,n);

 printf("%.0f\n",min);

 return 0;
 }
View Code

7.5.2 最 小 球 覆盖(求能覆盖所有点的最小球的半径。)

#include <iostream >
#include <cstdio >
#include <cmath >
#define oo 1e20
#define eps 1e-10
#define N 105
using namespace std;

struct Point
{
    double x,y,z;
};

double dis(Point a,Point b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)+(a.z-b.z)*(a.z-b.z));
}

int max_dis(Point pt[],int n,Point o)
{
    int i,res;
    double max,tmp;
    max=0;
    res=0;
    for (i=0; i<n; i++)
    {
        tmp=dis(pt[i],o);
        if (max<tmp)
        {
            max=tmp;
            res=i;
        }
    }
    return res;
}

int main()
{
    Point pt[N],o;
    int n,i,t;
    double dx,dy,dz,step,r,tmp;
    cin>>n;
    for (i=0; i<n; i++)
        cin>>pt[i].x>>pt[i].y>>pt[i].z;
    step=10000; //step选取最大的坐标范围
    r=oo;

    if (n==1)
    {
        o.x=pt[0].x;
        o.y=pt[0].y;
        o.z=pt[0].z;
    }
    else
    {
        o.x=o.y=o.z=0;
        while (step>eps)
        {
            t=max_dis(pt,n,o);
            tmp=dis(pt[t],o);

            if (r>tmp) r=tmp;

            dx=(pt[t].x-o.x)/tmp;
            dy=(pt[t].y-o.y)/tmp;
            dz=(pt[t].z-o.z)/tmp;

            o.x+=step*dx;
            o.y+=step*dy;
            o.z+=step*dz;

            step*=0.9993; //系数的选取根据具体精度调整

        }
    }
    printf("%.6f %.6f %.6f\n",o.x,o.y,o.z);
    return 0;
}
View Code

 

posted @ 2017-10-19 22:38  啦啦啦天啦噜  阅读(292)  评论(0编辑  收藏  举报