計算幾何筆記

 

本文中全部圖片及其餘引用均已得到原做者贊成

向量

常規操做

struct Vector {
    double x, y;
    Vector(double x = 0, double y = 0) : x(x), y(y) {};
    Vector operator + (const Vector &A) const {
        return Vector(x + A.x, y + A.y);
    }//向量的加法 
    Vector operator - (const Vector &A) const {
        return Vector(x - A.x, y - A.y);
    }//減法 
    Vector operator * (const double &A) const {
        return Vector(x * A, y * A);
    }//數乘 
    Vector operator / (const double &A) const {
        return Vector(x / A, y / A);
    }//數除 
    bool operator == (const Vector &A) const  {
        return dcmp(x - A.x) == 0 && dcmp(y - A.y) == 0;
    }//相等 
};

向量的極角

極角：這裏指從$x$軸旋轉到向量方向的弧度

double PolarAngle(Vector A) {
    return atan2(A.y, A.x);
}//向量的極角 

向量的旋轉

若向量$(x, y)$旋轉角度爲$a$，則旋轉後的向量爲$(xcosa - ysina, y cosa + xsina)$

證實：

設旋轉以前的向量的極角爲$t$，半徑爲$r$

那麼

$$x = rcost$$

$$y = rsint$$

旋轉時候的向量爲

$$x' = rcos(t + a) = r(costcosa - sintsina) = xcosa - ysina$$

$$y' = rsin(t + a) = r(sintcosa + costsina) = ycosa + xsina$$

Vector rotate(Vector &A, double a) {
    return A = Vector(A.x * cos(a) - A.y * sin(a), A.x * sin(a) + A.y * cos(a));
}//向量旋轉，a爲弧度 

向量的點積

$a \cdot b = |a||b|cos<a,b>$

兩向量的點積獲得的是標量，即一個向量的模長乘另外一個向量在該向量上正投影的數量

double Dot(Vector A, Vector B) {
    return A.x * B.x +  A.y * B.y;
}//兩向量點積 

向量的叉積

$a \times b = |a||b| sin<a,b>$

兩向量叉積獲得的是向量，在二維平面中獲得的是三維空間中與這兩個向量垂直的向量

在平面中，向量$v$和$w$的叉積等於$v$和$w$組成的三角形的有向面積的兩倍

記$cross(v,w)$表示兩向量的叉積，若$cross(v,w) > 0 $則說明$w$在$v$的左側，不然$w$在$v$的右側

double Cross(Vector A, Vector B) {
    return A.x * B.y - A.y * B.x; 
}//兩向量叉積 

計算三角形面積

直接利用叉積的定義

double Area(Point A, Point B, Point C) {
    return fabs(Cross(B - A, C - A) / 2);
}//計算三角形的面積 

計算向量的長度

直接利用點積的定義

double Length(Vector A) {
    return sqrt(Dot(A, A));
}//計算向量的長度 

計算向量的夾角

一樣直接利用點積的定義

double Angle(Vector A, Vector B) {
    return acos(Dot(A, B) / Length(A) / Length(B));
}//計算向量的夾角 

線段

判斷兩直線的相對位置

判斷$P_1P_0$與$P_1P_2$的相對位置關係，能夠轉化爲判斷$P_1P_0$與$P_2P_0$叉積的正負性

int Direction(Vector P1, Vector P2, Vector P0) {
    int opt = Cross(P1 - P0, P2 - P0);
    return dcmp(opt);
}//判斷P1-P0,P1-P2的相對位置關係，-1爲逆時針，1爲順時針(P1P0順時針旋轉到P1P2)，0爲共線 

判斷兩直線的交點

尼瑪看不懂

Point GetLineIntersection(Point P, Vector V, Point Q, Vector W) {
    Vector u = P - Q;
    double t = Cross(W, u) / Cross(V, W);
    return P + V * t;
}//判斷兩直線(P + tv,Q + tW)的交點(看不懂直接上y = kx + b吧) 

計算點到直線的距離

利用叉積算出他們圍城的平行四邊形的面積，再除以底，高即爲距離

double DistanceToLine(Point P, Point A, Point B) {
    Vector v1 = B - A, v2 = P - A;
    return fabs(Cross(v1, v2)) / Length(v1);
}//計算點P到直線AB的距離（平行四邊形面積 / 底）

計算點到線段的距離

這個要分三種狀況討論

double DistanceToSegment(Point P, Point A, Point B) {
    if(A == B) return Length(P - A);
    if(dcmp(Dot(P - A, B - A)) < 0) return Length(P - A);
    if(dcmp(Dot(P - B, A - B)) < 0) return Length(P - B);
    return DistanceToLine(P, A, B);
}//計算點P到線段AB的距離 

 

多邊形

計算外接圓的圓心

某次考試遇到的

    double x1, x2, x3, y1, y2, y3;
    x1  =  p[a].x;
    x2  =  p[b].x;
    x3  =  p[c].x;
    y1  =  p[a].y;
    y2  =  p[b].y;
    y3  =  p[c].y;
    //求外接圓圓心
    double t1 = x1 * x1 + y1 * y1;
    double t2 = x2 * x2 + y2 * y2;
    double t3 = x3 * x3 + y3 * y3;
    double temp = x1 * y2 + x2 * y3 + x3 * y1 - x1 * y3 - x2 * y1 - x3 * y2;
    double x = (t2 * y3 + t1 * y2 + t3 * y1 - t2 * y1 - t3 * y2 - t1 * y3) / temp / 2;
    double y = (t3 * x2 + t2 * x1 + t1 * x3 - t1 * x2 - t2 * x3 - t3 * x1) / temp / 2;
    Point ma;//中點 
    ma.x = x;
    ma.y = y;

 

計算多邊形的有向面積

將$n$邊形拆成三角形

double PolygonAread(Point *P, int N) {
    double area = 0;
    for(int i = 1; i <= N - 1; i++) 
        area += Cross(P[i] - P[0], P[i + 1] - P[0]);
    return area / 2;
}//計算多邊形的有向面積

判斷點是否在多邊形內部

基本思想：從點$P$向右作一條射線，判斷從無限遠處到點$P$，射線穿過了幾條邊

有兩種須要特判的狀況

1.射線與某條邊重合，該邊不統計入答案

2.射線與端點重合

此時，咱們欽定邊是由編號小的連向編號大的

若邊從上到下穿過了射線，包含終點不包含起點

若邊從下到上穿過了射線，包含起點不包含重點

int isPointInPolygon(Point P, Point Poly[], int N) {
    int wn = 0, k, d1, d2;
    for(int i = 0; i < N; i++) {
        if(!dcmp(DistanceToSegment(P, Poly[i], Poly[(i + 1) % N]))) return -1;
    //點在線段的邊界上
        k = dcmp(Cross(Poly[(i + 1) % N] - Poly[i], P - Poly[i]));
        d1 = dcmp(Poly[i].y - P.y);
        d2 = dcmp(Poly[(i + 1) % N].y - P.y);
        if(k > 0 && d1 <= 0 && d2 > 0) wn++;//點在左，下上穿
        if(k > 0 && d2 <= 0 && d1 > 0) wn++;//點在右，上下穿
        return wn & 1; // 1：內 2：外    
    }    
}//判斷點是否在多邊形內部

 

對踵點

定義：若點對$(a, b)$均爲多邊形上的點且存在過$a$點的切線與過$b$點的切線平行，則成$(a, b)$爲多邊形上的對踵點

計算方法：

設$p_{ymin}$表示$y$最小的點，$q_{ymax}$表示$y$最大的點，顯然它們是一對對踵點

接下來以相同的角速度逆時針旋轉兩條射線，當其中的一條穿過多邊形的下一個端點$p_{next}$時，用$p_{next}$做爲新的端點，同時與$q_{pre}$構成新的對踵點。

在這個算法中，咱們快速的找出兩條射線到底是哪條先穿過下一個端點，咱們能夠用叉積來優化這一過程。

 

求凸多邊形的直徑

定義：凸多邊形的直徑爲多邊形的上最遠的點對的距離

很顯然，直徑必定是在對踵點處取得，直接枚舉對踵點便可

double RotatingCaliper_diameter(Point Poly[], int N) {
    int p = 0, q = N - 1, fl;
    double ret = 0;
    for(int i = 0; i < N; i++) {
        if(Poly[i].y > Poly[q].y) q = i;
        if(Poly[i].y < Poly[p].y) p = i;
    }
    for(int i = 0; i < N * 3; i++) {
        ret = max(ret, Length(Poly[p % N] - Poly[q % N]));
        fl = dcmp(Cross(Poly[(p + 1) % N] - Poly[p % N], Poly[q % N] - Poly[(q + 1) % N]));
        if(!fl) {
            ret = max(ret, Length(Poly[p % N] - Poly[(q + 1) % N]));
            ret = max(ret, Length(Poly[q % N] - Poly[(p + 1) % N]));
            p++, q++;
        } else {
            if(fl > 0) p++;
            else q++;
        }
    }
}//計算多邊形直徑 

 

凸多邊形的寬度

凸多邊形最小面積外接矩形

凸包-Andrew算法

首先按照$x$爲第一關鍵字，$y$爲第二關鍵字從小到大排序，並刪除重複的點

用棧維護凸包內的點

一、把$p_1, p_2$放入棧中

二、若$p_{i{(i > 3)}}$在直線$p_{i - 1}, p_{i - 2}$的右側，則不斷的彈出棧頂，直到該點在直線左側

三、此時咱們已經獲得了下凸包，那麼反過來從$p_n$再作一次便可獲得下凸包

題目連接

// luogu-judger-enable-o2
#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
const int eps = 1e-10;
int dcmp(double x) {
    if(fabs(x) < eps) return 0;
    return x < 0 ? -1 : 1;
}
#define Point Vector 
struct Vector {
    double x, y;
    Vector(double x = 0, double y = 0) : x(x), y(y) {};
    bool operator < (const Vector &rhs) const {
        return dcmp(x - rhs.x) == 0 ? y < rhs.y : x < rhs.x;
    }
    Vector operator - (const Vector &rhs) const {
        return Vector(x - rhs.x, y - rhs.y);
    }
};
double Cross(Vector A, Vector B) {
    return A.x * B.y - A.y * B.x; 
}
double dis(Point a, Point b) {
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
int N; 
Point p[10001], q[10001];
int top;
void Push(Point p) {
    while(Cross(q[top] - q[top - 1], p - q[top - 1]) < 0) top--;
    q[++top] = p;
}
void Andrew() {
    q[0] = q[top = 1] = p[1];
    for(int i = 2; i <= N; i++) Push(p[i]);
    for(int i = N - 1; i; --i)  Push(p[i]);
}
int main() {    
    scanf("%d", &N);
    for(int i = 1; i <= N; i++) scanf("%lf%lf", &p[i].x, &p[i].y);
    sort(p + 1, p + N + 1);
    Andrew();
    double ans = 0;
    for(int i = 1; i < top; i++)
        ans += dis(q[i], q[i + 1]);
    printf("%.2lf", ans);
    return 0;
}

凸包

 

 

 

 

 

要作的題

POJ 2187 凸多邊形直徑

 

 

 

#include<cstdio>
#include<cmath>
using namespace std;
const int eps = 1e-10;
int dcmp(double x) {
    if(fabs(x) < eps) return 0;
    return x < 0 ? -1 : 1;
}
#define Point Vector 
struct Vector {
    double x, y;
    Vector(double x = 0, double y = 0) : x(x), y(y) {};
    Vector operator + (const Vector &A) const {
        return Vector(x + A.x, y + A.y);
    }//向量的加法 
    Vector operator - (const Vector &A) const {
        return Vector(x - A.x, y - A.y);
    }//減法 
    Vector operator * (const double &A) const {
        return Vector(x * A, y * A);
    }//數乘 
    Vector operator / (const double &A) const {
        return Vector(x / A, y / A);
    }//數除 
    bool operator == (const Vector &A) const  {
        return dcmp(x - A.x) == 0 && dcmp(y - A.y) == 0;
    }//相等 
};
double PolarAngle(Vector A) {
    return atan2(A.y, A.x);
}//向量的極角 
Vector rotate(Vector &A, double a) {
    return A = Vector(A.x * cos(a) - A.y * sin(a), A.x * sin(a) + A.y * cos(a));
}//向量旋轉，a爲弧度 
double Dot(Vector A, Vector B) {
    return A.x * B.x +  A.y * B.y;
}//兩向量點積 
double Cross(Vector A, Vector B) {
    return A.x * B.y - A.y * B.x; 
}//兩向量叉積 
double Area(Point A, Point B, Point C) {
    return fabs(Cross(B - A, C - A) / 2);
}//計算三角形的面積 
double Length(Vector A) {
    return sqrt(Dot(A, A));
}//計算向量的長度 
double Angle(Vector A, Vector B) {
    return acos(Dot(A, B) / Length(A) / Length(B));
}//計算向量的夾角 
int Direction(Vector P1, Vector P2, Vector P0) {
    int opt = Cross(P1 - P0, P2 - P0);
    return dcmp(opt);
}//判斷P1-P0,P1-P2的相對位置關係，-1爲逆時針，1爲順時針(P1P0順時針旋轉到P1P2)，0爲共線 
Point GetLineIntersection(Point P, Vector V, Point Q, Vector W) {
    Vector u = P - Q;
    double t = Cross(W, u) / Cross(V, W);
    return P + V * t;
}//判斷兩直線(P + tv,Q + tW)的交點(看不懂直接上y = kx + b吧) 
double DistanceToLine(Point P, Point A, Point B) {
    Vector v1 = B - A, v2 = P - A;
    return fabs(Cross(v1, v2)) / Length(v1);
}//計算點P到直線AB的距離（平行四邊形面積 / 底）
double DistanceToSegment(Point P, Point A, Point B) {
    if(A == B) return Length(P - A);
    if(dcmp(Dot(P - A, B - A)) < 0) return Length(P - A);
    if(dcmp(Dot(P - B, A - B)) < 0) return Length(P - B);
    return DistanceToLine(P, A, B);
}//計算點P到線段AB的距離
double PolygonAread(Point *P, int N) {
    double area = 0;
    for(int i = 1; i <= N - 1; i++) 
        area += Cross(P[i] - P[0], P[i + 1] - P[0]);
    return area / 2;
}//計算多邊形的有向面積
int isPointInPolygon(Point P, Point Poly[], int N) {
    int wn = 0, k, d1, d2;
    for(int i = 0; i < N; i++) {
        if(!dcmp(DistanceToSegment(P, Poly[i], Poly[(i + 1) % N]))) return -1;
    //點在線段的邊界上
        k = dcmp(Cross(Poly[(i + 1) % N] - Poly[i], P - Poly[i]));
        d1 = dcmp(Poly[i].y - P.y);
        d2 = dcmp(Poly[(i + 1) % N].y - P.y);
        if(k > 0 && d1 <= 0 && d2 > 0) wn++;//點在左，下上穿
        if(k > 0 && d2 <= 0 && d1 > 0) wn++;//點在右，上下穿
        return wn & 1; // 1：內 2：外    
    }    
}//判斷點是否在多邊形內部

int main() {    
        
    return 0;
}

 

 

參考資料

也許是史上最不良心的低階計算幾何講解和習題集？？