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

2016-04-11 19:32:06 Tabris_ 阅读数：1617

博客爬取于2020-06-14 22:39:20
以下为正文

版权声明：本文为Tabris原创文章，未经博主允许不得私自转载。
https://blog.csdn.net/qq_33184171/article/details/51124611

相当全的计算几何模板
http://paste.ubuntu.com/24458430/

精度控制

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

向量运算

二维空间

向量旋转矩阵

我们想将向量x,y→\overrightarrow {x, y} 以xx为轴点逆时针旋转，且旋转角为 $\alpha，， \left[\begin{array}{cc}\ x’ \ y’ \end{array}\right]=\left[\begin{array}{cc}\ \cos(\alpha) &&-\sin(\alpha) \ \sin(\alpha) && \cos(\alpha) \end{array}\right] \times \left[\begin{array}{cc}\ x \ y \end{array}\right]$

叉积
1
2
3
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));
}
两点求距离
1
2
3
double dis(Point p1,Point p2){
    return(sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)));
}

三维空间

点,线,面
1
2
3
4
5
6
7
8
9
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
1
2
3
4
5
6
7
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
1
2
3
double dmult(point u,point v){
    return u.x*v.x+u.y*v.y+u.z*v.z;
}
矢量差 U - V
1
2
3
4
5
6
7
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;
}
取平面法向量
1
2
3
4
5
6
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));
}
两点距离,单参数取向量大小
1
2
3
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));
}
向量大小
1
2
3
double vlen(point p){
    return sqrt(p.x*p.x+p.y*p.y+p.z*p.z);
}
判断四点是不是共面
1
2
3
4
bool judge(point a,point b,point c,point d){
    double tmp = dmult(prec(a,b,c),smult(d,a));
    return  ( abs(tmp) < eps );
}
点到平面距离
1
2
3
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));
}
三角形面积
1
2
3
4
5
6
7
8
9
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;
}

极角排序

//极角排序

1
2
3
4
5
6
7
8
9
10
11
bool dy(double x,double y)  {   return x > y + eps;} // x > y
bool xy(double x,double y)  {   return x < y - eps;} // x < y
bool dyd(double x,double y) {   return x > y - eps;} // x >= y
bool xyd(double x,double y) {   return x < y + eps;} // x <= y
bool dd(double x,double y)  {   return fabs( x - y ) < eps;}  // x == y

bool 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)O(nlog_{2}{n})复杂度寻求凸包

1

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
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;
}

线与线

线段相交
1
2
3
4
5
6
7
8
9
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);
}
求两线段交点
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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;
}

两圆相交

两圆相交面积
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
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。

1
2
3
4
5
6
7
8
9
10
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;
}