1. 程式人生 > 實用技巧 >【筆記】計算幾何

【筆記】計算幾何

本筆記正在更新! 歡迎催更!

想要更舒適的閱讀體驗?歡迎在作者部落格觀看>>

\(\sf Part\ -999\) 感謝列表

(排名不分先後)

\(\sf Part\ 1\) 前言

聽說 rui_er 神仙做夢學了二維凸包後,我也想學學了計算幾何,可是咱沒人家的天賦,可以通過睡覺學,只能看資料學了/dk

@rui_er 快寫篇文章教教我們怎麼做夢學/se(

迴歸正題,首先說明一下,本人是剛學 \(\mathsf{OI}\)

的萌新,本學習筆記如有錯誤,並非有意,但仍然歡迎在討論去狂 \(\sf D\) 她。

另:本文各個部分的講解順序 並不按照難度排序,而是按照我學習的順序,請選擇合適的內容觀看。

關於圖片:本文所有圖片均為作者純手畫 是 rui_er 做夢時託夢給我的(

\(\sf Part\ 2\) 何為計算幾何

計算幾何,就是利用計算機建立數學模型解決幾何問題。

要用電腦解幾何題?數學好的同學們笑了。

要用電腦解幾何題?rui_er 笑了,我做夢就解決了(不是

我們並不是用計算機算數學卷子上的幾何題去了,而是解決一些更加複雜的幾何相關問題。

為了解決複雜且抽象的問題,我們一定要選擇合適的研究方法。對於計算機來說,給它看幾何圖形……

\(\sf Part\ 3\) 凸包

\(\sf Part\ 3.1\) 二維凸包

\(\sf Part\ 3.1.1\) 凸多邊形

凸多邊形是指所有內角大小都在 \([0, \pi]\) 範圍內的 簡單多邊形

\(\sf Part\ 3.1.2\) 凸包

在平面上能包含所有給定點的最小凸多邊形叫做凸包。

其定義為:對於給定集合 \(X\) ,所有包含 \(X\) 的凸集的交集 \(S\) 被稱為 \(X\)凸包

\(\qquad\qquad\) —— OI-Wiki

其實我們可以把凸包看成一個拿橡皮筋圍成的一個圖形。

假設有一個佈滿小凸起的板子:

我們要把這些凸起都圍起來,怎麼圍呢?

顯然,最簡單方變的方法是這樣:

但是,我們知道,橡皮筋是有彈力的,所以橡皮筋會往裡縮,直到這樣:

最外圈的凸起撐起了橡皮筋。

此時橡皮筋圍成的多邊形的頂點就是最外圈凸起所在的位置。

由此,我們就定義橡皮筋圍成的圖形為一個平面凸包。

那麼,換一種定義,就為:

平面凸包是指覆蓋平面上 \(n\) 個點的最小的凸多邊形。

當然,我們發現在程式中卻無法模擬橡皮筋收縮的過程,於是有了下文的誕生。

\(\sf Part\ 3.1.3\) 二維凸包的求法

\(\sf Part\ 3.1.3.1\) 斜率逼近法

其實這也是一種容易想到的演算法,但是並不常用(程式碼複雜度高),所以我們省略帶過。

\(\sf Part\ 3.1.3.2\) Jarvis 演算法

這其實是一種數學構造法

此演算法的時間複雜度為 \(O(nm)\)

但是,這個複雜度會爆炸。

於是我們偉大的 \(\sf OIer\) 就想到了 Graham 演算法。

\(\sf Part\ 3.1.3.3\) Graham 演算法

Graham 演算法的本質:

Graham 掃描演算法維護一個凸殼,通過不斷在凸殼中加入新的點和去除影響凸性的點,最後形成凸包。

凸殼:凸包的一部分。

此演算法主要分為兩部分:

  • 排序
  • 掃描
\(\sf Part\ 3.1.3.3.1\) 排序

我們先選擇一個 \(y\) 最小的點(如 \(y\) 相同選 \(x\) 最小),記為 \(p_1\)

剩下的點,按照極角的大小逆時針排序,記為 \(p_2,p_3,\dots, p_m\)

我們按照排序結束時的順序列舉每一個點,依次連線,這裡可以使用一個棧來儲存,每次入棧,如果即將入棧的元素與棧頂兩個元素所構成了一個類似於凹殼的東西,那麼顯然處於頂點的那個點一定不在這個點集的凸包上,而他正好在棧頂,所以把它彈出棧,新點入棧。——節選自 數論小白都能看懂的平面凸包詳解 - ShineEternal的部落格

\(\sf Part\ 3.1.3.3.2\) 掃描

(下列所說的左右等是指以上一條連線為鉛垂線,新的連線偏移的方向)

剛開始,我們的點集是這樣的:

\(p_1\) 為起始點。

隨後,\(p_2\) 準備入棧,由於棧元素很少,所以可以入棧。

再看 \(p_3\),因為 \(p_3\) 向左,符合凸包條件,入棧。

隨後 \(p_4\) 也一切正常,依然向左,入棧。

\(p_5\) 依然向左,入棧。

\(p_6\) 時,我們發現了點問題,就是不再是向左了,而是向右了,所以我們此時要將 \(p_5\) 出棧,\(p_6\) 入棧。

入棧後,我們發現,相對於 \(p_4\)\(p_6\) 依然是向右的,所以我們還要把 \(p_4\) 出棧,\(p_6\) 入棧。

接下來 \(p_7\) 沒有問題。

\(p_8\) 時,我們發現,也是向右的,所以將 \(p_7\) 出棧,\(p_8\) 入棧。

接下來 \(p_9\) 正常,入棧。

最後,我們再把最後一個點和第一個點連起來。

此時,我們的 Graham 演算法的全過程就結束了。

掃描的時間複雜度:\(O(n)\)

但是不可能那麼優秀,還要加上排序的時間複雜度:\(O(n\log n)\)

所以總時間複雜度為 \(O(n \log n)\)

可見此演算法相比之前的演算法優秀的多。

\(\sf Part\ 3.1.3.4\) Andrew 演算法

假設我們有這些點:

首先把所有點以橫座標為第一關鍵字,縱座標為第二關鍵字排序。

相對於 Graham 演算法來說,Andrew 演算法排序更簡單,按 \(x, y\) 座標排序,時間複雜度也更低(一般的座標系中排序方法)。

首先將 \(p_1\) 入棧。

然後也將 \(p_2\) 入棧,\(p_2\) 可能在,也可能不在,等著之後判斷。

隨後,發現 \(p_3\) 偏右,所以我們將 \(p_2\) 出棧。

發現 \(p_4\) 依然偏右,\(p_3\) 出棧,\(p_4\) 入棧。

\(p_5\) 向右,\(p_4\) 出棧,\(p_5\) 入棧。

\(p_6\) 向左,入棧。

\(p_7\) 向右,\(p_6\) 出棧,\(p_7\) 入棧。

\(p_8\) 向右,\(p_7\) 出棧,繼續檢查發現相對於 \(p_5\) \(p_8\) 仍然向右,\(p_5\) 出棧,\(p_8\) 入棧。

此時,我們發現,凸包明明還空一半就到頭了???

然而這是意料之中,我們這種演算法必然會只算出一半的凸包。

所以我們需要再從排序末尾的點(也就是 \(p_8\))出發,按照一模一樣的方式再算一遍就行了。

當然如果我們走過的點就不許要再走了(除了 \(p_1\)

\(p_8\)\(p_7\),向左,\(p_7\) 入棧。

\(p_6\) 向右,\(p_7\) 出棧,\(p_6\) 入棧。

\(p_5\) 向左,入棧。

\(p_4\) 向左,入棧。

\(p_3\) 向右,\(p_4\) 出棧,對於 \(p_5\) \(p_3\) 依然向右,\(p_5\) 出棧,\(p_3\) 入棧。

\(p_2\) 向右,\(p_3\) 出棧,\(p_2\) 入棧。

最後將 \(p_2\)\(p_1\) 連起來。

至此,我們的 Andrew 演算法就完成了!

掃描的時間複雜度:\(O(n)\)(已過濾常數)

排序時間複雜度:\(O(n \log n)\)

總時間複雜度:\(O(n \log n)\)

\(\sf Part\ 3.1.4\) 實戰演練

\(\sf Part\ 3.1.4.1\) P2742 [USACO5.1]圈奶牛Fencing the Cows /【模板】二維凸包

先拿模板題練練手。

題目簡述:求一個點集凸包的邊長和。

拿 Graham 演算法做即可。

直接套模板即可,手寫棧應該更快(?

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<string>
#define line cout << endl
using namespace std;
const int NR = 1e5 + 5;
int n;
double ans;
struct point {
	double x, y;
};
point p[NR], ps[NR];
double dis (point pa, point pb) { //求兩點間距離
	return sqrt ((pb.x - pa.x) * (pb.x - pa.x) + (pb.y - pa.y) * (pb.y - pa.y));
}
double cp (point pa1, point pa2, point pb1, point pb2) { //求叉積
	return (pa2.x - pa1.x) * (pb2.y - pb1.y) - (pb2.x - pb1.x) * (pa2.y - pa1.y);
}
bool cmp (point px, point py) { //排序
	double num = cp (p[1], px, p[1], py);
	if (num > 0) return true;
	if (num == 0 && dis (p[0], px) < dis (p[0], py)) return true;
	return false;
}
int main () {
	cin >> n;
	for (int i = 1; i <= n; i++) {
		cin >> p[i].x >> p[i].y;
		if(i != 1 && p[i].y < p[1].y) { //去重
			swap (p[i].y, p[1].y);
			swap (p[i].x, p[1].x);
        }
	}
	sort (p + 2, p + n + 1, cmp);
	ps[1] = p[1]; //最低點是肯定在凸包裡的
	int h = 1;
	for (int i = 2; i <= n; i++) {
		while (h > 1 && cp (ps[h - 1], ps[h], ps[h], p[i]) <= 0) { //判斷是向左還是向右,如果向右就出棧
			h--;
		}
		h++;
		ps[h] = p[i];
	}
	ps[h + 1] = p[1]; //最後一個點跟第一個點相連
	for (int i = 1; i <= h; i++) {
		ans += dis (ps[i], ps[i + 1]); //累加
	}
	printf ("%.2lf\n", ans);
	return 0;
}

\(\sf Part\ 3.1.4.2\) UVA11626 Convex Hull

拿 Andrew 練練手吧qaq,也是道模板題,但是好像拿 Graham 會 TLE(?

/*
  Problem:UVA11626
  Dpte:03/07/20 14:05
*/
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<string>
#define line cout << endl
using namespace std;
const int NR = 1e5 + 5;
const double eps = 1e-7;
int n;
struct point {
    double x, y;
    point () {}
    point (double a, double b) : x (a), y (b) {}
    bool operator < (const point &b) const {
        if (x < b.x) return 1;
        if (x > b.x) return 0;
        return y < b.y;
    }
    point operator - (const point &b) {
        return point (x - b.x, y - b.y);
    }
};
point p[NR], sp[NR];
int cmp (double x) {
    if (fabs (x) < eps) return 0;
    return x > 0 ? 1 : -1;
}
double dis (point a, point b) {
    return sqrt ((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
double cp (point a, point b) {
    return a.x * b.y - a.y * b.x;
}
int Andrew () {
    sort (p + 1, p + 1 + n);
    int len = 0;
    for (int i = 1; i <= n; i++) {
        while (len > 1 && cmp (cp (sp[len] - sp[len - 1], p[i] - sp[len - 1])) < 0) 
            len--;
        sp[++len] = p[i];
    }
    int k = len;
    for (int i = n - 1; i >= 1; i--) {
        while (len > k && cmp (cp (sp[len] - sp[len - 1], p[i] - sp[len - 1])) < 0)
            len--;
        sp[++len] = p[i];
    }
    return len;
}
int main () {
    int t;
    cin >> t;
    while (t--) {
        cin >> n;
        char c;
        for (int i = 1; i <= n; i++)
            cin >> p[i].x >> p[i].y >> c;
        int t = Andrew();
        cout << t - 1 << endl;
        for (int i = 1; i < t; i++) 
            printf ("%.0lf %.0lf\n", sp[i].x, sp[i].y);
    }
    return 0;
}

隨時更新...