[Multiple line segment intersection] Bentley–Ottmann 演算法

 基本上這個問題就是給你一些線段(格式通常為兩個端點)，你要找出這些線段的交點。直觀的做法兩兩進行計算會花上$\ord{n^2}$的時間，但大多數的情況下交點不會很多。為了解決這個問題，修改自Shamos–Hoey演算法的Bentley–Ottmann演算法可以在$\ord{(n+k)\log n}$的時間內找出所有交點，其中$k$是交點數量。

這裡附上實作時需要用到的基本資料結構：

演算法使用掃描線進行。掃描線是一條垂直線從左邊掃到右邊(有些實作是水平線從上面掃到下面)，並且在遇到事件點的時候進行相關處理。

線段的兩端點以及交點都作為事件點被紀錄在最終結果中。對於每個事件點$P$，我們會計算三個集合：

• U集合：所有以$P$為起始點的線段集合
• C集合：所有包含$P$的線段集合
• L集合：所有以$P$為結束點的線段集合

當然要先保證每條線段的起始點移動會在結束點的左方，只要得到線段後稍微判斷一下就可以做到了。每個事件點找出這三個集合後就可以很容易的判斷相交資訊，但要注意的是會有以下的退化情形：

• 線段退化成點：這種情況該點的U和L都會包含該線段。
• 兩線段重合：只有重合處的兩端點會被紀錄為事件點，可以根據UCL判斷出是否線段重合
• 垂直線段：排序點和線段時如果x一樣就按照y來比較
  

最後是掃描線的資料結構，需要一棵平衡的BST根據當前掃描線和各個線段切點的y值進行排序，但這件事是可以用STL做到的！我們把當前事件點傳進比較函數裡面進行計算，因為在任何一個時刻BST中的資料都是根據當前的比較函數由小排到大的，應該不算undefined behavior。另外該演算法的浮點數誤差很大，建議使用時套上處理誤差的模板或是直接用分數計算：

最後是測試的部分，以下圖做為測試範例：

將該圖轉換成我們接受的input如下：

10

-2 7 2 0

-2 7 -2 0

-2 6 2 5

-2 6 2 2

-2 4 2 7

-2 4 2 2

-2 4 4 1

-2 0 2 2

0 1 0 1

0 3 4 1

最後附上測試程式碼，需要的話可以自己執行看看：

[ Delaunay Triangulation ] Delaunay 三角剖分

我高中的時候在morris的網站第一次看到這個東西，那個時候就很想試著寫出來看看，但那時候英文很不好所以不會找資料就沒有成功。直到今年的NPSC，聽說台大出了一題找Voronoi Diagram的題目沒有人寫出來，於是我就想挑戰看看，寒假開始之後就認真研讀相關資料，終於在前幾天完成了隨機增量法以及分治法的模板。

先來說說甚麼是三角剖分 Triangulation 吧！假設平面上散布很多不重複的點，現在用一些線段將這些點連成許多的三角形，滿足所有三角形的內部不會有其他點，而且三角形的數量越多越好，這就是一個三角剖分。

Delaunay Triangulation是一種特殊的三角剖分，滿足所有的三角形其外接圓中不會出現其他的點，這樣的性質會使得所有三角形中最小的角盡量大，而且最重要的是它是 Voronoi Diagram 的對偶圖，Voronoi Diagram 直接求的話很容易有浮點誤差而且程式也不好寫，因此從Delaunay Triangulation直接轉換過去會更加方便。

還有一個很重要的點，透過Delaunay Triangulation 的邊集合可以找出平面歐幾里得距離最小生成樹。

這裡我先給出點的資料結構，方便接下來的講解：

一般來說找出 Delaunay Triangulation 的其中一個操作是要判斷一個點是不是在一個三角形的外接圓內，但直接去求外接圓圓心會有浮點誤差，因此有一個作法是將所有點投影到三維的拋物線曲面中，也就是下圖的碗狀曲面，會發現把一個圓投影上去之後，該圓會落在一個平面上，圓內的點會在平面下面，圓外的點會在平面上面，這樣就可以利用一些處理三維幾何的技術就可以避免浮點誤差的判斷點是否在外接圓內了。

範例程式碼中點在圓內輸出1，圓外輸出-1，圓上輸出0：

當然實作時並不會真正的寫一個三維點的資料結構，這樣程式碼會比較長。

Delaunay Triangulation 我所知道的有兩種，分別是隨機增量法以及分治法。隨機增量法再點隨機加入的時候複雜度是$\ord{n\log n}$，分治法則是穩定$\ord{n\log n}$，但我的實作分治法的效能是隨機增量法的三倍快。由於兩種做法都已經有詳細教學了這裡就不用文字敘述。

1. 隨機增量法教學在《Computational Geometry - Algorithms and Applications》書中有提到，該書有中文版
2. 分治法教學請點該網站連結

這裡附上程式碼，順帶一提morris的分治法好像在某些測資會因為三點共線爛掉

隨機增量法：

分治法：

[ lower convex hull self-balancing binary search tree ] 下凸包平衡樹

下凸包平衡樹是一個動態維護下凸包的資料結構，他至少能在$\ord{\log n}$時間完成以下兩個操作：

1. 查詢 point(x,y) 是否被包含於下凸包中
2. 將 point(x,y) 加入至下凸包中，然後將不位於下凸包頂點上的點刪除

由於繼承自C++ STL map，因此所有能對map用的操作都能被使用，對於那些會改動map資料的操作請小心使用避免影響到凸包的正確性，例如operator[]就盡量不要使用。

對於某些動態規劃問題需要用凸包斜率優化時應該非常有用，有空會補上專門用來動態做斜率優化的版本

[ closest pair of points ] 平面上歐基里德距離最近點對

用分治法寫的，複雜度$\ord{n\log^2 n}$:

[ dynamic kd tree ] 動態kd樹模板

參考資料:

• k-d tree教學投影片
• http://andrewd.ces.clemson.edu/courses/cpsc805/references/nearest_search.pdf
• https://cw.fel.cvut.cz/wiki/_media/courses/a4m33pal/paska13.pdf
• 在二維空間的查找優化
• K-D Tree在信息学竞赛中的应用
• K-D tree的估价
• 从K近邻算法、距离度量谈到KD树、SIFT+BBF算法 
• 感謝Morris大大的模板

kd樹的資料在網路上其實滿多的，但很多code效率不高或是功能不全，所以這邊稍微講解一下kd樹基本的一些操作

1. 構造
   現在有一個point序列S[0~n-1]，每個point都包含一個序列d[0~kd-1]表示維度，S[i].d[j]表示第i個點的第j維。
   節點的資料型別定義為node，裡面有l,r兩個指標分別指向左右兩顆子樹，還有一個point pid表示所存的point。
   構造的過程是一個遞迴結構，大致代碼如下:
   當前這個節點的k稱為分裂維度
   
2. 查找
   kd樹支援兩種查找法:
   一、給定一個點p和一個數字k，求離p前k近的點有哪些
   二、給定一個範圍，求範圍內有哪些點
   兩種方法跟一般爆搜有點像，但是利用了kd樹可以做到有效的剪枝，有時候可以判斷被切分出來的範圍內有沒有機會存在前k近點或是在範圍內，直接看模板應該就懂了
3. 刪除
   模板裡沒有附刪除的code所以會寫在這裡。
   首先如果找到要刪除的節點是葉節點直接刪除就好了；如果不是葉節點，假設現在這個點的分裂維度是k，就拿他右子樹第k維最小節點mi去替代他，接著遞迴刪除右子樹的mi；如果沒有右子樹，就拿他左子樹第k維最小節點mi去替代他，然後把左右子樹互換，接著遞迴刪除右子樹的mi。
   
   找到要刪除的點用查找中的第一種方法就好了，這裡p=要刪除的點，k=1，查找的時候順便維護最近節點其位置與其對p的距離，若最近距離!=0則刪除失敗。
   
   那對一個子樹找第k維最小節點呢?方法很簡單，也是遞迴定義的:
   首先如果當前節點o的分裂維度剛好是第k維，則若o有右子樹的話答案必在o的右子樹，否則答案為o，如果o的分裂維度!=k，則遞迴搜尋左右子樹，把得到的答案和o自己進行比較，求最小。
   接下來附上只有刪除操作的模板:
   
   

插入刪除通常配合替罪羊樹使用，插入很簡單模板有就不說了，接下來講一下複雜度:
設N=size(root)

• 構造:
  很明顯就是\(\ord{N \; log \; N}\)就不說了
• 查找:
  已經有人證明了，假設現在的維度有k維，則查找的最差複雜度是\(\ord{N^{1-1/k}}\)
  但是平均狀況下為\(\ord{log \; N}\)
• 插入:
  複雜度為\(\ord{樹的高度}\)，因為是套替罪羊樹，而且重構的時間是\(\ord{N \; log \; N}\)，所以單次插入的均攤時間複雜度是\(\ord{log \; N \; * \; log \; N}\)
• 刪除:
  有三個步驟需要分析，假設現在的維度有k維
  1. findmin:
     最多會遍歷\(\alpha^{(k-1)*(log_{\alpha}N)/k}=N^{1-1/k}\)個節點，所以是\(\ord{N^{1-1/k}}\)
     但實際操作量為\(N^{1-1/k}-n^{1-1/k}\)，n是最小值的子樹大小
  2. nearest:
     就是查找操作，所以是\(\ord{N^{1-1/k}}\)
     但因為是找相同點，可以優化code，所以實際操作量為\(N^{1-1/k}-n^{1-1/k}\)，n是相同點的子樹大小
  3. 刪除操作本身:
     看code明顯就是重複findmin直到把刪除的點變成葉節點為止，會感覺做了很多操作，但是我們發現把所有操作加起來後會=\(\ord{N^{1-1/k}}\)
     設\(o_1,o_2,...,o_d\)=findmin(root)找到的點,findmin(\(o_1\))找到的點,\(...\),要刪除的葉節點，\(n_1,n_2,...,n_d\)=findmin(root)找到的點的子樹大小,findmin(\(o_1\))找到的點的子樹大小,findmin(\(o_2\))找到的點的子樹大小,\(...\),要刪除的葉節點的子樹大小，\(d\)為樹的高度
     複雜度為:
     $$ N^{1-1/k}-n1^{1-1/k}+n1^{1-1/k}-n2^{1-1/k}+n2^{1-1/k}...-nd^{1-1/k}=\ord{N^{1-1/k}} $$

模板裡的code找的是曼哈頓距離，如果要找歐基里德距離只需要修改point裡的dist跟kd tree裡的heuristic即可，以下提供修改方法:
這樣查找時回傳的就是歐基里德距離的平方

以下附只有插入操作的模板:

以下附支援插入刪除的模板:
模板的使用方法請參考:

日月卦長的解題紀錄 [ IOICAMP2016 ] 動態曼哈頓最短距離

[ round-off error eps ] 浮點數誤差模板

有些計算幾何的過程中會產生捨位誤差，這時候就需要決定一個很小的數，去避開這個捨位誤差問題，我們稱那個很小的數為EPS

• A < B ....... A - B < - EPS
• A > B ....... A - B > EPS
• A == B ....... abs(A - B) <= EPS
• A <= B ....... A - B <= EPS
• A >= B ....... A - B >= - EPS

因為我的計算幾何模板在實作時是用template實作的，所以如果需要處理浮點數誤差時，需要傳入一個新的資料結構，所以寫了這份模板