wys的个人博客

你有很多事放不下？做人要潇洒一点~

0%

计算几何模板(持续更新)

Posted on Edited on In

1.计算叉积

1
2
3
4
double corss(PDD a,PDD b){
    return a.x*b.y - a.y*b.x;
}

当 $a$ 在 $b$ 的逆时针方向时，叉积为正，当 $a$ 在 $b$ 的顺时针方向时叉积为负

2.计算三角形面积

1
2
3
double area(PDD a,PDD b,PDD c){
    return corss(b-a,c-a)/2;
}

这里计算的时三角形的有向面积，当 $ab$ 在 $ac$ 逆时针方向时面积为正，顺时针方向时面积为负。

3.符号函数

1
2
3
4
5
int sign(double x){
    if(fabs(x) < eps) return 0;
    else if(x < 0) return -1;
    else return 1;
}

当 $x$ 在误差范围内返回 $0$，其余根据情况返回 $-1$ 和 $1$。

4.求凸包

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
double andrew(){
    
    sort(q+1,q+n+1);
    int top = 0;
    for(int i=1;i<=n;i++){
        while(top >= 2 && area(q[stk[top-1]],q[stk[top]],q[i]) < 0){
            if(area(q[stk[top-1]],q[stk[top]],q[i]) < 0)
                used[stk[top -- ]] = false;
            else top--;
        }
        stk[++ top] = i;
        used[i] = true;
    }
    used[1] = false;
    int k = top;
    for(int i=n;i > 0;i--){
        if(used[i]) continue;
        while (top > top && area(q[stk[top - 1]], q[stk[top]], q[i]) < 0)
            top -- ;
        stk[ ++ top] = i;
    }
    double res = 0;
    for (int i = 2; i <= top; i ++ )
        res += get_dist(q[stk[i - 1]], q[stk[i]]);
    return res;
}

当保留凸包上的点时，area(q[stk[top - 1]], q[stk[top]], q[i]) < 0 中的 < 改为 <=

5.两点之间距离

1
2
3
4
5
double get_dist(PDD a,PDD b){
    double dx = a.x - b.x;
    double dy = a.y - b.y;
    return sqrt(dx*dx + dy*dy);
}

返回两点之间的距离

6.直线旋转

1
2
3
PDD rotate(PDD a,double b){
    return {a.x*cos(b) + a.y*sin(b),-a.x*sin(b) + a.y*cos(b)};
}

向量 $a$ 顺时针旋转弧度 $b$

7.直线单位法向量

1
2
3
4
PDD normal(PDD a){
	double L = length(a);
	return {-a.y/L,a.x/L};
}

返回直线的单位法向量

8.直线的角度

1
2
3
double get_angle(const Line& a){
    return atan2(a.ed.y - a.st.y,a.ed.x - a.st.x);
}

返回直线的角度，要注意的时 $atan2$ 的前面为 $y$ 坐标，后面为 $x$ 坐标。

9.两直线的交点

1
2
3
4
5
6
7
8
9
PDD get_line_intersection(PDD p,PDD v,PDD q,PDD w){
    auto u = p-q;
    double t = cross(w,u) / cross(v,w);
    return {p.x+v.x*t,p.y + v.y*t};
}

PDD get_line_intersection(Line a,Line b){
    return get_line_intersection(a.st,a.ed - a.st,b.st,b.ed-b.st);
}

返回两直线的交点

注意在精度要求比较高的时候可以自定义以下分数类来确保精度

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
typedef pair<Point,Point> PPP;

struct Point
{
    int x,y;
    Point(int a, int b){
        x = a/__gcd(a,b);
        y = b/__gcd(a,b);
    }
    bool operator < (const Point &o) const{
        return x*o.y < y*o.x;
    }
};

PPP get_line_intersection(PDD p,PDD v,PDD q,PDD w){
    PDD u = p-q;
    int a = cross(u,w),b = cross(w,v);
    Fraction t = (cross(u,w),cross(w,v));
    return {Point(a*v.x + p.x*b,b),Point(a*v.y + p.y*b,b)};
}

10.半平面交

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
bool cmp(const Line&a,const Line&b){
    double A = get_angle(a),B = get_angle(b);
    if(!dcmp(A,B)) return area(a.st,a.ed,b.ed) < 0;
    return A < B;
}

bool on_right(Line& a,Line b,Line& c){
    auto o = get_line_intersection(b,c);
    return sign(area(a.st,a.ed,o)) <= 0;
}

double half_plane_intersection(){
    sort(line,line+cnt,cmp);
    int hh = 0,tt = -1;
    for(int i=0;i<cnt;i++){
        if(i && !dcmp(get_angle(line[i]),get_angle(line[i-1]))) continue;
        while(hh + 1 <= tt && on_right(line[i],line[q[tt-1]],line[q[tt]])) tt--;
        while(hh + 1 <= tt && on_right(line[i],line[q[hh]],line[q[hh+1]])) hh++;
        q[++ tt] = i;
    }
    while (hh + 1 <= tt && on_right(line[q[hh]], line[q[tt - 1]], line[q[tt]])) tt -- ;
    while (hh + 1 <= tt && on_right(line[q[tt]], line[q[hh]], line[q[hh + 1]])) hh ++ ;
    q[ ++ tt] = q[hh];
    int k = 0;
    for(int i=hh;i<tt;i++){
        ans[k++] = get_line_intersection(line[q[i]],line[q[i+1]]);
    }
    double res  = 0;
    for(int i=1;i+1<k;i++){
        res += area(ans[0],ans[i],ans[i+1]);
    }
    return res;
}

第一个函数用来对直线进行排序，当直线的角度相同时，保留左侧的直线。

第二个函数用来判断 $bc$ 的交点是否在直线的右侧。

第三个函数用来得到半平面交上的点，并返回半平面交的面积。

11.旋转卡壳

1
2
3
4
5
6
7
8
9
10
int rotating_calipers(){
    if(top <= 2) return get_dist(q[0],q[n-1]);
    int res = 0;
    for(int i=0,j=2;i<top;i++){
        auto d = q[stk[i]], e = q[stk[i+1]];
        while(area(d,e,q[stk[j]]) < area(d,e,q[stk[j+1]])) j = (j+1)%top;
        res = max(res,max(get_dist(d,q[stk[j]]),get_dist(e,q[stk[j]])));
    }
    return res;
}

对于旋转卡壳应该先求凸包，第一句话是为了特判当所有点在一条直线上的情况。

12.判断两直线规范相交

1
2
3
4
5
6
7
8
9
bool on_segment(PDD a1,PDD a2,PDD b1,PDD b2){
    double c1 = cross(a2-a1,b1-a1),c2 = cross(a2-a1,b2-a1);
    double c3 = cross(b2-b1,a1-b1),c4 = cross(b2-b1,a2-b1);
    return sign(c1) * sign(c2) < 0 && sign(c3) * sign(c4) < 0 ;
}

bool on_segment(Line a,Line b){
    return on_segment(a.st,a.ed,b.st,b.ed);
}

我们定义“规范相交”为两线段恰好有一个公共点，且不在任何一条线段的端点。线段规范相交的充要条件是：每条线段的两个端点，都在另一条线段的两侧（这里的“两侧”是指叉积的符号不同）。

以上方法我们称之为跨立实验。

我认为在确定了两线段不平行以后，可以先求出直线交点，然后用以下代码判断是否相交。

1
2
3
bool OnSegment(PDD p,PDD a1,PDD a2){
	return dcmp(Dot(a1-p,a2-p)) < 0;
}

三维几何相关

13.三维叉积

大小为 $\abs{a}\abs{b}sin(\theta)$

$ A\times B = \left|\begin{array}{cccc} i &j & k \\ x_1 &y_1 & z_1 \\ x_2 & y_2 &z_2 \\ \end{array}\right| = (y_1z_2-y_2z_1)\cdot i+(z_1x_2-z_2x_1)\cdot j+(x_1y_2-x_2y_1)\cdot k$

方向遵守右手定则 及对于 $a \times b$ 右手四指开始指向$a$，向$b$的方向弯曲，大拇指的方向即为叉积的方向。同时与平面垂直，即为法向量的方向。

1
2
3
Point cross(Point a,Point b){
    return {a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x};
}

两个三维向量的叉积

1
2
3
double area(Point a,Point b,Point c){
    return Length(cross((b-a),(c-a)))/2;
}

利用叉积计算三角形的面积

14.空间三角形相关

1
2
3
4
5
bool PointinTri(Point p,Point p0,Point p1,Point p2){
    double a1 = area(p,p0,p1),a2 = area(p,p1,p2),a3 = area(p,p0,p2);
    double s = area(p0,p1,p2);
    return !sign(a1+a2+a3 - s);
}

判断一个点是否在三角形内部

1
2
3
4
5
6
7
8
9
10
bool TriSegInersection(Point p0,Point p1,Point p2,Point a,Point b,Point& p){
    Point n = cross(p1-p0,p2-p0);
    if(!sign(n&(b-a))) return false;
    else{
        double t = (n&(p0-a))/(n&(b-a));
        if(sign(t) < 0 || sign(t-1) > 0) return false;
        p = a + (b-a)*t;
        return PointinTri(p,p0,p1,p2);
    }
}

判断一条线段是否与三角形相交，并通过引用返回交点坐标

1
2
3
4
5
6
7
8
bool TriTriIntersection(Point* t1,Point* t2){
    Point p;
    for(int i=0;i<3;i++){
        if(TriSegInersection(t1[0],t1[1],t1[2],t2[i],t2[(i+1)%3],p)) return true;
        if(TriSegInersection(t2[0],t2[1],t2[2],t1[i],t1[(i+1)%3],p)) return true;
    }
    return false;
}

判断两个三角形是否有交点

15.判断点是否在多边形内部

1
2
3
4
5
6
7
8
9
10
11
12
13
bool check_point(PDD p,PDD* poly,int n){
    int wn = 0;
    for(int i=0;i<n;i++){
    if(is_point_on_seg(p,poly[i],poly[(i+1)%n])) return 1;
    int k = sign(cross(poly[(i+1)%n]-poly[i],p-poly[i]));
    int d1 = sign(poly[i].y - p.y);
    int d2 = sign(poly[(i+1)%n].y - p.y);
    if(k > 0 && d1 <= 0 && d2 > 0) wn++;
    if(k < 0 && d2 <= 0 && d1 > 0) wn--;
    }
    if(wn != 0) return 1;
    return 0;
}

用射线法判断点是否在多边形内部