演算法筆記 - Convex Hull

Convex Hull

中譯「凸包」或「凸殼」。在多維空間中有一群散佈各處的點，「凸包」是包覆這群點的所有外殼當中，表面積暨容積最小的一個外殼，而最小的外殼一定會是凸的。「凸包」圍住的區域，是各點做線性內插後的範圍。

至於「凸」的定義是：圖形內任意兩點的連線不會經過圖形外部，http://mathworld.wolfram.com/Convex.html。這裡指的「凸」並不是表面弧狀凸起之意，事實上凸包是由許多平坦表面組成的。

以下討論比較簡單的情況：從二維平面上散佈的點當中找出凸包。二維平面上的凸包會是一個凸多邊形，在所有點的外圍繞一圈即得凸包。另外，最頂端、最底端、最左端、最右端的點，一定是凸包上的點。

計算凸包時需考慮一些特殊情況：一、凸包上多點重疊；二、凸包上多點共線；三、凸包呈一條線段、一個點、沒有點。通常我們會簡化資訊，以最少的點來記錄凸包，去掉重疊、共線的點。

UVa 109 132 218 361 675 681 811 819 10002 10065 10078 10135 10173 10256 10625 11168 11626

UVa 802 10089

演算法

從一個凸包上的頂點開始，順著外圍繞一圈，順時針或逆時針都可以。

每當尋找下一個要被包覆的點，則窮舉平面上所有點，找出最外圍的一點來包覆（可利用外積運算來判斷）。

時間複雜度為O(N*M)，N為所有點的數目，M為凸包的頂點數目。

// P為平面上散佈的點。設定為10點。
// CH為凸包上的頂點。設定為逆時針方向排列。
struct Point {int x, y;} P[10], CH[10];

// 小於。用以找出最低最左邊的點。
bool compare(const Point& a, const Point& b)
{
	return (a.y < b.y) || (a.y == b.y && a.x < b.x);
}

// 向量OA外積向量OB。大於零表示從OA到OB為逆時針旋轉。
int cross(Point& o, Point& a, Point& b)
{
	return (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x);
}

// 兩點距離的平方
int length2(Point& a, Point& b)
{
	return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}

// 以o點作為中心點，a點比較遠，b點比較近。
bool far(Point& o, Point& a, Point& b)
{
	return length2(o, a) > length2(o, b)
}

void Jarvis_march()
{
	/* 起點必須是凸包上的頂點。這裡用最低最左邊的點作為起點。 */

int start = 0;
	for (int i=0; i<10; ++i)
		if (compare(P[i], P[start]))
			start = i;

/* 包禮物，逆時針方向。 */
	int m = 0;			// m 為凸包頂點數目
	CH[m] = P[start];	// 紀錄起點

for (int m=1; true; ++m)
	{
		/* 找出位於最外圍的一點，若凸包上有多點共線則找最遠的一點。 */

int next = start;
		for (int i=0; i<10; ++i)
		{
			int c = cross(CH[m-1], P[i], P[next]);
			if (c > 0 ||
				c == 0 && far(CH[m-1], P[i], P[next]))
				next = i;
		}

if (next == start) break;	// 繞一圈後回到起點了
		CH[m] = P[next];			// 紀錄方才所找到的點
	}
}

實作時盡量避免：將一筆規模較大的資料，複製多份，存放多處。上面的程式碼是不好的示範，實作時請多多運用指標或索引值。

如果想找出凸包上重疊的點與共線的點，則改為尋找離上一點最近的點，並且標記目前已包過的點。

演算法

由前面段落可知：凸包上的頂點很有順序的沿著外圍繞行一圈，這個順序是時針順序。

Graham's Scan嘗試先將所有點按照時針順序排好，再來做繞一圈的動作，繞行途中順便掃除凸包內的點，如此就不必以窮舉所有點的方式來尋找最外圍的點。

要讓各點依時針順序排好，只要將中心點設定在凸包內部或設定在凸包上面，然後讓各點依角度排序即可。如果把中心點設定在凸包外部，結果就不見得是時針順序了。

一般來說，選擇凸包上面的頂點當作角度排序的中心點是比較好的，因為兩點間最大的夾角必會小於180度，可以使用外積運算來排序（外積在大於180度時會得負值、等於180度時會等於零，導致排序錯誤。）另一個好處是，中心點也可以作為包覆的起點。

包覆時必須按照時針順序，才能確保結果正確。錯誤例子是：給定一個簡單多邊形（simple polygon），直接照連接順序包覆，忘了做角度排序。正確例子是：給定一個星狀多邊形（star-shaped polygon），可以直接照連接順序包覆，不必做角度排序。

如果想找出凸包上重疊的點與共線的點，便要小心處理剛開始包、快要包好時產生共線的情形。相當麻煩。

有一個解決方法是，當角度排序遇到角度相同的情況，則讓距離中心點較近的點排前面，然後包覆時由起點開始同時往左右兩邊包。下面這段程式碼寫出一些特別要注意的地方：

至於凸包退化成線段或點的情況，則更難解決，此處不討論。

時間複雜度

尋找起點的時間，加上排序的時間，加上包覆的時間。

尋找起點需時O(N)。排序需時O(NlogN)，用Counting Sort之類的排序法可達到O(N)。包覆需時O(N)：總共前進N次，最多倒退N次，共為2N次。

演算法

前面段落採用的是時針順序，此處改為把所有點依座標大小排序：http://www.algorithmist.com/index.php/Monotone_Chain_Convex_Hull。我也找到了有趣的Applet：http://wind.lcs.mit.edu/~aklmiu/6.838/convexhull/index.html。一口氣解決了凸包有重疊的點、共線的點、退化成線段和點的問題，是相當優美的演算法。

以幾何的角度來看，這個演算法可以想做是Graham's Scan當中，做為排序基準的點，被放在無限遠處，於是夾角大小變成了水平距離差，時針順序變成了水平座標順序。

時間複雜度為排序的時間加上包覆的時間，為O(NlogN)。

演算法

這是online演算法，隨時維護一個凸包。每當輸入一點，如果此點在凸包內部就直接捨棄；不然就計算此點與當前凸包的兩條切線，然後更新凸包。

要找切線，窮舉切點即可。切點的左右鄰點都在切線同側，反之則否，判斷僅需O(1)時間。要小心切線與邊重疊的情況。

凸包的資料結構採用circular linked list，找到兩個切點後，移除其間的連續凹點僅需O(1)時間。總時間複雜度是O(N^2)。

改進

換個角度想，想要找到新凸包，直接窮舉新凸包的頂點不就好了？

以試誤法嘗試舊凸包的每個頂點。以當前輸入點為基準，若為凸面、切點，則留下；若為凹面，則捨棄。最後就得到新凸包。

如此一來就不需要特殊資料結構了。時間複雜度是O(N^2)。

加速

以當前輸入點為基準，切點的一側是凸面，另一側是凹面，切點恰為凹凸分際。故可用Binary Search找切點。

凸包的資料結構可以採用Splay Tree，找切點、移除連續頂點僅需O(logN)均攤時間。總時間複雜度是O(NlogN)。

Splay Tree作排序時，可以參考凸包最左下點，以角度排序。

預先排序

如果預先按照XY座標排序所有點（平移的掃描線），此演算法便與Andrew's Monotone Chain大同小異，時間複雜度變成O(NlogN)。

另外，預先排序之後，當前輸入點必定都在凸包外部（或是凸包頂點上），必定擁有兩條切線。

演算法

一開始將所有點以X座標位置排序。
Divide：將所有點分成左半部和右半部。
Conquer：左半部和右半部分別求凸包。
Merge：將兩個凸包合併成一個凸包。

合併凸包，找到兩條外公切線即可。從左凸包最右點、右凸包最左點開始，固定左端順時針轉、固定右端逆時針轉，輪流前進直到卡死，就得到下公切線，時間複雜度為O(N)。

注意到，下公切線並不見得是左右兩凸包的最低點連線，所以就算一開始抓的是左右凸包的最低點，運氣不好時，仍需要O(N)時間才能找到下公切線。況且抓最低點有可能已經衝過頭了。

附帶一提，內公切線也可如法炮製，改為從左凸包最左點、右凸包最右點開始。如果一個取內側、一點取外側，找公切線有可能衝過頭。

至於正確性證明：找外公切線的過程中絕對不會與兩凸包相交、找內公切線的過程中一定會與兩凸包相交，卡死時即是相交與不相交的分際，分際即是公切線。

也許可以仿照前述Incremental Method的凹凸分際概念，加速至O(logN)？這就留給讀者思考了。

時間複雜度為下述兩項總和：一、一次排序的時間，通常為O(NlogN)；二、Divide and Conquer向下遞迴O(logN)深度，合併凸包需要O(N)時間，總共為O(NlogN)。

不預先排序

一開始可以不必排序，只要把所有點分成兩等份即可。兩個凸包合併成一個凸包時，兩個凸包可能會重疊，仍然可以用O(N)時間解決，不過演算法較複雜，此處省略之。

http://westhoffswelt.de/blog/0040_quickhull_introduction_and_php_implementation.html

http://taibif.org.tw/informatics/?p=439

時間複雜度是O(N^2)。當點呈隨機分布時，時間複雜度是O(N)。是現實中最有效率的演算法。

一、找出左端點與右端點，連線，所有點分為上半部與下半部。
　　上半部與下半部分開求解。
二、處理上半部，找出上半凸包：
　甲、上端點、左端點、右端點形成一個三角形區域。
　　　三角形內部的點，不會是凸包上的點，可剔除之。
　乙、運用外積運算，找出三角形外的點，分為左、右兩份。
　丙、左、右兩份分別遞迴求解，直到找出上半凸包。
三、處理下半部，找出下半凸包。

演算法

原理是Trial and Error與Divide and Conquer。這是一個複合式的演算法，使用到了Graham's Scan、Jarvis' March，另外也用到了Divide and Conquer的技巧。

一、假設凸包有M個點。
二、使用試誤法，依序嘗試2^(2^0)、2^(2^1)、2^(2^2)、……這些數值做為M。
　甲、每M個點為一組，所有點被分作⌈N/M⌉組。
　　　O(N)。
　乙、每一組各自求凸包，一共得到⌈N/M⌉個凸包。《Graham's Scan》
　　　O(MlogM * ⌈N/M⌉) = O(NlogM)。
　丙、嘗試求出這些凸包的凸包。《Jarvis' March》
　　　O(3 * ⌈N/M⌉ * M) = O(N)。
　　子、每個凸包各用一個指標，指著各自的最低點。
　　　　O(N)。
　　丑、找出所有凸包的最低點，從最低點開始包覆。
　　　　O(N)。
　　寅、每當尋找下一個要被包覆的點：
　　　回、各凸包各自找出最靠外面的一點，一共得到⌈N/M⌉個點。
　　　　　由指標處繼續往後找，指標只進不退。要預覽下一點。
　　　　　O(3 * ⌈N/M⌉)，均攤。
　　　回、再從這⌈N/M⌉個點當中，看看哪一點最靠外面。
　　　　　O(⌈N/M⌉)。
　　卯、若包了M個點還未形成凸包，則馬上停止，回到步驟二！
　　　　如果途中形成凸包，即是正解。演算法結束。

時間複雜度

每次執行步驟二的時間複雜度是O(NlogM)，M是會變動的。對M進行試誤時，謹慎的選擇M的數值，可以將整體的時間複雜度強壓在O(NlogM)以內。

對M進行試誤時，
用0、1、2、……，整體的時間複雜度為：
O(Nlog0 + Nlog1 + Nlog2 + ... + NlogM) = O(N * logM * M)

用2^0、2^1、2^2、……，整體的時間複雜度為：
O(N*0 + N*1 + N*2 + ... + NlogM) = O(N * logM * logM)

用2^(2^0)、2^(2^1)、2^(2^2)、……，整體的時間複雜度為：
O(N*1 + N*2 + N*4 + ... + NlogM) = O(N * logM)

直接用N，馬上能找出解，不過整體的時間複雜度，不是我們所要的：
O(NlogN)

選擇M的原理，其實就是猜數字的第一個步驟，請參考「Binary Search延伸閱讀：猜數字」。

整體的時間複雜度是O(NlogM)，N為所有點的數目，M為凸包的頂點數目，是目前時間複雜度最低的演算法。然而實際執行起來，通常會比先前那些演算法還慢。

演算法

求出一簡單多邊形的凸包。

http://cgm.cs.mcgill.ca/~athens/cs601/Melkman.html

時間複雜度為O(N)。是相當優美的演算法。

Convex Layers（Onion）

由外往內，一層一層的凸包，每一層都是一個Convex Layer。全部的Convex Layer合稱為一個「洋蔥」。

最內部的Convex Layer可能退化成一個點或是一條線段。

演算法

使用Jarvis' March一直繞圈，時間複雜度為O(N^2)。

排序所有點之後，不斷找出剩餘諸點的凸包，最多找O(N)次凸包。可以採用Graham's Scan或者Andrew's Monotone Chain。總時間複雜度為O(NlogN) + O(N^2) = O(N^2)。

ICPC 3655