【计算几何各种小模板总结贴】[不定期更新]

精度控制

Ps：尽量不要用除法，三角函数，强制类型转换等操作，否则精度损失比较大
const double PI = acos(-1.0);
const double EPS = 1e-8;
任何double的比较运算都要用eps

向量运算

二维空间

向量旋转矩阵

我们想将向量 x,y−→− 以x为轴点逆时针旋转，且旋转角为 α，
[ x′y′]=[ cos(α)sin(α)−sin(α)cos(α)]×[ xy]

叉积

double mul(Point p1,Point p2,Point p0){    return((p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y));}

两点求距离

double dis(Point p1,Point p2){    return(sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)));}

三维空间

点,线,面

struct point{    double x,y,z;};struct line{    point a,b;};struct plane{    point a,b,c;};

计算 cross product product product product U x V

point xmult(point u,point v){    point ret;    ret.x=u.y*v.z-v.y*u.z;    ret.y=u.z*v.x-u.x*v.z;    ret.z=u.x*v.y-u.y*v.x;    return ret;}

计算 dot product product product product U . V

double dmult(point u,point v){    return u.x*v.x+u.y*v.y+u.z*v.z;}

矢量差 U - V

point subt(point u,point v){    point ret;    ret.x=u.x-v.x;    ret.y=u.y-v.y;    ret.z=u.z-v.z;    return ret;}

取平面法向量

point pvec(plane s){    return xmult(subt(s.a,s.b),subt(s.b,s.c));}point pvec(point s1,point s2,point s3){    return xmult(subt(s1,s2),subt(s2,s3));}

两点距离,单参数取向量大小

double dis(point p1,point p2){    return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)+(p1.z-p2.z)*(p1.z-p2.z));}

向量大小

double vlen(point p){    return sqrt(p.x*p.x+p.y*p.y+p.z*p.z);}

判断四点是不是共面

bool judge(point a,point b,point c,point d){    double tmp = dmult(prec(a,b,c),smult(d,a));    return  ( abs(tmp) < eps );}

点到平面距离

double ptoplane(point p,point s1,point s2,point s3){    return fabs(dmult(pvec(s1,s2,s3),subt(p,s1)))/vlen(pvec(s1,s2,s3));}

三角形面积

double Area_triangle(point a,point b,point c){    double ab=dis(a,b),bc=dis(b,c),ac=dis(a,c);    double p=(ab+bc+ac)/2;    return sqrt(p* (p-ab) * (p-bc)*(p-ac));}double Area_triangle(point b,point p1,point p2){    point a=xmult(smult(b,p1),smult(b,p2));    return sqrt(a.x*a.x+a.y*a.y+a.z*a.z)/2.0;}

极角排序

//极角排序

bool dy(double x,double y)  {   return x > y + eps;} // x > ybool xy(double x,double y)  {   return x < y - eps;} // x < ybool dyd(double x,double y) {   return x > y - eps;} // x >= ybool xyd(double x,double y) {   return x < y + eps;} // x <= ybool dd(double x,double y)  {   return fabs( x - y ) < eps;}  // x == ybool cmp(point a,point b)   // 第一次排序{    if( dd(a.y ,b.y ) )  return xy(a.x, b.x);    return xy(a.y,b.y);}

凸包

grabam 扫描法

O(nlog2n)复杂度寻求凸包

void Graham_scan(Point p[],Point ch[],int n,int &len){    int i,j,k=0,top=2;    Point tmp;    //找到最下且偏左的那个点    for(i=1; i<n; i++)        if ((p[i].y<p[k].y)||((p[i].y==p[k].y)&&(p[i].x<p[k].x)))            k=i;    //将这个点指定为PointSet[0]    tmp=p[0];    p[0]=p[k];    p[k]=tmp;    //按极角从小到大,距离偏短进行排序    for (i=1; i<n-1; i++)    {        k=i;        for (j=i+1; j<n; j++)            if( (mul(p[j],p[k],p[0])>0)                    ||((mul(p[j],p[k],p[0])==0)                       &&(dis(p[0],p[j])<dis(p[0],p[k]))) )                k=j;//k保存极角最小的那个点,或者相同距离原点最近        tmp=p[i];        p[i]=p[k];        p[k]=tmp;    }    //第三个点先入栈    ch[0]=p[0];    ch[1]=p[1];    ch[2]=p[2];    //判断与其余所有点的关系    for (i=3; i<n; i++)    {        //不满足向左转的关系,栈顶元素出栈        while(mul(p[i],ch[top],ch[top-1])>=0) top--;        //当前点与栈内所有点满足向左关系,因此入栈.        ch[++top]=p[i];    }    len=top+1;}

凸包运算（周长+面积+最小园覆盖）

#include <iostream>#include <cmath>#include <cstdio>using namespace std;/*p[]：输入的点集ch[]：输出的凸包上的点集，按照逆时针方向排列n：PointSet中的点的数目len：输出的凸包上的点的个数*/struct Point{    double x,y;};//小于0,说明向量p0p1的极角大于p0p2的极角double mul(Point p1,Point p2,Point p0){    return((p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y));}double dis(Point p1,Point p2){    return(sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)));}void Graham_scan(Point p[],Point ch[],int n,int &len){    int i,j,k=0,top=2;    Point tmp;    //找到最下且偏左的那个点    for(i=1; i<n; i++)        if ((p[i].y<p[k].y)||((p[i].y==p[k].y)&&(p[i].x<p[k].x)))            k=i;    //将这个点指定为PointSet[0]    tmp=p[0];    p[0]=p[k];    p[k]=tmp;    //按极角从小到大,距离偏短进行排序    for (i=1; i<n-1; i++)    {        k=i;        for (j=i+1; j<n; j++)            if( (mul(p[j],p[k],p[0])>0)                    ||((mul(p[j],p[k],p[0])==0)                       &&(dis(p[0],p[j])<dis(p[0],p[k]))) )                k=j;//k保存极角最小的那个点,或者相同距离原点最近        tmp=p[i];        p[i]=p[k];        p[k]=tmp;    }    //第三个点先入栈    ch[0]=p[0];    ch[1]=p[1];    ch[2]=p[2];    //判断与其余所有点的关系    for (i=3; i<n; i++)    {        //不满足向左转的关系,栈顶元素出栈        while(mul(p[i],ch[top],ch[top-1])>=0) top--;        //当前点与栈内所有点满足向左关系,因此入栈.        ch[++top]=p[i];    }    len=top+1;}const int maxN=50000;  Point p[maxN];Point ch[maxN];int n;int len;int main(){    int t,l;    //scanf("%d",&t);    while(~scanf("%d",&n))    {        double distence=0.0;        for(int i=0; i<n; i++)        {            scanf("%lf%lf",&p[i].x,&p[i].y);        }        if(n==0) break;        //这是一些值得特判的情况        //if(n==0) if(n==1) if(n==2)        if(n<=1)        {            printf("0.50\n");            continue;        }        if (n == 2)        {            printf("%.2lf\n",dis(p[0],p[1])*0.50+0.50);            continue;        }        Graham_scan(p,ch,n,len);        for(int i=0; i<len; i++) //凸包 各个顶点的顺序            cout<<ch[i].x<<" "<<ch[i].y<<endl;        //求凸包的周长        double perimeters=0.0;        for(int i=1;i<len;i++)        {            perimeters+=dis(ch[i],ch[i-1]);        }        perimeters+=dis(ch[0],ch[len-1]);        printf("%.0lf\n",perimeters);        //求凸包的面积        double area=0.0;        for(int i=2; i<len; i++)        {            area+=mul(ch[i-1],ch[i],ch[0]);        }        printf("%d\n",(int )area/100);        //求图包最小圆覆盖        double maxr = -1;        double a, b, c, r, s;        for (int i=0; i<len; i++)  //枚举凸包上的点        {            for (int j=i+1; j<len; j++)            {                for (int k=j+1; k<len; k++)                {                    a = dis(ch[i], ch[j]);                    b = dis(ch[i], ch[k]);                    c = dis(ch[j], ch[k]);                    if (a*a+b*b<c*c || a*a+c*c<b*b || b*b+c*c<a*a)                        r = max(max(a, b), c) / 2.0;//钝角时  半径为最长边的一半                    else                    {                        s = fabs(mul(ch[i], ch[j], ch[k])) / 2.0;                        r = a * b * c / (s * 4.0);//三角形外接圆公式                    }                    if (maxr < r) maxr = r;                }            }        }        printf ("%.2lf\n", maxr+0.50);        //if(t) printf("\n");  //这是控制两组输出数据之间有一个空格的    }    return 0;}

线与线

线段相交

bool IsIntersected(point s1,point e1,point s2,point e2)//两个线段相交{    return(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,e1.y))&&           (multi(s1,s2,e1)*multi(s1,e1,e2)>0)&&           (multi(s2,s1,e2)*multi(s2,e2,e1)>0);}

求两线段交点

point intersection(point &A,point &B,point &C,point &D){    /* AB与CD交点*/    point p;    double a1=A.y-B.y;    double b1=B.x-A.x;    double c1=A.x*B.y-B.x*A.y;    double a2=C.y-D.y;    double b2=D.x-C.x;    double c2=C.x*D.y-D.x*C.y;    p.x=(b1*c2-b2*c1)/(a1*b2-a2*b1);    p.y=(a2*c1-a1*c2)/(a1*b2-a2*b1);    return p;}

两圆相交

两圆相交面积

struct Round {      double x, y;      double r;  };  double solve(Round a, Round b){    double d = dis(a, b);    if (d >= a.r + b.r)        return 0;    if (d <= fabs(a.r - b.r))    {        double r = a.r < b.r ? a.r : b.r;        return PI * r * r;    }    double ang1 = acos((a.r * a.r + d * d - b.r * b.r) / 2. / a.r / d);    double ang2 = acos((b.r * b.r + d * d - a.r * a.r) / 2. / b.r / d);    double ret = ang1 * a.r * a.r + ang2 * b.r * b.r - d * a.r * sin(ang1);    return ret;}

四点共面

可以由4个点构成3个向量，3个向量共面的充要条件是向量为a，b，c，存在实数x,y,z不全为零，使得xa+yb+zc=0。转化为线性代数的3个向量线性相关的行列式为0。

double is_coplanar(point a, point b, point c, point d){    point ab, ac, ad;    ab.x=a.x-b.x, ab.y=a.y-b.y, ab.z=a.z-b.z;    ac.x=a.x-c.x, ac.y=a.y-c.y, ac.z=a.z-c.z;    ad.x=a.x-d.x, ad.y=a.y-d.y, ad.z=a.z-d.z;    double gg=ab.x*ac.y*ad.z+ab.y*ac.z*ad.x+ac.x*ad.y*ab.z;    gg-=ab.z*ac.y*ad.x+ab.x*ac.z*ad.y+ab.y*ac.x*ad.z;    return gg;}