Google S2 中的四叉樹求 LCA 最近公共祖先

一. 尋找父親節點和孩子節點

首先需要回顧一下希爾伯特曲線的生成方式，具體程式碼見筆者上篇文章的分析，在這個分析中，有4個方向比較重要，接下來的分析需要，所以把這4個方向的圖搬過來。

在舉例之前還需要說明一點，有些網站提供的二進位制轉換，並沒有標明有符號還是無符號的轉換，這樣就會導致使用者的一些誤解。筆者開始並沒有發現這個問題，導致掉入了這個坑，好一會才轉過彎來。筆者在網上查詢了很多線上轉換計算器的工具，都發現了這個問題。比如常見的線上進位制轉換http://tool.oschina.net/hexconvert，隨便找兩個64位的二進位制數，有符號的和無符號的分別轉換成十進位制，或者反過來轉換，你會驚喜的發現，兩次結果居然相同！例如你輸入 3932700003016900608 和 3932700003016900600，你會發現轉換成二進位制以後結果都是 11011010010011110000011101000100000000000000000000000000000000。但是很明顯這兩個數不同。

假如 3932700003016900608 是無符號的，3932700003016900600 是有符號的，正確的結果應該如下：

// 3932700003016900608
11011010010011110000011101000100000000000000000000000000000000

// 3932700003016900600
11011010010011110000011101000011111111111111111111111111111000複製程式碼

差距明顯很大。這種例子其實還有很多，隨便再舉出幾組：無符號的 3932700011606835200 和有符號的 3932700011606835000；無符號的 3932700020196769792 和有符號的 3932700020196770000；無符號的 3932700028786704384 和有符號的 3932700028786704400……可以舉的例子很多，這裡就不再舉了。

利用網上的這個工具，十進位制轉二進位制是無符號的轉換，二進位制轉十進位制就會變成有符號的轉換了。而 Google S2 預設是無符號的 CellID，所以用有符號的 CellID 會出現錯誤。所以轉換中需要注意一下。筆者之前沒有發現這一點的時候出現了一些問題，後來突然發現了這一點，茅塞頓開。

好了，進入正題。接下來直接看一個例子，筆者用例子來說明每個 Cell 之間的關係，以及如何查詢父親節點的。

假設有上圖這4個連在一起的 Cell。先根據經緯度把 CellID 計算出來。

對應的4個 CellID 分別是：

由於前兩位都是0，所以其實可以省略，但是為了讀者看的更清晰，筆者還是補全了64位，前面2個0還是補上

// 3932700003016900608 右上
0011011010010011110000011101000100000000000000000000000000000000      

// 3932700011606835200 左上
0011011010010011110000011101001100000000000000000000000000000000 

// 3932700020196769792 左下
0011011010010011110000011101010100000000000000000000000000000000

// 3932700028786704384 右下
0011011010010011110000011101011100000000000000000000000000000000複製程式碼

在前篇文章裡面我們也分析了 Cell 64位的結構，這裡是4個 Level 14的 Cell，所以末尾有 64 - 3 - 1 - 14 * 2 = 32 個 0 。從末尾往前的第33位是一個1，第34位，第35位是我們重點需要關注的。可以看到分別是00，01，10，11 。正好是連續的4個二進位制。

根據這個順序，我們可以匹配到當前這4個 Level 14 的 Cell 對應的順序是上圖圖中的圖0 。只不過當前方向旋轉了45°左右。

右上的 CellID 是 3932700003016900608 ，從右往左的第34位和第35位是00，從右上這個 Cell 想找到希爾伯特曲線的下一個 Cell，由於它目前方向是圖0的方向，所以右上的 Cell 的下一個 Cell 是左上那個 Cell，那麼第34位和第35位就應該是01，變換成01以後，就3932700011606835200，這也就是對應的 CellID。右數第34位增加了一個1，對應十進位制就增加了 2^33^ = 8589934592，算一下兩個 CellID 對應十進位制的差值，也正好是這個數。目前一切都是正確的。

同理可得，左上的 Cell 的下一個 Cell 就是左下的 Cell，也是相同的第34位和第35位上變成10，對應十進位制增加 2^33^ = 8589934592，得到左下的 CellID 為 3932700020196769792。繼續同理可以得到最後一個 CellID，右下的 CellID，為 3932700028786704384。

看到這裡，讀者應該對查詢同 Level 的兄弟節點的方法清清楚楚了。可能有人有疑問了，要查詢父親節點和孩子節點和兄弟有啥關係？他們之間的聯絡就在這一章節開頭說的4個方向圖上面。

回顧一下希爾伯特曲線的生成方式，在遞迴生成希爾伯特曲線的時候，儲存了一個 posToIJ 陣列，這個陣列裡面記錄著4個方向。希爾伯特曲線生成的形式和這4個方向是密切相關的。如果忘記了這部分，還請回看之前筆者的文章分析。

所以同一級的 Level 想查詢孩子節點，首先要找到在這一級中，當前 Cell 所處的位置以及當前 Cell 所在的4個方向中的哪一個方向。這個方向決定了要查詢孩子位於下一級或者下幾級的哪個位置。因為希爾伯特曲線相當於是一個顆四叉樹，每個根節點有4個孩子，雖然按層可以很輕鬆的遍歷到孩子所在的層級，但是同一個根節點的孩子有4個，究竟要選哪一個就需要父親節點的方向一級級的來判斷了。

舉個例子來說明這個問題：

func testS2() {

    latlng := s2.LatLngFromDegrees(29.323773, 107.727194)
    cellID := s2.CellIDFromLatLng(latlng)
    cell := s2.CellFromCellID(cellID) //9279882742634381312

    // cell.Level()
    fmt.Println("latlng = ", latlng)
    fmt.Println("cell level = ", cellID.Level())
    fmt.Printf("cell = %d %b\n", cellID, cellID)
    smallCell := s2.CellFromCellID(cellID.Parent(13))
    fmt.Printf("smallCell level = %d\n", smallCell.Level())
    fmt.Printf("smallCell id = %d / cellID = %d /smallCell(14).Parent = %d (level = %d)/smallCell(15).Parent = %d (level = %d)/cellID(13).Parent = %d (level = %d)/cellID(14).Parent = %d (level = %d)/cellID(15).Parent = %d (level = %d)/ %b \n", smallCell.ID(), cellID, smallCell.ID().Parent(14), (smallCell.ID().Parent(14).Level()), smallCell.ID().Parent(15), (smallCell.ID().Parent(15).Level()), cellID.Parent(13), (cellID.Parent(13)).Level(), cellID.Parent(14), (cellID.Parent(14)).Level(), cellID.Parent(15), (cellID.Parent(15)).Level(), smallCell.ID())

}複製程式碼

讀者考慮一下上述的程式會輸出什麼呢？或者這樣問，同一個節點，先找它的 Level 13 的父親節點，再通過 Level 13 的這個節點再找它的 Level 15 的父親節點得到的 Level 15 的節點，和，直接找它的 Level 15 的父親節點，最終結果一樣麼？

當然，這裡先找 Level 13，再找 Level 15，得到的結果不是 Level 28，這裡不是相加的關係，結果還是 Level 15 的父親節點。所以上面兩種方式得到的結果都是 Level 15的。那麼兩種做法得到的結果是一樣的麼？讀者可以先猜一猜。

實際得到的結果如下：

latlng =  [29.3237730, 107.7271940]
cell level =  30
cell = 3932700032807325499 11011010010011110000011101011111101111101001011100111100111011
smallCell level = 13
smallCell id = 3932700015901802496 / cellID = 3932700032807325499 /smallCell(14).Parent = 3932700020196769792 (level = 14)/smallCell(15).Parent = 3932700016975544320 (level = 15)/cellID(13).Parent = 3932700015901802496 (level = 13)/cellID(14).Parent = 3932700028786704384 (level = 14)/cellID(15).Parent = 3932700032007929856 (level = 15)/ 11011010010011110000011101010000000000000000000000000000000000複製程式碼

可以看到，兩種做法得到的結果是不同的。但是究竟哪裡不同呢？直接看 CellID 不夠直觀，那把它們倆畫出來看看。

可以看到兩者雖然 Level 是相同的，但是位置是不同的，為何會這樣呢？原因就是之前說的，四叉樹4個孩子，究竟選擇哪一個，是由於父親節點所在方向決定的。

1. 查詢孩子節點

還是繼續按照上面程式舉的例子，看看如何查詢孩子節點的。

我們把 Level 13 - Level 17 的節點都列印出來。

smallCell id = 3932700015901802496 (level = 13)/
smallCell(14).Parent = 3932700020196769792 (level = 14)/ 
smallCell(15).Parent = 3932700016975544320 (level = 15)/
smallCell(16).Parent = 3932700016170237952 (level = 16)/
smallCell(17).Parent = 3932700015968911360 (level = 17)/複製程式碼

畫在圖上，

從 Level 13 是如何變換到 Level 14 的呢？我們知道當前選擇的是圖0的方向。那麼當前 Level 14是圖0 中的2號的位置。

所以從右往左數 64 - 3 - 1 - 13 * 2 = 34位和35位，應該分別填上01，從前往後就是10，對應的就是上圖中的2的位置，並且末尾的那個標誌位1往後再挪2位。

11011010010011110000011101010000000000000000000000000000000000   13
11011010010011110000011101010100000000000000000000000000000000   14複製程式碼

即可從 Level 13 變換到 Level 14 。這就是從父親節點到孩子節點的變換方法。

同理在看看 Level 15，Level 16，Level 17 的變換方法。

前面說過，檢視孩子節點的時候需要知道當前節點的其他3個兄弟節點的方向。

根據下圖，Level 14 對應的是圖0，並且當前選擇了2號位置，從下圖中可以看到圖0中的2號位置的下一級的圖是“U”形的，說明還是圖0的樣子。

所以可以知道當前 Level 14 所處的方向依舊是圖0 。按照方向標識在圖上，如下圖。

所以如果還是選擇上圖中0號的位置，那麼 Level 15 的從右往左數 64 - 3 - 1 - 14 * 2 = 32位和第33位上填入00 。

11011010010011110000011101010100000000000000000000000000000000   14 
11011010010011110000011101010001000000000000000000000000000000   15複製程式碼

由於選擇了圖0的0號位置，所以下一級的方向對應的是圖1 。(注意整個圖的方向是向左旋轉了90°) 。

圖1中繼續選擇0號的位置，所以 Level 16 的從右往左數 64 - 3 - 1 - 15 * 2 = 30位和31位填上00 。那麼就可以得到 Level 16 。

11011010010011110000011101010001000000000000000000000000000000   15
11011010010011110000011101010000010000000000000000000000000000   16複製程式碼

由於選擇了圖1的0號位置，所以下一級的方向對應的是還是圖0。

圖0中繼續選擇0號的位置，所以 Level 17 的從右往左數 64 - 3 - 1 - 16 * 2 = 28位和第29位填上00 。那麼就可以得到 Level 17 。

11011010010011110000011101010000010000000000000000000000000000   16
11011010010011110000011101010000000100000000000000000000000000   17複製程式碼

同理，其他的孩子節點都可以按照這個方法推算得到。

從 Level 13 開始，由於 Level 13 對應的方向是圖0，當前選擇3號位置，就可以得到 Level 14，所以 Level 14 的末尾標誌位1前面的兩位是 11 。於是就可以從 Level 13 變換到 Level 14 。

11011010010011110000011101010000000000000000000000000000000000   13 3932700015901802496
11011010010011110000011101011100000000000000000000000000000000   14 3932700028786704384複製程式碼

由於圖0選擇了3號位置，那麼 Level 14 的方向就是圖3 。

Level 14 對應的方向是圖3，當前選擇3號位置，就可以得到 Level 15，所以 Level 15 的末尾標誌位1前面的兩位是 11 。於是就可以從 Level 14 變換到 Level 15 。

11011010010011110000011101011100000000000000000000000000000000   14 3932700028786704384
11011010010011110000011101011111000000000000000000000000000000   15 3932700032007929856複製程式碼

由於圖3選擇了3號位置，那麼 Level 15 的方向就是圖0。

Level 15 對應的方向是圖0，當前選擇0號位置，就可以得到 Level 16，所以 Level 16 的末尾標誌位1前面的兩位是 00 。於是就可以從 Level 15 變換到 Level 16 。

11011010010011110000011101011111000000000000000000000000000000   15 3932700032007929856
11011010010011110000011101011110010000000000000000000000000000   16 3932700031202623488複製程式碼

由於圖0選擇了0號位置，那麼 Level 16 的方向就是圖1。

Level 16 對應的方向是圖1，當前選擇1號位置，就可以得到 Level 17，所以 Level 17 的末尾標誌位1前面的兩位是 01 。於是就可以從 Level 16 變換到 Level 17 。

11011010010011110000011101011110010000000000000000000000000000   16 3932700031202623488
11011010010011110000011101011110001100000000000000000000000000   17 3932700031135514624複製程式碼

由於圖1選擇了1號位置，那麼 Level 18 的方向還是圖1。

到此讀者應該對查詢 CellID 孩子節點的流程瞭然於心了。在 Google S2 中，查詢孩子節點的具體實現程式碼如下。

func (ci CellID) Children() [4]CellID {
    var ch [4]CellID
    lsb := CellID(ci.lsb())
    ch[0] = ci - lsb + lsb>>2
    lsb >>= 1
    ch[1] = ch[0] + lsb
    ch[2] = ch[1] + lsb
    ch[3] = ch[2] + lsb
    return ch
}複製程式碼

現在在來看這段程式碼應該毫無壓力了吧。

這裡比較重要的位運算的操作就是 lsb 了。從名字上其實也可以知道它是做什麼的。

// lsb 返回最低有效位
func (ci CellID) lsb() uint64 { return uint64(ci) & -uint64(ci) }複製程式碼

這裡需要注意的一點就是負數的儲存方式是以原碼的補碼，即符號位不變，每位取反再加1 。

舉個例子，Level 16 的某個 CellID 如下：

11011010010011110000011101011110010000000000000000000000000000   16 3932700031202623488複製程式碼

對它進行 lsb 計算：

得到的結果就是最低有效位為1，其他每位都為0 。

ch[0] = ci - lsb + lsb>>2複製程式碼

這一行實際是把 Level 對應的下一級 Level 的末尾標誌位1移動到位。即往後挪2位。並且標誌位前面2位都為0，所以這步操作完成以後就是0號的孩子。

0號孩子找到以後接下來就很好辦了。lsb 往右移動一位以後，不斷的加上這個值，就可以得到剩下的4個孩子了。如下圖：

這樣就可以得到4個孩子，上面這一小段程式挺簡單的，比前面從地圖上解釋的更簡單，原因是因為沒有視覺化的4個孩子的相互位置關係，這個關係需要從當前所在的方向來決定。前面地圖上也一再強調每一級的方向位置關係也是為了視覺化展現在地圖上是符合希爾伯特曲線的相對位置。

2. 判斷是否是葉子節點

如果對 CellID 的資料結構很瞭解，這個判斷就很簡單了。

func (ci CellID) IsLeaf() bool { return uint64(ci)&1 != 0 }複製程式碼

由於 CellID 是64位的，末尾有一個1的標誌位，如果這個標誌位到了最後一位，那麼就肯定是葉子節點了，也就是 Level 30 的 Cell。

3. 查詢當前孩子位置關係

在前面講解查詢孩子節點的時候，由於是四叉樹，每個父親下面對應4個孩子，00，01，10，11，所以判斷4個孩子之間相對的位置關係只需要判斷這兩個二進位制位就可以了。

func (ci CellID) ChildPosition(level int) int {
    return int(uint64(ci)>>uint64(2*(maxLevel-level)+1)) & 3
}複製程式碼

上面這個函式入參是一個父親節點的 Level 等級，返回的是這個父親節點下面孩子節點的位置資訊。即是 00，01，10，11 中的一個。

4. 查詢父親節點

在 Google S2 中，由於預設生成出來的 Cell 就是 Level 30 的，也就是 Level 最低的，位於樹的最下層的葉子節點。所以生成 Level 比較低的 Cell 必須只能查詢父親節點。

由於前面講解了如何查詢孩子節點，查詢父親節點就是逆向的過程。

func lsbForLevel(level int) uint64 { return 1 << uint64(2*(maxLevel-level)) }複製程式碼

第一步就是先找到最右邊的標誌位，它決定了 Level 的值。

(uint64(ci) & -lsb)複製程式碼

第二步是保留住標誌位前面所有的二進位制位上的值。這裡對第一步的 lsb 的相反數進行按位與操作就可以實現。lsb 的相反數其實就是 lsb 低位往左的高位都為1 ，相當於一個掩碼。

最後一步將標誌位1放好就可以了。

func (ci CellID) Parent(level int) CellID {
    lsb := lsbForLevel(level)
    return CellID((uint64(ci) & -lsb) | lsb)
}複製程式碼

以上就是查詢父親節點的具體實現。

二. LCA 查詢最近公共祖先

關於 CellID 的計算，還有很關鍵的一部分就是查詢最近公共祖先的問題。問題背景：給定一棵四叉樹中任意兩個 Level 的 CellID ，如何查詢兩者的最近公共祖先。

由 CellID 的資料結構我們知道，想查詢兩個 Level 的最近公共祖先的問題可以轉化為從左往右查詢兩個二進位制串最長公共序列，最長的即是從根節點開始最遠的公共祖先，也就是最近公共祖先。

那麼現在問題就轉換成從左往右找到第一個不相同的二進位制位，或者從右往左找到最後一個不相同的二進位制位。

查詢過程中存在一個特殊情況，那就是要查詢公共祖先的兩個節點本身就在一個分支上，即其中一個 CellID 本來就是另外一個 CellID 的祖先，那麼他們倆的公共祖先就直接是 CellID 大的那個。

那麼到此就可以確定出接下來查詢的流程。

第一步，先對兩個 CellID 進行異或，找到不同的二進位制位分別在那些位置上。

    bits := uint64(ci ^ other)複製程式碼

第二步，判斷是否存在特殊情況：兩個存在祖先關係。

    if bits < ci.lsb() {
        bits = ci.lsb()
    }
    if bits < other.lsb() {
        bits = other.lsb()
    }複製程式碼

第三步，查詢左邊最高位不同的位置。

func findMSBSetNonZero64(bits uint64) int {
    val := []uint64{0x2, 0xC, 0xF0, 0xFF00, 0xFFFF0000, 0xFFFFFFFF00000000}
    shift := []uint64{1, 2, 4, 8, 16, 32}
    var msbPos uint64
    for i := 5; i >= 0; i-- {
        if bits&val[i] != 0 {
            bits >>= shift[i]
            msbPos |= shift[i]
        }
    }
    return int(msbPos)
}複製程式碼

由於 CellID 是 64 位的，所以兩者不同的位數可能的範圍是 [0，63]。分別準備6種掩碼，對應的分別是0x2, 0xC, 0xF0, 0xFF00, 0xFFFF0000, 0xFFFFFFFF00000000。如下圖。

查詢的過程就是利用的二分的思想。64位先查詢高位32位，如果對高32位的掩碼進行按位與運算以後結果不為0，那麼就說明高32位是存在不相同的位數的，那麼最終結果 msbPos 加上32位，並把數字右移32位。因為高32位存在不同的數，由於我們需要求最左邊的，所以低32位就可以直接捨去了，直接右移去掉。

同理，繼續二分，16位，8位，4位，2位，1位，這樣迴圈完，就一定能把最左邊的不同的位數找到，並且結果位即為 msbPos。

第四步，判斷 msbPos 的合法性，並輸出最終結果。

如果 msbPos 比60還要大，那麼就是非法值，直接返回 false 即可。

    msbPos := findMSBSetNonZero64(bits)
    if msbPos > 60 {
        return 0, false
    }
    return (60 - msbPos) >> 1, true複製程式碼

最終輸出的為最近公共祖先的 Level 值，所以 60 - msbPos 以後還需要再除以2 。

完整的演算法實現如下：

func (ci CellID) CommonAncestorLevel(other CellID) (level int, ok bool) {
    bits := uint64(ci ^ other)
    if bits < ci.lsb() {
        bits = ci.lsb()
    }
    if bits < other.lsb() {
        bits = other.lsb()
    }

    msbPos := findMSBSetNonZero64(bits)
    if msbPos > 60 {
        return 0, false
    }
    return (60 - msbPos) >> 1, true
}複製程式碼

舉個例子：

在上面的例子中，我們挑選不存在祖先關係的兩個 Level 的 CellID。

11011010010011110000011101010000000100000000000000000000000000   17
11011010010011110000011101011111000000000000000000000000000000   15複製程式碼

如果從這串二進位制裡面直接找最近公共祖先，一定可以發現，從左往右最長的公共二進位制串是：

1101101001001111000001110101複製程式碼

那麼他們倆的最近公共祖先就是：

11011010010011110000011101010000000000000000000000000000000000複製程式碼

對應的 Level 是13，CellID 是 3932700015901802496。

---

空間搜尋系列文章：

如何理解 n 維空間和 n 維時空
高效的多維空間點索引演算法 — Geohash 和 Google S2
Google S2 中的四叉樹求 LCA 最近公共祖先

GitHub Repo：Halfrost-Field

Follow: halfrost · GitHub

Source: halfrost.com/go_s2_lowes…