作者|一縷殤流化隱半邊冰霜
關於鄰居的定義，相鄰即為鄰居，那麼鄰居分為2種，邊相鄰和點相鄰。邊相鄰的有4個方向，上下左右。點相鄰的也有4個方向，即4個頂點相鄰的。
如上圖，綠色的區域是一顆四叉樹表示的範圍，四叉樹上面有一個點，圖中黃色區域標明的點。現在想求四叉樹上黃色的點的希爾伯特曲線鄰居。圖中黑色的線就是一顆穿過四叉樹的希爾伯特曲線。希爾伯特曲線的起點0在左上角的方格中，終點63在右上角的方格中。
紅色的四個格子是黃色格子邊相鄰鄰居，藍色的四個格子是黃色格子的頂點相鄰的鄰居，所以黃色格子的鄰居為8個格子，分別表示的點是8，9，54，11，53，30，31，32 。可以看出來這些鄰居在表示的點上面並不是相鄰的。
那麼怎麼求四叉樹上任意一點的希爾伯特曲線鄰居呢？
邊鄰居最直接的想法就是 先拿到中心點的坐標 (i，j) ，然後通過坐標系的關係，拿到與它邊相鄰的 Cell 的坐標 (i + 1，j) ， (i - 1，j) ， (i，j - 1) ， (i，j + 1) 。
實際做法也是如此。不過這裡涉及到需要轉換的地方。這裡需要把希爾伯特曲線上的點轉換成坐標以後才能按照上面的思路來計算邊鄰居。
關於 CellID 的生成與數據結構，見筆者這篇《Google S2 中的 CellID 是如何生成的 ？》
按照上述的思路，實現出來的代碼如下：
func (ci CellID) EdgeNeighbors() [4]CellID {
level := ci.Level()
size := sizeIJ(level)
f, i, j, _ := ci.faceIJOrientation()
return [4]CellID{
cellIDFromFaceIJWrap(f, i, j-size).Parent(level),
cellIDFromFaceIJWrap(f, i+size, j).Parent(level),
cellIDFromFaceIJWrap(f, i, j+size).Parent(level),
cellIDFromFaceIJWrap(f, i-size, j).Parent(level),
}
}
邊按照，下邊，右邊，上邊，左邊，逆時針的方向依次編號0，1，2，3 。
接下來具體分析一下裡面的實現。
func sizeIJ(level int) int {
return 1 << uint(maxLevel-level)
}
sizeIJ 保存的是當前 Level 的格子邊長大小。這個大小是相對於 Level 30 來說的。比如 level = 29，那麼它的 sizeIJ 就是2，代表 Level 29 的一個格子邊長是由2個 Level 30 的格子組成的，那麼也就是22=4個小格子組成的。如果是 level = 28，那麼邊長就是4，由16個小格子組成。其他都以此類推。
func (ci CellID) faceIJOrientation() (f, i, j, orientation int) {
f = ci.Face()
orientation = f & swapMask
nbits := maxLevel - 7*lookupBits // first iteration
for k := 7; k >= 0; k-- {
orientation += (int(uint64(ci)>>uint64(k*2*lookupBits+1)) & ((1 << uint((2 * nbits))) - 1)) << 2
orientation = lookupIJ[orientation]
i += (orientation >> (lookupBits + 2)) << uint(k*lookupBits)
j += ((orientation >> 2) & ((1 << lookupBits) - 1)) << uint(k*lookupBits)
orientation &= (swapMask | invertMask)
nbits = lookupBits // following iterations
}
// 下面這個判斷詳細解釋
if ci.lsb()&0x1111111111111110 != 0 {
orientation ^= swapMask
}
return
}
這個方法就是把 CellID 再分解回原來的 i 和 j。這裡具體的過程在筆者這篇《Google S2 中的 CellID 是如何生成的 ？》裡面的 cellIDFromFaceIJ 方法裡面有詳細的敘述，這裡就不再贅述了。cellIDFromFaceIJ 方法和 faceIJOrientation 方法是互為逆方法。cellIDFromFaceIJ 是把 face，i，j 這個當入參傳進去，返回值是 CellID，faceIJOrientation 是把 CellID 分解成 face，i，j，orientation。faceIJOrientation 比 cellIDFromFaceIJ 分解出來多一個 orientation。
這裡需要重點解釋的是 orientation 怎麼計算出來的。
我們知道 CellID 的數據結構是 3位 face + 60位 position + 1位標誌位。那麼對於 Level - n 的非葉子節點，3位 face 之後，一定是有 2 * n 位二進位位，然後緊接著 2*(maxLevel - n) + 1 位以1開頭的，末尾都是0的二進位位。maxLevel = 30 。
例如 Level - 16，中間一定是有32位二進位位，然後緊接著 2*(30 - 16) + 1 = 29位。這29位是首位為1，末尾為0組成的。3 + 32 + 29 = 64 位。64位 CellID 就這樣組成的。
當 n = 30，3 + 60 + 1 = 64，所以末尾的1並沒有起任何作用。當 n = 29，3 + 58 + 3 = 64，於是末尾一定是 100 組成的。10對方向並不起任何作用，最後多的一個0也對方向不起任何作用。關鍵就是看10和0之間有多少個00 。當 n = 28，3 + 56 + 5 = 64，末尾5位是 10000，在10和0之間有一個「00」。「00」是會對方向產生影響，初始的方向應該再異或 01 才能得到。
關於 「00」 會對原始的方向產生影響，這點其實比較好理解。CellID 從最先開始的方向進行四分，每次四分都將帶來一次方向的變換。直到變換到最後一個4個小格子的時候，方向就不會變化了，因為在4個小格子之間就可以唯一確定是哪一個 Cell 被選中。所以這也是上面看到了， Level - 30 和 Level - 29 的方向是不變的，除此以外的 Level 是需要再異或一次 01 ，變換以後得到原始的 orientation。
最後進行轉換，具體代碼實現如下：
func cellIDFromFaceIJWrap(f, i, j int) CellID {
// 1.
i = clamp(i, -1, maxSize)
j = clamp(j, -1, maxSize)
// 2.
const scale = 1.0 / maxSize
limit := math.Nextafter(1, 2)
u := math.Max(-limit, math.Min(limit, scale*float64((i<<1)+1-maxSize)))
v := math.Max(-limit, math.Min(limit, scale*float64((j<<1)+1-maxSize)))
// 3.
f, u, v = xyzToFaceUV(faceUVToXYZ(f, u, v))
return cellIDFromFaceIJ(f, stToIJ(0.5*(u+1)), stToIJ(0.5*(v+1)))
}
轉換過程總共分為三步。第一步先處理 i，j 邊界的問題。第二步，將 i，j 轉換成 u，v 。第三步，u，v 轉 xyz，再轉回 u，v，最後轉回 CellID 。
第一步：
func clamp(x, min, max int) int {
if x < min {
return min
}
if x > max {
return max
}
return x
}
clamp 函數就是用來限定 i ， j 的範圍的。i，j 的範圍始終限定在 [-1，maxSize] 之間。
第二步：
最簡單的想法是將（i，j）坐標轉換為（x，y，z）（這個點不在邊界上），然後調用 xyzToFaceUV 方法投影到對應的 face 上。
我們知道在生成 CellID 的時候，stToUV 的時候，用的是一個二次變換：
func stToUV(s float64) float64 {
if s >= 0.5 {
return (1 / 3.) * (4*s*s - 1)
}
return (1 / 3.) * (1 - 4*(1-s)*(1-s))
}
但是此處，我們用的變換就簡單一點，用的是線性變換。
u = 2 * s - 1
v = 2 * t - 1
u，v 的取值範圍都被限定在 [-1，1] 之間。具體代碼實現：
const scale = 1.0 / maxSize
limit := math.Nextafter(1, 2)
u := math.Max(-limit, math.Min(limit, scale*float64((i<<1)+1-maxSize)))
v := math.Max(-limit, math.Min(limit, scale*float64((j<<1)+1-maxSize)))
第三步：找到葉子節點，把 u，v 轉成 對應 Level 的 CellID。
f, u, v = xyzToFaceUV(faceUVToXYZ(f, u, v))
return cellIDFromFaceIJ(f, stToIJ(0.5*(u+1)), stToIJ(0.5*(v+1)))
這樣就求得了一個 CellID 。
由於邊有4條邊，所以邊鄰居有4個。
return [4]CellID{
cellIDFromFaceIJWrap(f, i, j-size).Parent(level),
cellIDFromFaceIJWrap(f, i+size, j).Parent(level),
cellIDFromFaceIJWrap(f, i, j+size).Parent(level),
cellIDFromFaceIJWrap(f, i-size, j).Parent(level),
}
上面數組裡面分別會裝入當前 CellID 的下邊鄰居，右邊鄰居，上邊鄰居，左邊鄰居。
如果在地圖上顯示出來的話，就是下圖的這樣子。
中間方格的 CellID = 3958610196388904960 , Level 10 。按照上面的方法求出來的邊鄰居，分別是：
3958603599319138304 // 下邊鄰居
3958607997365649408 // 右邊鄰居
3958612395412160512 // 上邊鄰居
3958599201272627200 // 左邊鄰居
在地圖上展示出來：
這裡的共頂點鄰居和文章開始講的頂點鄰居有點區別。並且下面還會有一些看似奇怪的例子，也是筆者在實際編碼中踩過的坑，分享一下。
這裡先說明一種特殊情況，即 Cell 正好在地球的外切立方體的8個頂點上。那麼這個點的頂點鄰居只有3個，而不是4個。因為這8個頂點每個點只有3個面與其連接，所以每個面上有且只有一個 Cell 是它們的頂點鄰居。除去這8個點以外的 Cell 的頂點鄰居都有4個！
j
|
| (0,1) (1,1)
| (0,0) (1,0)
|
> i
在上述的坐標軸中，i 軸方向如果為1，就落在4個象限的右邊一列上。如果 j 軸方向如果為，就落在4個象限的上面一行上。
假設 Cell Level 不等於 30，即末尾標誌位1後面還有0，那麼這種 Cell 轉換成 i，j 以後，i，j 的末尾就都是1 。
上面的結論可以證明的，因為在 faceIJOrientation 函數拆分 Cell 的時候，如果遇到了都是0的情況，比如 orientation = 11，Cell 末尾都是0，那麼取出末尾8位加上orientation，00000000 11，經過 lookupIJ 轉換以後得到 1111111111 ，於是 i = 1111，j = 1111 ，方向還是 11。Cell 末尾的00還是繼續循環上述的過程，於是 i，j 末尾全是1111 了。
所以我們只需要根據 i，j 判斷入參給的 Level 在哪個象限，就可以把共頂點的鄰居都找到。
假設入參給定的 Level 小，即 Cell 的面積大，那麼就需要判斷當前 Cell (函數調用者) 的共頂點是位於入參 Cell 的4個頂點的哪個頂點上。Cell 是一個矩形，有4個頂點。當前 Cell (函數調用者) 離哪個頂點近，就選那個頂點為公共頂點。再依次求出以公共頂點周圍的4個 Cell 即可。
假設入參給定的 Level 大，即 Cell 的面積小，那麼也需要判斷入參 Cell 的共頂點是位於當前 Cell (函數調用者)的4個頂點的哪個頂點上。Cell 是一個矩形，有4個頂點。入參 Cell 離哪個頂點近，就選那個頂點為公共頂點。再依次求出以公共頂點周圍的4個 Cell 即可。
由於需要判斷位於一個 Cell 的四等分的哪一個，所以需要判斷它的4個孩子的位置情況。即判斷 Level - 1 的孩子的相對位置情況。
halfSize := sizeIJ(level + 1)
size := halfSize << 1
f, i, j, _ := ci.faceIJOrientation()
var isame, jsame bool
var ioffset, joffset int
這裡需要拿到 halfSize ，halfSize 其實就是入參 Cell 的孩子的格子的 size 。
if i&halfSize != 0 {
// 位於後邊一列，所以偏移量要加上一個格子
ioffset = size
isame = (i + size) < maxSize
} else {
// 位於左邊一列，所以偏移量要減去一個格子
ioffset = -size
isame = (i - size) >= 0
}
這裡我們根據 halfSize 那一位是否為1來判斷距離矩形的4個頂點哪個頂點近。這裡還需要注意的是，如果 i + size 不能超過 maxSize，如果超過了，就不在同一個 face 上了。同理， i - size 也不能小於 0，小於0頁不在同一個 face 上了。
j 軸判斷原理和 i 完全一致。
if j&halfSize != 0 {
// 位於上邊一行，所以偏移量要加上一個格子
joffset = size
jsame = (j + size) < maxSize
} else {
// 位於下邊一行，所以偏移量要減去一個格子
joffset = -size
jsame = (j - size) >= 0
}
最後計算結果，先把入參的 Cell 先計算出來，然後在把它周圍2個軸上的 Cell 計算出來。
results := []CellID{
ci.Parent(level),
cellIDFromFaceIJSame(f, i+ioffset, j, isame).Parent(level),
cellIDFromFaceIJSame(f, i, j+joffset, jsame).Parent(level),
}
如果 i，j 都在同一個 face 上，那麼共頂點就肯定不是位於外切立方體的8個頂點上了。那麼就可以再把第四個共頂點的 Cell 計算出來。
if isame || jsame {
results = append(results, cellIDFromFaceIJSame(f, i+ioffset, j+joffset, isame && jsame).Parent(level))
}
綜上，完整的計算共頂點鄰居的代碼實現如下：
func (ci CellID) VertexNeighbors(level int) []CellID {
halfSize := sizeIJ(level + 1)
size := halfSize << 1
f, i, j, _ := ci.faceIJOrientation()
fmt.Printf("halfsize 原始的值 = %v-%b", halfSize, halfSize)
var isame, jsame bool
var ioffset, joffset int
if i&halfSize != 0 {
// 位於後邊一列，所以偏移量要加上一個格子
ioffset = size
isame = (i + size) < maxSize
} else {
// 位於左邊一列，所以偏移量要減去一個格子
ioffset = -size
isame = (i - size) >= 0
}
if j&halfSize != 0 {
// 位於上邊一行，所以偏移量要加上一個格子
joffset = size
jsame = (j + size) < maxSize
} else {
// 位於下邊一行，所以偏移量要減去一個格子
joffset = -size
jsame = (j - size) >= 0
}
results := []CellID{
ci.Parent(level),
cellIDFromFaceIJSame(f, i+ioffset, j, isame).Parent(level),
cellIDFromFaceIJSame(f, i, j+joffset, jsame).Parent(level),
}
if isame || jsame {
results = append(results, cellIDFromFaceIJSame(f, i+ioffset, j+joffset, isame && jsame).Parent(level))
}
return results
}
下面來舉幾個例子。
第一個例子是相同大小 Cell 。入參和調用者 Cell 都是相同 Level - 10 的。
VertexNeighbors := cellID.Parent(10).VertexNeighbors(10)
// 11011011101111110011110000000000000000000000000000000000000000
3958610196388904960 // 右上角
3958599201272627200 // 左上角
3958603599319138304 // 右下角
3958601400295882752 // 左下角
第二個例子是不是大小的 Cell 。調用者 Cell 是默認 Level - 30 的。
VertexNeighbors := cellID.VertexNeighbors(10)
// 11011011101111110011110000000000000000000000000000000000000000
3958610196388904960 // 右下角
3958599201272627200 // 左下角
3958612395412160512 // 右上角
3958623390528438272 // 左上角
上面兩個例子可以說明一個問題，同樣是調用 VertexNeighbors(10) 得到的 Cell 都是 Level - 10 的，但是方向和位置是不同的。本質在它們共的頂點是不同的，所以生成出來的4個Cell生成方向也就不同。
在 C++ 的版本中，查找頂點鄰居有一個限制：
DCHECK_LT(level, this->level());
入參的 Level 必須嚴格的比要找的 Cell 的 Level 小才行。也就是說入參的 Cell 的格子面積大小要比 Cell 格子大小更小才行。但是在 Go 的版本實現中並沒有這個要求，入參或大或小都可以。
下面這個舉例，入參比 Cell 的 Level 小。（可以看到成都市已經小成一個點了）
VertexNeighbors := cellID.Parent(10).VertexNeighbors(5)
3957538172551823360 // 右下角
3955286372738138112 // 左下角
3959789972365508608 // 右上角
3962041772179193856 // 左上角
下面這個舉例，入參比 Cell 的 Level 大。（可以看到 Level 15 的面積已經很小了）
VertexNeighbors := cellID.Parent(10).VertexNeighbors(15)
3958610197462646784 // 左下角
3958610195315163136 // 右下角
3958610929754570752 // 左上角
3958609463023239168 // 右上角
最後回來文章開頭問的那個問題中。如何在四叉樹上如何求希爾伯特曲線的鄰居 ？經過前文的一些鋪墊，再來看這個問題，也許讀者心裡已經明白該怎麼做了。
查找全鄰居有一個要求，就是入參的 Level 的面積必須要比調用者 Cell 的小或者相等。即入參 Level 值不能比調用者的 Level 小。因為一旦小了以後，鄰居的 Cell 的面積變得巨大，很可能一個鄰居的 Cell 裡面就裝滿了原來 Cell 的所有鄰居，那這樣的查找並沒有任何意義。
舉個例子，如果入參的 Level 比調用者 Cell 的 Level 小。那麼查找它的全鄰居的時候，查出來會出現如下的情況：
AllNeighbors := cellID.Parent(10).AllNeighbors(5)
這個時候是可以查找到全鄰居的，但是可能會出現重疊 Cell 的情況，為何會出現這樣的現象，下面再分析。
如果入參和調用者 Cell 的 Level 是相同的，那麼查找到的全鄰居就是文章開頭說到的問題了。理想狀態如下：
具體實現如下：
func (ci CellID) AllNeighbors(level int) []CellID {
var neighbors []CellID
face, i, j, _ := ci.faceIJOrientation()
// 尋找最左下角的葉子節點的坐標。我們需要規範 i，j 的坐標。因為入參 Level 有可能比調用者 Cell 的 Level 要大。
size := sizeIJ(ci.Level())
i &= -size
j &= -size
nbrSize := sizeIJ(level)
for k := -nbrSize; ; k += nbrSize {
var sameFace bool
if k < 0 {
sameFace = (j+k >= 0)
} else if k >= size {
sameFace = (j+k < maxSize)
} else {
sameFace = true
// 上邊鄰居 和 下邊鄰居
neighbors = append(neighbors, cellIDFromFaceIJSame(face, i+k, j-nbrSize,
j-size >= 0).Parent(level))
neighbors = append(neighbors, cellIDFromFaceIJSame(face, i+k, j+size,
j+size < maxSize).Parent(level))
}
// 左邊鄰居，右邊鄰居，以及2個對角線上的頂點鄰居
neighbors = append(neighbors, cellIDFromFaceIJSame(face, i-nbrSize, j+k,
sameFace && i-size >= 0).Parent(level))
neighbors = append(neighbors, cellIDFromFaceIJSame(face, i+size, j+k,
sameFace && i+size < maxSize).Parent(level))
// 這裡的判斷條件有2個用途，一是防止32-bit溢出，二是循環的退出條件，大於size以後也就不用再找了
if k >= size {
break
}
}
return neighbors
}
上述代碼簡單的思路用注釋寫了。需要講解的部分現在來講解。
首先需要理解的是 nbrSize 和 size 的關係。為何會有 nbrSize ？ 因為入參 Level 是可以和調用者 Cell 的 Level 不一樣的。入參的 Level 代表的 Cell 可大可小也可能相等。最終結果是以 nbrSize 格子大小來表示的，所以循環中需要用 nbrSize 來控制格子的大小。而 size 只是原來調用者 Cell 的格子大小。
循環中 K 的變化，當 K = -nbrSize 的時候，這個時候循環只會計算左邊和右邊的鄰居。對角線上的頂點鄰居其實也是左邊鄰居和右邊鄰居的一種特殊情況。接下來 K = 0，就會開始計算上邊鄰居和下邊鄰居了。K 不斷增加，直到最後 K >= size ，最後一次循環內，會先計算一次左邊和右邊鄰居，再 break 退出。
調用者的 Cell 在中間位置，所以想要跳過這個 Cell 到達另外一邊(上下，或者左右)，那麼就需要跳過 size 的大小。具體代碼實現是 i + size 和 j + size 。
先看左右鄰居的循環掃描方式。
左鄰居是 i - nbrSize，j + k，k 在循環。這表示的就是左鄰居的生成方式。它生成了左鄰居一列。從左下角開始生成，一直往上生成到左上角。
右鄰居是 i + size，j + k，k 在循環。這表示的就是右鄰居的生成方式。它生成了右鄰居一列。從右下角開始生成，一直往上生成到右上角。
再看上下鄰居的循環掃描方式。
下鄰居是 i + k，j - nbrSize，k 在循環。這表示的就是下鄰居的生成方式。它生成了下鄰居一行。從下鄰居最左邊開始生成，一直往上生成到下鄰居最右邊。
上鄰居是 i + k，j + size，k 在循環。這表示的就是上鄰居的生成方式。它生成了上鄰居一行。從上鄰居最左邊開始生成，一直往上生成到上鄰居最右邊。
舉例：
中間 Cell 的周圍的全鄰居是上圖的 8 的相同 Level 的 Cell。
生成順序用需要標識出來了。1，2，5，6，7，8 都是左右鄰居生成出來的。3，4 是上下鄰居生成出來的。
上面這個例子是都是 Level - 10 的 Cell 生成出來的。全鄰居正好是8個。
AllNeighbors := cellID.Parent(10).AllNeighbors(10)
3958601400295882752,
3958605798342393856,
3958603599319138304,
3958612395412160512,
3958599201272627200,
3958607997365649408,
3958623390528438272,
3958614594435416064
再舉一個 Level 比調用者 Cell 的 Level 大的例子。
AllNeighbors := cellID.Parent(10).AllNeighbors(11)
3958600575662161920,
3958606622976114688,
3958603324441231360,
3958611570778439680,
3958600025906348032,
3958607172731928576,
3958603874197045248,
3958613220045881344,
3958599476150534144,
3958608821999370240,
3958623115650531328,
3958613769801695232
它的全鄰居生成順序如下：
1，2，5，6，9，10，11，12 都是左右鄰居，3，4，7，8 是上下鄰居。我們可以看到左右鄰居是從下往上生成的。上下鄰居是從左往右生成的。
如果 Level 更大，比如 Level - 15 ，就會生成更多的鄰居：
現在再解釋一下如果入參 Level 比調用者 Cell 的 Level 小的情況。
舉例，入參 Level = 9 。
AllNeighbors := cellID.Parent(10).AllNeighbors(9)
3958589305667977216,
3958580509574955008,
3958580509574955008,
3958615693947043840,
3958598101760999424,
3958606897854021632,
3958624490040066048,
3958615693947043840
生成的全鄰居如下：
可以看到本來有8個鄰居的，現在只有6個了。其實生成出來的還是8個，只不過有2個重複了。重複的見圖中深紅色的兩個 Cell。
為何會重疊？
中間調用者的 Level - 10 的 Cell 先畫出來。
因為是 Level - 9 的，所以它是中間那個 Cell 的四分之一。
我們把 Level - 10 的兩個上鄰居也畫出來。
可以看到上鄰居 Up 和頂點鄰居 up-right 都位於同一個 Level - 9 的 Cell 內了。所以上鄰居和右上角的頂點鄰居就都是同一個 Level - 9 的 Cell 。所以重疊了。同理，下鄰居和右下的頂點鄰居也重疊了。所以就會出現2個 Cell 重疊的情況。
而且中間也沒有空出調用者 Cell 的位置。因為 i + size 以後，範圍還在同一個 Level - 9 的 Cell 內。
如果 Level 更小，重疊情況又會發生變化。比如 Level - 5 。
AllNeighbors := cellID.Parent(10).AllNeighbors(5)
3953034572924452864,
3946279173483397120,
3946279173483397120,
3957538172551823360,
3955286372738138112,
3957538172551823360,
3962041772179193856,
3959789972365508608
畫在地圖上就是
重疊的位置也發生了變化。
至此，查找鄰居相關的算法都介紹完了。
空間搜索系列文章：
GitHub Repo：Halfrost-Field
Follow: halfrost · GitHub
Source: https://halfrost.com/go_s2_Hilbert_neighbor/