Convex Hull

http://www.csie.ntnu.edu.tw/~u91029/ConvexHull.html#4

Convex Hull

Convex Hull

中譯「凸包」或「凸殼」。在多維空間中有一群散佈各處的點，「凸包」是包覆這群點的所有外殼當中，表面積暨容積最小的一個外殼，而最小的外殼一定是凸的。

至於「凸」的定義是：圖形內任意兩點的連線不會經過圖形外部。「凸」並不是指表面呈弧狀隆起，事實上凸包是由許多平坦表面組成的。

以下討論比較簡單的情況：替二維平面上散佈的點，找到凸包。

二維平面上的凸包是一個凸多邊形，在所有點的外圍繞一圈即得凸包。另外，最頂端、最底端、最左端、最右端的點，一定是凸包上的點。

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

UVa 109 132 218 361 596 675 681 811 81910002 10065 10078 10135 10173 10256 10625 1116811626 ICPC 4450

UVa 802 10089 ICPC 7585

任意圖形的凸包

任意圖形都能求出凸包。例如一個多邊形的凸包、大量三角形的凸包、大量線段的凸包。這些問題都可以簡化為一群點的凸包。

圓形、曲線的凸包，我們以下不討論。

Convex Hull: 
Jarvis' March 
（ Gift Wrapping Algorithm ）

演算法

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

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

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

 

1. // P為平面上散佈的點。設定為10點。
2. // CH為凸包上的頂點。設定為逆時針方向排列。
3. struct Point {int x, y;} P[10], CH[10];
4.  
5. // 小於。用以找出最低最左邊的點。
6. bool compare(const Point& a, const Point& b)
7. {
8.     return (a.y < b.y) || (a.y == b.y && a.x < b.x);
9. }
10.  
11. // 向量OA叉積向量OB。大於零表示從OA到OB為逆時針旋轉。
12. int cross(Point& o, Point& a, Point& b)
13. {
14.     return (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x);
15. }
16.  
17. // 兩點距離的平方
18. int length2(Point& a, Point& b)
19. {
20.     return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
21. }
22.  
23. // 以o點作為中心點，a點比較遠，b點比較近。
24. bool far(Point& o, Point& a, Point& b)
25. {
26.     return length2(o, a) > length2(o, b)
27. }
28.  
29. void Jarvis_march()
30. {
31.     /* 起點必須是凸包上的頂點。這裡用最低最左邊的點作為起點。 */
32.  
33.     int start = 0;
34.     for (int i=0; i<10; ++i)
35.         if (compare(P[i], P[start]))
36.             start = i;
37.  
38.     /* 包禮物，逆時針方向。 */
39.     int m = 0;          // m 為凸包頂點數目
40.     CH[m] = P[start];   // 記錄起點
41.  
42.     for (int m=1; true; ++m)
43.     {
44.         /* 找出位於最外圍的一點，若凸包上有多點共線則找最遠的一點。 */
45.  
46.         int next = start;
47.         for (int i=0; i<10; ++i)
48.         {
49.             int c = cross(CH[m-1], P[i], P[next]);
50.             if (c > 0 ||
51.                 c == 0 && far(CH[m-1], P[i], P[next]))
52.                 next = i;
53.         }
54.  
55.         if (next == start) break;   // 繞一圈後回到起點了
56.         CH[m] = P[next];            // 記錄方才所找到的點
57.     }
58. }

實作時請多多運用指標、索引值，儘量避免搬動原本資料。此處的程式碼是不良示範，僅供參考。

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

Convex Hull: 
Graham's Scan

演算法

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

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

要讓所有點依照時針順序排好，只要將中心點設定在凸包內部或者凸包上，然後讓各點依照角度排序即可。如果把中心點設定在凸包外部，結果就不見得是時針順序了。包覆時必須按照時針順序，才能確保結果正確。

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

角度排序時，遇到角度相同的情況，要小心排序。通常是讓距離中心點較近的點排前面。也可以排後面，但是不能亂排。

掃除的過程當中，經常株連許多點。使用 stack 資料結構來儲存凸包，逐一判斷 stack 頂端的點，逐一彈出凹陷的點。凹陷的點必定不是凸包上的頂點。

 

1. // P為平面上散佈的點。設定為10點。
2. // CH為凸包上的頂點。設定為逆時針方向排列。可以視作一個stack。
3. struct Point {int x, y;} P[10+1], CH[10+1];
4.  
5. // 向量OA叉積向量OB。大於零表示從OA到OB為逆時針旋轉。
6. int cross(Point& o, Point& a, Point& b)
7. {
8.     return (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x);
9. }
10.  
11. // 小於。用以找出最低最左邊的點。
12. bool compare_position(const Point& a, const Point& b)
13. {
14.     return (a.y < b.y) || (a.y == b.y && a.x < b.x);
15. }
16.  
17. /*
18. // 小於。以P[0]當中心點做角度排序，角度由小排到大（即逆時針方向）。
19. // 角度相同時，順序隨便。
20. bool compare_angle(const Point& a, const Point& b)
21. {
22.     return cross(P[0], a, b) > 0;
23. }
24. */
25.  
26. // 小於。以P[0]當中心點做角度排序，角度由小排到大（即逆時針方向）。
27. // 角度相同時，距離較離中心點較近的點排前面。
28. bool compare_angle(Point& a, Point& b)
29. {
30.     // 加入角度相同時，距離長度的判斷。
31.     int c = cross(P[0], a, b);
32.     return c > 0 ||
33.             (c == 0 && length2(P[0], a) < length2(P[0], b));
34. }
35.  
36. void Graham_scan()
37. {
38.     // 這裡用最低最左邊的點作為起點，並將中心點換到第零點。O(N)
39.     swap(P[0], *min_element(P, P+10, compare_position));
40.  
41.     // 其餘各點依角度排序。O(NlogN)
42.     sort(P+1, P+10, compare_angle);
43.  
44.     // 直接把中心點作為起點，開始包覆，逆時針方向。O(N)
45.     P[N] = P[0];    // 讓程式方便處理包至最後一點的情況。
46.  
47.     int m = 0;      // m 為凸包頂點數目
48.     for (int i=0; i<=10; ++i)
49.     {
50.         // 掃除凹陷的點
51.         while (m >= 2 && cross(CH[m-2], CH[m-1], P[i]) <= 0) m--;
52.         // 添加新的點
53.         CH[m++] = P[i];
54.     }
55.  
56.     m--;    // 最後一個點是重複出現的起點，故要減一。
57. }

如果想找出凸包上重疊的點、共線的點，必須特別小心處理剛開始包、快包好時那些重疊的點、共線的點。一種解決方式是：從最低最左點開始，分頭往左右兩邊包，包至最高最右點。相當麻煩。

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

總而言之，此演算法無法漂亮的解決多點共線的問題。

時間複雜度

尋找起點的時間 O(N) 。加上排序的時間 O(NlogN) 。加上包覆的時間 O(N) ：總共前進 N 次，最多倒退 N 次，共為 2N 次。

全部加起來是 O(NlogN) ，主要取決於排序的時間。

星狀多邊形

星狀多邊形（ star-shaped polygon ）的定義是：多邊形內部存在一個點，可以看到整個多邊形內部。也就是說，星狀多邊形可以直接依照連接順序包覆，不必再做角度排序。

Convex Hull: 
Andrew's Monotone Chain

演算法

前面章節採用時針順序，此處改為座標順序。以 X 座標大小排序，當 X 座標相同則以 Y 座標大小排序。

先從起點開始，按照順序掃描，找到下半凸包。再從終點開始，按照相反順序掃描，找到上半凸包。合起來就是完整的凸包。

一口氣解決了凸包有重疊的點、共線的點、退化成線段和點的問題，是相當優美的演算法。

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

 

1. // P為平面上散佈的點。設定為10點。
2. // CH為凸包上的頂點。設定為逆時針方向排列。可以視作一個stack。
3. struct Point {int x, y;} P[10], CH[10+1];
4.  
5. // 向量OA叉積向量OB。大於零表示從OA到OB為逆時針旋轉。
6. double cross(Point& o, Point& a, Point& b)
7. {
8.     return (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x);
9. }
10.  
11. // 小於。依座標大小排序，先排 x 再排 y。
12. bool compare(Point& a, Point& b)
13. {
14.     return (a.x < b.x) || (a.x == b.x && a.y < b.y);
15. }
16.  
17. void Andrew_monotone_chain()
18. {
19.     // 將所有點依照座標大小排序
20.     sort(P, P+10, compare);
21.  
22.     int m = 0;  // m 為凸包頂點數目
23.  
24.     // 包下半部
25.     for (int i=0; i<10; ++i)
26.     {
27.         while (m >= 2 && cross(CH[m-2], CH[m-1], P[i]) <= 0) m--;
28.         CH[m++] = P[i];
29.     }
30.  
31.     // 包上半部，不用再包入方才包過的終點，但會再包一次起點
32.     for (int i=10-2, t=m+1; i>=0; --i)
33.     {
34.         while (m >= t && cross(CH[m-2], CH[m-1], P[i]) <= 0) m--;
35.         CH[m++] = P[i];
36.     }
37.  
38.     m--;    // 最後一個點是重複出現兩次的起點，故要減一。
39. }

 

1. // 點、下凸包、上凸包
2. struct Point {int x, y;} P[10], L[10], U[10];
3.  
4. void Andrew_monotone_chain()
5. {
6.     sort(P, P+10, cmp);
7.  
8.     int l = 0, u = 0;   // 下凸包的點數、上凸包的點數
9.     for (int i=0; i<10; ++i)
10.     {
11.         while (l >= 2 && cross(L[l-2], L[l-1], P[i]) <= 0) l--;
12.         while (u >= 2 && cross(U[u-2], U[u-1], P[i]) >= 0) u--;
13.         L[l++] = P[i];
14.         U[u++] = P[i];
15.     }
16. }

Convex Hull: 
Quick Hull Algorithm

演算法

一、找出最左點與最右點，連線，所有點分為上半部與下半部。
　　上半部與下半部分開求解。
二、處理上半部，找出上半凸包：
　甲、找到距離底線最遠的點，形成一個三角形區域。
　　　（如果有許多點，就找最靠近左端點的那一點。避免凸包上多點共線。）
　乙、運用叉積運算，找出三角形外部的點，分為左、右兩份。
　　　至於三角形內部、三角形上的點，不會是凸包頂點，盡數剔除之。
　丙、左、右兩份分別遞迴求解，直到找出上半凸包。
　　　三角形的兩翼，分別做為左、右兩份的底線。
三、處理下半部，找出下半凸包。

儘管時間複雜度是 O(N^2) ，卻是實務上最有效率的演算法。此演算法的好處是：預先剔除凸包內部大部分的點，而且不必預先排序所有點。

排序的時間複雜度是 O(NlogN) ，空間複雜度是 O(N) 。當點的數量很大時，例如一千萬，時間約是一千萬的 23 倍。又由於記憶體不足，必須採用外部排序，需要更多時間。

此演算法一開始掃描一次所有點，找到三角形外部的點，時間約是一千萬的 1 倍。如果三角形外部的點很少，例如一千點，那麼接下來的步驟得以節省許多時間。因此，總時間通常遠低於一千萬的 23 倍。

遞迴呼叫函數很花時間。當遞迴一兩次後，點數變少，此時可以直接套用其他凸包演算法，節省時間，例如 Andrew's Monotone Chain 。

Convex Hull: 
Divide and Conquer

演算法

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

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

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

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

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

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

不預先排序

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

Convex Hull: 
Incremental Method

演算法

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

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

凸包的資料結構採用 circular 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) 。

預先排序之後，當前輸入點必在凸包外部（點不重複時）、必有兩條切線。

預先隨機排列

隨機排列、排序，兩者概念相反。

預先隨機排列所有點，則第 i 回合固定新增兩點、平均掃描 N/i 點，平均時間複雜度為 O(NlogN) 。推廣到三維空間仍然如此。

Convex Hull: 
Chan's Algorithm

演算法

原理是 Trial and Error 與 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、……、M，整體的時間複雜度為：
O(Nlog0 + Nlog1 + Nlog2 + ... + NlogM) = O(N * logM * M)

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

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

用N，直接就找到答案，不過整體的時間複雜度，不是我們所要的：
O(NlogN)

選擇 M 的原理，其實就是「倍增搜尋」。

總時間複雜度是 O(NlogM) ， N 為所有點的數目， M 為凸包的頂點數目，是目前時間複雜度最低的演算法。然而實際執行起來，比先前介紹的演算法來得慢。

Convex Hull of Simple Polygon: 
Melkman's Algorithm

演算法

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

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

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

Convex Layers

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

◀ Index

Convex Hull

Convex Hull:
Jarvis' March
（Gift Wrapping Algorithm）

Convex Hull:
Graham's Scan

Convex Hull:
Andrew's Monotone Chain

Convex Hull:
Quick Hull Algorithm

Convex Hull:
Divide and Conquer

Convex Hull:
Incremental Method

Convex Hull:
Chan's Algorithm

Convex Hull of Simple Polygon:
Melkman's Algorithm

Convex Layers