码迷,mamicode.com
首页 > 其他好文 > 详细

计算几何 综合模板

时间:2017-10-20 10:21:21      阅读:197      评论:0      收藏:0      [点我收藏+]

标签:ret   int()   cross   print   i++   end   通过   ***   ios   

1、基本函数
  1.1 Point 定义

const double eps=1e-8;
const double PI = acos(-1.0);
int sgn(double x)
{
    return x<-eps?-1:x>eps;
}
struct point
{
    double x,y;
    point() {}
    point(double _x,double _y)
    {
        x = _x;
        y = _y;
    }
    point operator -(const point &b)const
    {
        return point(x - b.x,y - b.y);
    }
//叉积
    double operator ^(const point &b)const
    {
        return x*b.y - y*b.x;
    }
//点积
    double operator *(const point &b)const
    {
        return x*b.x + y*b.y;
    }
//绕原点旋转角度B(弧度值),后x,y的变化
    void transXY(double B)
    {
        double tx = x,ty = y;
        x = tx*cos(B) - ty*sin(B);
        y = tx*sin(B) + ty*cos(B);
    }
};

1.2 Line 定义

struct Line
{
    point s,e;
    Line() {}
    Line(point _s,point _e)
    {
        s = _s;
        e = _e;
    }
//两直线相交求交点
//第一个值为0表示直线重合,为1表示平行,为0表示相交,为2是相交
//只有第一个值为2时,交点才有意义
    pair<int,point> operator &(const Line &b)const
    {
        point res = s;
        if(sgn((s-e)^(b.s-b.e)) == 0)
        {
            if(sgn((s-b.e)^(b.s-b.e)) == 0)
                return make_pair(0,res);//重合
            else return make_pair(1,res);//平行
        }
        double t = ((s-b.s)^(b.s-b.e))/((s-e)^(b.s-b.e));
        res.x += (e.x-s.x)*t;
        res.y += (e.y-s.y)*t;
        return make_pair(2,res);
    }
};

1.4 判断:线段相交

bool inter(Line l1,Line l2)
{
    return
        max(l1.s.x,l1.e.x) >= min(l2.s.x,l2.e.x) &&
        max(l2.s.x,l2.e.x) >= min(l1.s.x,l1.e.x) &&
        max(l1.s.y,l1.e.y) >= min(l2.s.y,l2.e.y) &&
        max(l2.s.y,l2.e.y) >= min(l1.s.y,l1.e.y) &&
        sgn((l2.s-l1.e)^(l1.s-l1.e))*sgn((l2.e-l1.e)^(l1.s-l1.e)) <= 0 &&
        sgn((l1.s-l2.e)^(l2.s-l2.e))*sgn((l1.e-l2.e)^(l2.s-l2.e)) <= 0;
}

1.5 判断:直线和线段相交

//判断直线和线段相交
bool Seg_inter_line(Line l1,Line l2) //判断直线l1和线段l2是否相交
{
    return sgn((l2.s-l1.e)^(l1.s-l1.e))*sgn((l2.e-l1.e)^(l1.s-l1.e)) <= 0;
}

1.6 点到直线距离

//点到直线距离
//返回为result,是点到直线最近的点
Point PointToLine(Point P,Line L)
{
    Point result;
    double t = ((P-L.s)*(L.e-L.s))/((L.e-L.s)*(L.e-L.s));
    result.x = L.s.x + (L.e.x-L.s.x)*t;
    result.y = L.s.y + (L.e.y-L.s.y)*t;
    return result;
}

1.7 点到线段距离

Point NearestPointToLineSeg(Point P,Line L)
{
    Point result;
    double t = ((P-L.s)*(L.e-L.s))/((L.e-L.s)*(L.e-L.s));
    if(t >= 0 && t <= 1)
    {
        result.x = L.s.x + (L.e.x - L.s.x)*t;
        result.y = L.s.y + (L.e.y - L.s.y)*t;
    }
    else
    {
        if(dist(P,L.s) < dist(P,L.e))
            result = L.s;
        else result = L.e;
    }
    return result;
}

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));
}

 

1.8 计算多边形面积

//计算多边形面积
//点的编号从0~n-1
double CalcArea(Point p[],int n)
{
    double res = 0;
    for(int i = 0; i < n; i++)
        res += (p[i]^p[(i+1)%n])/2;
    return fabs(res);
}

1.9 判断点在线段上

//*判断点在线段上
bool OnSeg(point P,Line L)
{
    return
        sgn((L.s-P)^(L.e-P)) == 0 &&
        sgn((P.x - L.s.x) * (P.x - L.e.x)) <= 0 &&
        sgn((P.y - L.s.y) * (P.y - L.e.y)) <= 0;
}

1.11 判断点在任意多边形内

//*判断点在任意多边形内
//射线法,poly[]的顶点数要大于等于3,点的编号0~n-1
//返回值
//-1:点在凸多边形外
//0:点在凸多边形边界上
//1:点在凸多边形内
int inPoly(point p,point poly[],int n)
{
    int cnt;
    Line ray,side;
    cnt = 0;
    ray.s = p;
    ray.e.y = p.y;
    ray.e.x = -100000000000.0;//-INF,注意取值防止越界
    for(int i = 0; i < n; i++)
    {
        side.s = poly[i];
        side.e = poly[(i+1)%n];
        if(OnSeg(p,side))return 0;
//如果平行轴则不考虑
        if(sgn(side.s.y - side.e.y) == 0)
            continue;
        if(OnSeg(side.s,ray))
        {
            if(sgn(side.s.y - side.e.y) > 0)cnt++;
        }
        else if(OnSeg(side.e,ray))
        {
            if(sgn(side.e.y - side.s.y) > 0)cnt++;
        }
        else if(inter(ray,side))
            cnt++;
    }
    if(cnt % 2 == 1)return 1;
    else return -1;
}

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;
}

 

1.12 判断凸多边形

//判断凸多边形
//允许共线边
//点可以是顺时针给出也可以是逆时针给出
//点的编号1~n-1
bool isconvex(Point poly[],int n)
{
    bool s[3];
    memset(s,false,sizeof(s));
    for(int i = 0; i < n; i++)
    {
        s[sgn( (poly[(i+1)%n]-poly[i])^(poly[(i+2)%n]-poly[i]) )+1] = true;
        if(s[0] && s[2])return false;
    }
    return true;
}

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;
}

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;
}

 

简单的极角排序(以第一个点为基准点)

const int maxn=55;
point List[maxn];
double dist(point p0,point p1)
{
    return (double)sqrt((p1.x-p0.x)*(p1.x-p0.x)+(p1.y-p0.y)*(p1.y-p0.y));
}
bool _cmp(point p1,point p2)
{
    double tmp = (p1-List[0])^(p2-List[0]);
    if(sgn(tmp) > 0)return true;
    else if(sgn(tmp) == 0 && sgn(dist(p1,List[0]) - dist(p2,List[0])) <= 0)
        return true;
    else return false;
}

sort(List+1,List+n,_cmp);

2、凸包(ps:这个算法前面的的排序是先找到最左下角的一个点,算法复杂度O(nlogn)

/*
* 求凸包,Graham算法
* 点的编号0~n-1
* 返回凸包结果Stack[0~top-1]为凸包的编号
*/
const int MAXN = 1010;
point List[MAXN];
int Stack[MAXN],top;
//相对于List[0]的极角排序
bool _cmp(point p1,point p2)
{
    double tmp = (p1-List[0])^(p2-List[0]);
    if(sgn(tmp) > 0)return true;
    else if(sgn(tmp) == 0 && sgn(dist(p1,List[0]) - dist(p2,List[0])) <= 0)
        return true;
    else return false;
}
void Graham(int n)
{
    point p0;
    int k = 0;
    p0 = List[0];
//找最下边的一个点
    for(int i = 1; i < n; i++)
    {
        if( (p0.y > List[i].y) || (p0.y == List[i].y && p0.x > List[i].x) )
        {
            p0 = List[i];
            k = i;
        }
    }
    swap(List[k],List[0]);
    sort(List+1,List+n,_cmp);
    if(n == 1)
    {
        top = 1;
        Stack[0] = 0;
        return;
    }
    if(n == 2)
    {
        top = 2;
        Stack[0] = 0;
        Stack[1] = 1;
        return ;
    }
    Stack[0] = 0;
    Stack[1] = 1;
    top = 2;
    for(int i = 2; i < n; i++)
    {
        while(top > 1 &&
                sgn((List[Stack[top-1]]-List[Stack[top-2]])^(List[i]-List[Stack[top-2]])) <=
                0)top--;
        Stack[top++] = i;
    }
}

3、平面最近点对(HDU 1007)

#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <iostream>
#include <math.h>
using namespace std;
const double eps = 1e-6;
const int MAXN = 100010;
const double INF = 1e20;
struct Point
{
    double x,y;
};
double dist(Point a,Point b)
{
    return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y));
}
Point p[MAXN];
Point tmpt[MAXN];
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;
}
double Closest_Pair(int left,int right)
{
    double d = INF;
    if(left == right)return d;
    if(left + 1 == right)
        return dist(p[left],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(p[mid].x - p[i].x) <= d)
            tmpt[k++] = p[i];
    }
    sort(tmpt,tmpt+k,cmpy);
    for(int i = 0; i <k; i++)
    {
        for(int j = i+1; j < k && tmpt[j].y - tmpt[i].y < d; j++)
        {
            d = min(d,dist(tmpt[i],tmpt[j]));
        }
    }
    return d;
}
int main()
{
    int n;
    while(scanf("%d",&n)==1 && n)
    {
        for(int i = 0; i < n; i++)
            scanf("%lf%lf",&p[i].x,&p[i].y);
        sort(p,p+n,cmpxy);
        printf("%.2lf\n",Closest_Pair(0,n-1)/2);
    }
    return 0;
}

4、旋转卡壳

  4.1 求解平面最远点对(POJ 2187 Beauty Contest)

struct Point
{
    int x,y;
    Point(int _x = 0,int _y = 0)
    {
        x = _x;
        y = _y;
    }
    Point operator -(const Point &b)const
    {
        return Point(x - b.x, y - b.y);
    }
    int operator ^(const Point &b)const
    {
        return x*b.y - y*b.x;
    }
    int operator *(const Point &b)const
    {
        return x*b.x + y*b.y;
    }
    void input()
    {
        scanf("%d%d",&x,&y);
    }
};
//距离的平方
int dist2(Point a,Point b)
{
    return (a-b)*(a-b);
}
//******二维凸包,int***********
const int MAXN = 50010;
Point list[MAXN];
int Stack[MAXN],top;
bool _cmp(Point p1,Point p2)
{
    int tmp = (p1-list[0])^(p2-list[0]);
    if(tmp > 0)return true;
    else if(tmp == 0 && dist2(p1,list[0]) <= dist2(p2,list[0]))
        return true;
    else return false;
}
void Graham(int n)
{
    Point p0;
    int k = 0;
    p0 = list[0];
    for(int i = 1; i < n; i++)
        if(p0.y > list[i].y || (p0.y == list[i].y && p0.x > list[i].x))
        {
            p0 = list[i];
            k = i;
        }
    swap(list[k],list[0]);
    sort(list+1,list+n,_cmp);
    if(n == 1)
    {
        top = 1;
        Stack[0] = 0;
        return;
    }
    if(n == 2)
    {
        top = 2;
        Stack[0] = 0;
        Stack[1] = 1;
        return;
    }
    Stack[0] = 0;
    Stack[1] = 1;
    top = 2;
    for(int i = 2; i < n; i++)
    {
        while(top > 1 &&
                ((list[Stack[top-1]]-list[Stack[top-2]])^(list[i]-list[Stack[top-2]])) <= 0)
            top--;
        Stack[top++] = i;
    }
}
//旋转卡壳,求两点间距离平方的最大值
int rotating_calipers(Point p[],int n)
{
    int ans = 0;
    Point v;
    int cur = 1;
    for(int i = 0; i < n; i++)
    {
        v = p[i]-p[(i+1)%n];
        while((v^(p[(cur+1)%n]-p[cur])) < 0)
            cur = (cur+1)%n;
        ans = max(ans,max(dist2(p[i],p[cur]),dist2(p[(i+1)%n],p[(cur+1)%n])));
    }
    return ans;
}
Point p[MAXN];
int main()
{
    int n;
    while(scanf("%d",&n) == 1)
    {
        for(int i = 0; i < n; i++)list[i].input();
        Graham(n);
        for(int i = 0; i < top; i++)p[i] = list[Stack[i]];
        printf("%d\n",rotating_calipers(p,top));
    }
    return 0;
}

4.2 求解平面点集最大三角形

//旋转卡壳计算平面点集最大三角形面积
int rotating_calipers(Point p[],int n)
{
    int ans = 0;
    Point v;
    for(int i = 0; i < n; i++)
    {
        int j = (i+1)%n;
        int k = (j+1)%n;
        while(j != i && k != i)
        {
            ans = max(ans,abs((p[j]-p[i])^(p[k]-p[i])));
            while( ((p[i]-p[j])^(p[(k+1)%n]-p[k])) < 0 )
                k = (k+1)%n;
            j = (j+1)%n;
        }
    }
    return ans;
}
Point p[MAXN];
int main()
{
    int n;
    while(scanf("%d",&n) == 1)
    {
        if(n == -1)break;
        for(int i = 0; i < n; i++)list[i].input();
        Graham(n);
        for(int i = 0; i < top; i++)p[i] = list[Stack[i]];
        printf("%.2f\n",(double)rotating_calipers(p,top)/2);
    }
    return 0;
}

4.3 求解两凸包最小距离(POJ 3608)

const double eps = 1e-8;
int sgn(double x)
{
    if(fabs(x) < eps)return 0;
    if(x < 0)return -1;
    else return 1;
}
struct Point
{
    double x,y;
    Point(double _x = 0,double _y = 0)
    {
        x = _x;
        y = _y;
    }
    Point operator -(const Point &b)const
    {
        return Point(x - b.x, y - b.y);
    }
    double operator ^(const Point &b)const
    {
        return x*b.y - y*b.x;
    }
    double operator *(const Point &b)const
    {
        return x*b.x + y*b.y;
    }
    void input()
    {
        scanf("%lf%lf",&x,&y);
    }
};
struct Line
{
    Point s,e;
    Line() {}
    Line(Point _s,Point _e)
    {
        s = _s;
        e = _e;
    }
};
//两点间距离
double dist(Point a,Point b)
{
    return sqrt((a-b)*(a-b));
}
//点到线段的距离,返回点到线段最近的点
Point NearestPointToLineSeg(Point P,Line L)
{
    Point result;
    double t = ((P-L.s)*(L.e-L.s))/((L.e-L.s)*(L.e-L.s));
    if(t >=0 && t <= 1)
    {
        result.x = L.s.x + (L.e.x - L.s.x)*t;
        result.y = L.s.y + (L.e.y - L.s.y)*t;
    }
    else
    {
        if(dist(P,L.s) < dist(P,L.e))
            result = L.s;
        else result = L.e;
    }
    return result;
}
/*
* 求凸包,Graham算法
* 点的编号0~n-1
* 返回凸包结果Stack[0~top-1]为凸包的编号
*/const int MAXN = 10010;
Point list[MAXN];
int Stack[MAXN],top;
//相对于list[0]的极角排序
bool _cmp(Point p1,Point p2)
{
    double tmp = (p1-list[0])^(p2-list[0]);
    if(sgn(tmp) > 0)return true;
    else if(sgn(tmp) == 0 && sgn(dist(p1,list[0]) - dist(p2,list[0])) <= 0)
        return true;
    else return false;
}
void Graham(int n)
{
    Point p0;
    int k = 0;
    p0 = list[0];
//找最下边的一个点
    for(int i = 1; i < n; i++)
    {
        if( (p0.y > list[i].y) || (p0.y == list[i].y && p0.x > list[i].x) )
        {
            p0 = list[i];
            k = i;
        }
    }
    swap(list[k],list[0]);
    sort(list+1,list+n,_cmp);
    if(n == 1)
    {
        top = 1;
        Stack[0] = 0;
        return;
    }
    if(n == 2)
    {
        top = 2;
        Stack[0] = 0;
        Stack[1] = 1;
        return ;
    }
    Stack[0] = 0;
    Stack[1] = 1;
    top = 2;
    for(int i = 2; i < n; i++)
    {
        while(top > 1 &&
                sgn((list[Stack[top-1]]-list[Stack[top-2]])^(list[i]-list[Stack[top-2]])) <=
                0)
            top--;
        Stack[top++] = i;
    }
}
//点p0到线段p1p2的距离
double pointtoseg(Point p0,Point p1,Point p2)
{
    return dist(p0,NearestPointToLineSeg(p0,Line(p1,p2)));
}//平行线段p0p1和p2p3的距离
double dispallseg(Point p0,Point p1,Point p2,Point p3)
{
    double ans1 = min(pointtoseg(p0,p2,p3),pointtoseg(p1,p2,p3));
    double ans2 = min(pointtoseg(p2,p0,p1),pointtoseg(p3,p0,p1));
    return min(ans1,ans2);
}
//得到向量a1a2和b1b2的位置关系
double Get_angle(Point a1,Point a2,Point b1,Point b2)
{
    return (a2-a1)^(b1-b2);
}
double rotating_calipers(Point p[],int np,Point q[],int nq)
{
    int sp = 0, sq = 0;
    for(int i = 0; i < np; i++)
        if(sgn(p[i].y - p[sp].y) < 0)
            sp = i;
    for(int i = 0; i < nq; i++)
        if(sgn(q[i].y - q[sq].y) > 0)
            sq = i;
    double tmp;
    double ans = dist(p[sp],q[sq]);
    for(int i = 0; i < np; i++)
    {
        while(sgn(tmp = Get_angle(p[sp],p[(sp+1)%np],q[sq],q[(sq+1)%nq])) < 0)
            sq = (sq+1)%nq;
        if(sgn(tmp) == 0)
            ans = min(ans,dispallseg(p[sp],p[(sp+1)%np],q[sq],q[(sq+1)%nq]));
        else ans = min(ans,pointtoseg(q[sq],p[sp],p[(sp+1)%np]));
        sp = (sp+1)%np;
    }
    return ans;
}
double solve(Point p[],int n,Point q[],int m)
{
    return min(rotating_calipers(p,n,q,m),rotating_calipers(q,m,p,n));
}
Point p[MAXN],q[MAXN];
int main()
{
    int n,m;
    while(scanf("%d%d",&n,&m) == 2)
    {
        if(n == 0 && m == 0)break;
        for(int i = 0; i < n; i++)
            list[i].input();
        Graham(n);
        for(int i = 0; i < top; i++)
            p[i] = list[i];
        n = top;
        for(int i = 0; i < m; i++)
            list[i].input();
        Graham(m);
        for(int i = 0; i < top; i++)
            q[i] = list[i];
        m = top;
        printf("%.4f\n",solve(p,n,q,m));
    }
    return 0;
}

5、半平面交  

  5.1 半平面交模板(from UESTC)  

const double eps = 1e-8;
const double PI = acos(-1.0);
int sgn(double x)
{
    if(fabs(x) < eps) return 0;
    if(x < 0) return -1;
    else return 1;
}
struct point
{
    double x,y;
    point() {}
    point(double _x,double _y)
    {
        x = _x;
        y = _y;
    }
    point operator -(const point &b)const
    {
        return point(x - b.x, y - b.y);
    }
    double operator ^(const point &b)const
    {
        return x*b.y - y*b.x;
    }
    double operator *(const point &b)const
    {
        return x*b.x + y*b.y;
    }
};
struct Line
{
    point s,e;
    double k;
    Line() {}
    Line(point _s,point _e)
    {
        s = _s;
        e = _e;
        k = atan2(e.y - s.y,e.x - s.x);
    }
    point operator &(const Line &b)const
    {
        point res = s;
        double t = ((s - b.s)^(b.s - b.e))/((s - e)^(b.s - b.e));
        res.x += (e.x - s.x)*t;
        res.y += (e.y - s.y)*t;
        return res;
    }
};
//半平面交,直线的左边代表有效区域
//这个好像和给出点的顺序有关
bool HPIcmp(Line a,Line b)
{
    if(fabs(a.k - b.k) > eps)return a.k < b.k;
    return ((a.s - b.s)^(b.e - b.s)) < 0;
}
Line Q[110];
//第一个位代表半平面交的直线,第二个参数代表直线条数,第三个参数是相交以后把
//所得点压栈,第四个参数是栈有多少个元素
void HPI(Line line[], int n, point res[], int &resn)
{
    int tot = n;
    sort(line,line+n,HPIcmp);
    tot = 1;
    for(int i = 1; i < n; i++)
        if(fabs(line[i].k - line[i-1].k) > eps)
            line[tot++] = line[i];
    int head = 0, tail = 1;
    Q[0] = line[0];
    Q[1] = line[1];
    resn = 0;
    for(int i = 2; i < tot; i++)
    {
        if(fabs((Q[tail].e-Q[tail].s)^(Q[tail-1].e-Q[tail-1].s)) < eps ||
                fabs((Q[head].e-Q[head].s)^(Q[head+1].e-Q[head+1].s)) < eps)
            return;
        while(head < tail && (((Q[tail]&Q[tail-1]) -
                               line[i].s)^(line[i].e-line[i].s)) > eps)
            tail--;
        while(head < tail && (((Q[head]&Q[head+1]) -
                               line[i].s)^(line[i].e-line[i].s)) > eps)
            head++;
        Q[++tail] = line[i];
    }
    while(head < tail && (((Q[tail]&Q[tail-1]) -
                           Q[head].s)^(Q[head].e-Q[head].s)) > eps)
        tail--;
    while(head < tail && (((Q[head]&Q[head-1]) -
                           Q[tail].s)^(Q[tail].e-Q[tail].e)) > eps)
        head++;
    if(tail <= head + 1)return;
    for(int i = head; i < tail; i++)
        res[resn++] = Q[i]&Q[i+1];
    if(head < tail - 1)
        res[resn++] = Q[head]&Q[tail];
}

5.2 普通半平面交写法

POJ 1750
const double eps = 1e-18;
int sgn(double x)
{
    if(fabs(x) < eps)return 0;
    if(x < 0)return -1;
    else return 1;
}
struct Point
{
    double x,y;
    Point() {} Point(double _x,double _y)
    {
        x = _x;
        y = _y;
    }
    Point operator -(const Point &b)const
    {
        return Point(x - b.x, y - b.y);
    }
    double operator ^(const Point &b)const
    {
        return x*b.y - y*b.x;
    }
    double operator *(const Point &b)const
    {
        return x*b.x + y*b.y;
    }
};
//计算多边形面积
double CalcArea(Point p[],int n)
{
    double res = 0;
    for(int i = 0; i < n; i++)
        res += (p[i]^p[(i+1)%n]);
    return fabs(res/2);
}
//通过两点,确定直线方程
void Get_equation(Point p1,Point p2,double &a,double &b,double &c)
{
    a = p2.y - p1.y;
    b = p1.x - p2.x;
    c = p2.x*p1.y - p1.x*p2.y;
}
//求交点
Point Intersection(Point p1,Point p2,double a,double b,double c)
{
    double u = fabs(a*p1.x + b*p1.y + c);
    double v = fabs(a*p2.x + b*p2.y + c);
    Point t;
    t.x = (p1.x*v + p2.x*u)/(u+v);
    t.y = (p1.y*v + p2.y*u)/(u+v);
    return t;
}
Point tp[110];
void Cut(double a,double b,double c,Point p[],int &cnt)
{
    int tmp = 0;
    for(int i = 1; i <= cnt; i++)
    {
//当前点在左侧,逆时针的点
        if(a*p[i].x + b*p[i].y + c < eps)tp[++tmp] = p[i];
        else
        {
            if(a*p[i-1].x + b*p[i-1].y + c < -eps)
                tp[++tmp] = Intersection(p[i-1],p[i],a,b,c);
            if(a*p[i+1].x + b*p[i+1].y + c < -eps)
                tp[++tmp] = Intersection(p[i],p[i+1],a,b,c);
        }
    }
    for(int i = 1; i <= tmp; i++)
        p[i] = tp[i];
    p[0] = p[tmp];
    p[tmp+1] = p[1];
    cnt = tmp;
}
double V[110],U[110],W[110];
int n;
const double INF = 100000000000.0;
Point p[110];
bool solve(int id)
{
    p[1] = Point(0,0);
    p[2] = Point(INF,0);
    p[3] = Point(INF,INF);
    p[4] = Point(0,INF);
    p[0] = p[4];
    p[5] = p[1];
    int cnt = 4;
    for(int i = 0; i < n; i++)
        if(i != id)
        {
            double a = (V[i] - V[id])/(V[i]*V[id]);
            double b = (U[i] - U[id])/(U[i]*U[id]);
            double c = (W[i] - W[id])/(W[i]*W[id]);
            if(sgn(a) == 0 && sgn(b) == 0)
            {
                if(sgn(c) >= 0)return false;
                else continue;
            }
            Cut(a,b,c,p,cnt);
        }
    if(sgn(CalcArea(p,cnt)) == 0)return false;
    else return true;
}
int main()
{
    while(scanf("%d",&n) == 1)
    {
        for(int i = 0; i < n; i++)
            scanf("%lf%lf%lf",&V[i],&U[i],&W[i]);
        for(int i = 0; i < n; i++)
        {
            if(solve(i))printf("Yes\n");
            else printf("No\n");
        }
    }
    return 0;
}

6、三点求圆心坐标(三角形外心)

//过三点求圆心坐标
Point waixin(Point a,Point b,Point c)
{
    double a1 = b.x - a.x, b1 = b.y - a.y, c1 = (a1*a1 + b1*b1)/2;
    double a2 = c.x - a.x, b2 = c.y - a.y, c2 = (a2*a2 + b2*b2)/2;
    double d = a1*b2 - a2*b1;
    return Point(a.x + (c1*b2 - c2*b1)/d, a.y + (a1*c2 -a2*c1)/d);
}

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;
}

 

7、求两圆相交的面积

//两个圆的公共部分面积
double Area_of_overlap(Point c1,double r1,Point c2,double r2)
{
    double d = dist(c1,c2);
    if(r1 + r2 < d + eps)return 0;
    if(d < fabs(r1 - r2) + eps)
    {
        double r = min(r1,r2);
        return PI*r*r;
    }
    double x = (d*d + r1*r1 - r2*r2)/(2*d);
    double t1 = acos(x / r1);
    double t2 = acos((d - x)/r2);
    return r1*r1*t1 + r2*r2*t2 - d*r1*sin(t1);
}

8、Pick 公式
  顶点坐标均是整点的简单多边形:面积=内部格点数目+边上格点数目/2-1

 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;
}

 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));
}

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;
}

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;
}

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;
}

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;
}

 

计算几何 综合模板

标签:ret   int()   cross   print   i++   end   通过   ***   ios   

原文地址:http://www.cnblogs.com/lalalatianlalu/p/7696549.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!