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…