超酷演算法：用四叉樹和希爾伯特曲線做空間索引

隨著越來越多的資料和應用和地理空間相關，空間索引變得愈加重要。然而，有效地查詢地理空間資料是相當大的挑戰，因為資料是二維的（有時候更高），不能用標準的索引技術來查詢位置。空間索引通過各種各樣的技術來解決這個問題。在這篇博文中，我將介紹幾種：四叉樹，geohash（不要和geohashing混淆）以及空間填充曲線，並揭示它們是怎樣相互關聯的。

四叉樹

四叉樹是種很直接的空間索引技術。在四叉樹中，每個節點表示覆蓋了部分進行索引的空間的邊界框，根節點覆蓋了整個區域。每個節點要麼是葉節點，有包含一個或多個索引點的列表，沒有孩子。要麼是內部節點，有四個孩子，每個孩子對應將區域沿兩根軸對半分得到的四個象限中的一個，四叉樹也因此得名。

圖1    展示四叉樹是怎樣劃分索引區域的 來源：維基百科

將資料插入四叉樹很簡單：從根節點開始，判斷你的資料點屬於哪個象限。遞迴到相應的節點，重複步驟，直到到達葉節點，然後將該點加入節點的索引點列表中。如果列表中的元素個數超出了預設的最大數目，則將節點分裂，將其中的索引點移動到相應的子節點中去。

圖2    四叉樹的內部結構

查詢四叉樹時從根節點開始，檢查每個子節點看是否與查詢的區域相交。如果是，則遞迴進入該子節點。當到達葉節點時，檢查點列表中的每一個項看是否與查詢區域相交，如果是則返回此項。

注意四叉樹是非常規則的，事實上它是一種字典樹，因為樹節點的值不依賴於插入的資料。因此我們可以用直接的方式給節點編號：用二進位制給每個象限編號（左上是00，右上是10等等 譯者注：第一個位元位為0表示在左半平面，為1在右半平面。第二個位元位為0表示在上半平面，為1在下半平面），任一節點的編號是由從根開始，它的各祖先的象限號碼串接而成的。在這個編號系統中，圖2中右下角節點的編號是1101。

如果我們定義了樹的最大深度，不需通過樹就可以計算資料點所在節點的編號：只要把節點的座標標準化到適當的整數區間中（比如32位整數），然後把轉化後x, y座標的位元位交錯組合。每對位元指定了假想的四叉樹中的一個象限。（譯者注：不瞭解的讀者可看看Z-order，它和下文的希爾伯特曲線都是將二維的點對映到一維的方法）

Geohash

上述編號系統可能看起來有些熟悉，沒錯，就是geohash！此刻，你可以把四叉樹扔掉了。節點編號，或者說geohash，包含了對於節點在樹中位置我們需要的全部資訊。全高樹中的每個葉節點是個完整的geohash，每個內部節點代表從它最小的葉節點到最大的葉節點的區間。因此，通過查詢所需的節點覆蓋的數值區間中的一切（在geohash上索引），你可以有效地定位任意內部節點下的所有資料點。

一旦我們丟掉了四叉樹，查詢就變得複雜一點了。我們需要事先構建搜尋集合而不是在樹中遞迴地精煉搜尋集合。首先，找到完全覆蓋查詢區域的最小字首（或者說四叉樹節點  譯者注：注意在我們的編號系統中節點由位元串表示）。在最壞情況下，這可能遠大於實際的查詢區域，比如對於在索引區域中心、和四個象限都相交的小塊地方，查詢將要從根節點開始。

現在的目標是構建一組完全包含查詢區域的字首，並且儘可能少包含區域外的部分。如果沒有其他約束，我們可以簡單地選擇與查詢區域相交的葉節點，但這會造成大量的查詢。所以要加一個約束：使得要查詢的不同區間最少。

一種達到這個目的的方法是先設定我們願意承受的查詢區間的最大數目。構建一組區間，最開始都設為我們之前指定的字首。從中選擇可以再分裂而不超出最大區間數並將從查詢區域刪除最不受歡迎區域的節點。重複這個過程直到集合中再沒有區間可以細分。最後，檢查得到的集合，如果可能的話合併相鄰的區間。下面的圖說明了這對於查詢一個圓形區域且限制最大5個查詢區間是如何工作的。

圖3    一個對區域的查詢是怎樣分解成一連串geohash字首/區間的

這個方法工作地很好，它使我們避免了遞迴查詢。我們執行的一整套區間查詢都可以並行完成。由於每次查詢都預期要一次硬碟搜尋，將查詢並行化大大減少了返回結果需要的時間。

然而，我們還可以做得更好。你可能注意到上圖中我們要查詢的所有區域都是相鄰的，但我們卻只能將其中兩個合併（選擇區域的右下角的兩個）成一個單獨的查詢，進而只要4次單獨查詢。（譯者注：這兩個區域可以合併是因為它們在geohash以Z字形遍歷區域的路徑上是相鄰的）這個後果部分是由於geohash訪問子區域的順序，在每個象限中從左到右，從上到下。從右上角象限到左下角象限的不連續性使得我們不得不將本可以使之連續的區間分裂。如果以不同的順序訪問區域，可能我們就可以最小化或者消除這些不連續性，使得更多的區域可以被看做是相鄰的，一次查詢就可得到結果。通過這樣效率上的提升，對於同樣的覆蓋區域，我們可以做更少的查詢，或者相反地，同樣的查詢次數的情況下包含更少的無關區域。

圖4    geohash訪問象限的順序

希爾伯特曲線

現在假設我們以U字形來訪問區域。在每個象限中，我們同樣以U字形來訪問子象限，但是要調整好U字形的朝向使得和相鄰的象限銜接起來。如果我們正確地組織了這些U字形的朝向，我們就能完全消除不連續性，不管我們選擇了什麼解析度，都能連續地訪問整個區域，可以在完全地探訪了一個區域後才移動到下一個。這個方案不僅消除了不連續性，而且提高了總體的局域性。按照這個方案得到的圖案看起來有些熟悉，沒錯，就是希爾伯特曲線。

希爾伯特曲線屬於一類被稱為空間填充曲線的一維分形，因為它們雖然是一維的線，卻可以填充固定區域的所有空間。它們相當有名，部分是由於XKCD把它們用於網際網路地圖。如你所見，對於空間索引它們也是有用的，因為它們展現的正是我們需要的局域性和連續性。再看看之前用一組查詢來覆蓋圓的例子，我們發現（應用希爾伯特曲線）還可以減少一次查詢：左下方的小區域現在和它右邊的區域連起來了（減少一次），雖然底部的兩塊區域不再連續了（增加一次），右下角的區域現在卻和它上方的連續了（減少一次）。

圖5    希爾伯特曲線訪問象限的順序

到目前為止，我們優雅的系統還缺一樣東西：將(x,y)座標轉換為希爾伯特曲線上相應位置的方法。對於geohash，這是簡單而明顯的–只需將x, y座標交錯，但沒有明顯的方法修改這個方案使之對希爾伯特曲線也適用。在網上搜尋，你很可能遇到很多關於希爾伯特曲線是怎樣畫出來的描述，但很少有關於找到任意點（在曲線上）位置的。為了搞定它，我們需要更仔細看看希爾伯特曲線是怎麼遞迴構建的。

首先要注意到雖然大多數關於希爾伯特曲線的文獻都關注曲線是怎麼畫出來的，卻容易讓我們忽略曲線的本質屬性以及其重要性：曲線規定了平面上點的順序。如果我們用這順序來表達希爾伯特曲線，畫曲線就不值一提了：僅僅是把點連起來。忘記怎麼把子曲線連起來吧，把注意力集中在怎麼遞迴地列舉點上。

圖6    希爾伯特曲線規定了二維平面上的點的順序

在根這一層，列舉點很簡單：選定一個方向和一個起始點，環繞四個象限，用0到3給他們編號。當我們要確定訪問子象限的順序同時維護總體的鄰接屬性，困難就來了。通過檢查我們發現，子象限的曲線是原曲線的簡單變換，而且只有四種變換。自然地，這個結論也適用於子子象限，等等。對於一個給定的象限，我們在其中畫出的曲線是由象限所在大的方形的曲線以及該象限的位置決定的。只需要費一點力，我們就能構建出如下概況所有情況的表。

圖7

假設我們想用這個表來確定某個點在第三層希爾伯特曲線上的位置。在這個例子中，假設點的座標是(5,2)。（譯者注：請參照圖8）從上圖的第一個方形開始，找到你的點所在的象限。在這個例子中，是在右上方的象限。那麼點在希爾伯特曲線上的位置的第一部分是3（二進位制是11）。接著我們進入象限3裡面的方塊，在這個例子中，它是（圖7中的）第二個方塊。重複剛才的過程：我們的點落在哪個子象限？這次是左下角，意味著位置的下一部分是1（二進位制01），我們將進入的小方塊又是第二個。最後一次重複這個過程，發現點落在右上角的子子象限，因此位置的最後部分是3（二進位制11）。把這些位置連線起來，我們得到點在曲線上的位置是二進位制的110111，或者十進位制的55。

圖8  三階希爾伯特曲線

讓我們更系統一些，寫出從x, y座標到希爾伯特曲線位置轉換的方法。首先，我們要以計算機看得懂的形式表達圖7：

上面的程式碼中，每個hilbert_map的元素對應圖7四個方形中的一個。為了容易區分，我用一個字母來標識每個方塊：’a’是第一個方塊，’b’是第二個，等等。每個方塊的值是個字典，將(子)象限的x, y座標對映到曲線上的位置（元組值的第一部分）以及下一個用到的方塊（元組值的第二部分）。下面的程式碼展示了怎麼用這個來將x, y座標轉換成希爾伯特曲線上的位置：

函式的輸入是為整數的x, y座標和曲線的階。一階曲線填充2×2的格子，二階曲線填充4×4的格子，等等。我們的x, y座標應該先標準化到0到2^order-1的區間。這個函式從最高位開始，逐步處理x, y座標的每個位元位。在每個階段中，通過測試對應的位元位，可以確定座標處於哪個（子）象限，還可以從我們之前定義的hilbert_map中取得在曲線上的位置以及下一個要用的方塊。在這階段取得的位置，加入到目前總的位置的最低兩位。在下一次迴圈的開頭，總的位置左移兩位以便給下一個位置騰出地方。

讓我們執行一下之前的例子來檢驗一下函式寫對了沒有：

對了！為了進一步測試，我們可以用這個函式生成一條希爾伯特曲線的有序點的完整列表，然後用電子製表軟體把它們畫出來看我們是否真的得到了一條希爾伯特曲線。在Python互動直譯器中輸入如下程式碼：

將輸出的文字貼上到檔案中，儲存為hilbert.csv，用你最喜歡的電子製表軟體開啟，將資料畫成一個散點圖。結果當然是一條漂亮的希爾伯特曲線！

將hilbert_map做簡單的反轉就能實現point_to_hilbert的逆向功能（將希爾伯特曲線上的位置轉換為x, y座標），把這個留給讀者作為練習吧。

結論

空間索引，從四叉樹到geohash到希爾伯特曲線，到這就結束了。最後一點說明：如果你將一條希爾伯特曲線上的x, y座標的有序序列寫成二進位制形式，對於順序你注意到什麼有趣的東西嗎？你想到了什麼？

結束前的一點警告：我在這裡描述的全部索引方法都只適用於索引點。如果你想索引線、折線或者多邊形，這些方法可能就不管用了。據我所知，已知的唯一能有效索引形體的演算法是R-tree，這是一種完全不同且更復雜的方法。

打賞支援我翻譯更多好文章，謝謝！

打賞譯者