空間點像素索引（三）

1. 點與坐標軸點相互轉換

在 S2 算法中，默認劃分 Cell 的等級是30，也就是說把一個正方形劃分為 2^30 * 2^30個小的正方形。那么上一步的s，t映射到這個正方形上面來，對應該如何轉換呢？

s，t的值域是[0,1]，現在值域要擴大到[0,2³⁰-1]。

// stToIJ converts value in ST coordinates to a value in IJ
coordinates.

func stToIJ(s float64) int {

return clamp(int(math.Floor(maxSize*s)),
0, maxSize-1)

}

C ++ 的實現版本也一樣

inline int S2CellId::STtoIJ(double s) {

// Converting from floating-point to integers via
static_cast is very slow

// on Intel processors because it requires changing the rounding
mode.

// Rounding to the nearest integer using FastIntRound() is
much faster.

// 這里減去0.5是為了四舍五入

return max(0, min(kMaxSize - 1, MathUtil::FastIntRound(kMaxSize

• s - 0.5)));

}

到這一步，是h(face,s,t) -> H(face,i,j)。

坐標軸點與希爾伯特曲線Cell ID相互轉換

最后一步，如何把 i，j 和希爾伯特曲線上的點關聯起來呢？

const (

lookupBits = 4

swapMask
= 0x01

invertMask = 0x02

)

var (

ijToPos = [4][4]int{

{0, 1, 3, 2}, // canonical order

{0, 3, 1, 2}, // axes swapped

{2, 3, 1, 0}, // bits inverted

{2, 1, 3, 0}, // swapped &
inverted

}

posToIJ = [4][4]int{

{0, 1, 3, 2}, // canonical order: (0,0), (0,1), (1,1), (1,0)

{0, 2, 3, 1}, // axes swapped: (0,0), (1,0), (1,1), (0,1)

{3, 2, 0, 1}, // bits inverted: (1,1), (1,0), (0,0), (0,1)

{3, 1, 0, 2}, // swapped &
inverted: (1,1), (0,1), (0,0), (1,0)

}

posToOrientation = [4]int{swapMask, 0, 0, invertMask | swapMask}

lookupIJ [1 << (2*lookupBits + 2)]int

lookupPos [1 << (2*lookupBits + 2)]int

)

在變換之前，先來解釋一下定義的一些變量。posToIJ 代表的是一個矩陣，里面記錄了一些單元希爾伯特曲線的位置信息。把 posToIJ 數組里面的信息用圖表示出來，如下圖：

同理，把 ijToPos 數組里面的信息用圖表示出來，如下圖：

posToOrientation 數組里面裝了4個數字，分別是1,0,0,3。lookupIJ 和 lookupPos 分別是兩個容量為1024的數組。這里面分別對應的就是希爾伯特曲線 ID 轉換成坐標軸 IJ 的轉換表，和坐標軸 IJ 轉換成希爾伯特曲線 ID 的轉換表。

func init() {

initLookupCell(0, 0, 0, 0, 0, 0)

initLookupCell(0, 0, 0, swapMask, 0, swapMask)

initLookupCell(0, 0, 0, invertMask, 0, invertMask)

initLookupCell(0, 0, 0, swapMask|invertMask, 0, swapMask|invertMask)

}

這里是初始化的遞歸函數，在希爾伯特曲線的標準順序中可以看到是有4個格子，并且格子都有順序的，所以初始化要遍歷滿所有順序。入參的第4個參數，就是從0 - 3 。

// initLookupCell initializes the lookupIJ table at init
time.

func initLookupCell(level,
i, j, origOrientation, pos, orientation int) {

if level == lookupBits {

ij := (i << lookupBits) + j

lookupPos[(ij<<2)+origOrientation] = (pos <<
2) + orientation

lookupIJ[(pos<<2)+origOrientation] = (ij << 2) + orientation

return

}

level++

i <<= 1

j <<= 1

pos <<= 2

r := posToIJ[orientation]

initLookupCell(level, i+(r[0]>>1), j+(r[0]&1), origOrientation, pos, orientation^posToOrientation[0])

initLookupCell(level, i+(r[1]>>1), j+(r[1]&1), origOrientation, pos+1, orientation^posToOrientation[1])

initLookupCell(level, i+(r[2]>>1), j+(r[2]&1), origOrientation, pos+2, orientation^posToOrientation[2])

initLookupCell(level, i+(r[3]>>1), j+(r[3]&1), origOrientation, pos+3, orientation^posToOrientation[3])

}

上面這個函數是生成希爾伯特曲線的。我們可以看到有一處對pos << 2的操作，這里是把位置變換到第一個4個小格子中，所以位置乘以了4。由于初始設置的lookupBits = 4，所以i，j的變化范圍是從[0,15]，總共有1616=256個，然后i，j坐標是表示的4個格子，再細分，lookupBits = 4這種情況下能表示的點的個數就是2564=1024個。這也正好是 lookupIJ 和 lookupPos 的總容量。畫一個局部的圖，i，j從0-7變化。

上圖是一個4階希爾伯特曲線。初始化的實際過程就是初始化4階希爾伯特上的1024個點的坐標與坐標軸上的x，y軸的對應關系表。舉個例子，下表是i，j在遞歸過程中產生的中間過程。下表是 lookupPos 表計算過程。

取出一行詳細分析一下計算過程。假設當前(i,j)=(0,2)，ij 的計算過程是把 i 左移4位再加上 j，整體結果再左移2位。目的是為了留出2位的方向位置。ij的前4位是i，接著4位是j，最后2位是方向。這樣計算出ij的值就是8 。接著計算lookupPos[i j]的值。從上圖中可以看到(0,2)代表的單元格的4個數字是16，17，18，19 。計算到這一步，pos的值為4（pos是專門記錄生成格子到第幾個了，總共pos的值會循環0-255）。pos代表的是當前是第幾個格子(4個小格子組成)，當前是第4個，每個格子里面有4個小格子。所以44就可以偏移到當前格子的第一個數字，也就是16 。posToIJ 數組里面會記錄下當前格子的形狀。從這里我們從中取出 orientation 。看上圖，16，17，18，19對應的是 posToIJ 數組軸旋轉的情況，所以17是位于軸旋轉圖的數字1代表的格子中。這時 orientation = 1 。這樣 lookupPos[i j] 表示的數字就計算出來了，44+1=17 。這里就完成了i，j與希爾伯特曲線上數字的對應。那如何由希爾伯特曲線上的數字對應到實際的坐標呢？lookupIJ 數組里面記錄了反向的信息。lookupIJ 數組
和 lookupPos 數組存儲的信息正好是反向的。lookupIJ 數組
下表存的值是 lookupPos 數組
的下表。我們查 lookupIJ 數組
，lookupIJ[17]的值就是8，對應算出來(i,j)=(0,2)。這個時候的i，j還是大格子。還是需要借助 posToIJ 數組
里面描述的形狀信息。當前形狀是軸旋轉，之前也知道 orientation = 1，由于每個坐標里面有4個小格子，所以一個i，j代表的是2個小格子，所以需要乘以2，再加上形狀信息里面的方向，可以計算出實際的坐標 (0 * 2 + 1 , 2 * 2 + 0) = ( 1，4) 。至此，整個球面坐標的坐標映射就已經完成了。球面上的點S(lat,lng) -> f(x,y,z) -> g(face,u,v)
-> h(face,s,t) -> H(face,i,j) -> CellID。目前總共轉換了6步，球面經緯度坐標轉換成球面xyz坐標，再轉換成外切正方體投影面上的坐標，最后變換成修正后的坐標，再坐標系變換，映射到 [0,2³⁰-1]區間，最后一步就是把坐標系上的點都映射到希爾伯特曲線上。

S2 Cell ID 數據結構

最后需要來談談 S2 Cell ID 數據結構，這個數據結構直接關系到不同 Level 對應精度的問題。

上圖左圖中對應的是 Level 30 的情況，右圖對應的是 Level 24 的情況。(2的多少次方，角標對應的也就是 Level 的值)在 S2 中，每個 CellID 是由64位的組成的。可以用一個 uint64 存儲。開頭的3位表示正方體6個面中的一個，取值范圍[0,5]。3位可以表示0-7，但是6，7是無效值。64位的最后一位是1，這一位是特意留出來的。用來快速查找中間有多少位。從末尾最后一位向前查找，找到第一個不為0的位置，即找到第一個1。這一位的前一位到開頭的第4位（因為前3位被占用）都是可用數字。綠色格子有多少個就能表示劃分多少格。上圖左圖，綠色的有60個格子，于是可以表示[0,2³⁰ -1] * [0,2³⁰ -1]這么多個格子。上圖右圖中，綠色格子只有48個，那么就只能表示[0,2²⁴ -1]*[0,2²⁴ -1]這么多個格子。那么不同 level 可以代表的網格的面積究竟是多大呢？

由上一章我們知道，由于投影的原因，所以導致投影之后的面積依舊有大小差別。這里推算的公式比較復雜，就不證明了，具體的可以看文檔。

MinAreaMetric = Metric{2, 8 * math.Sqrt2 / 9}

AvgAreaMetric = Metric{2, 4 * math.Pi / 6}

MaxAreaMetric = Metric{2, 2.635799256963161491}

這就是最大最小面積和平均面積的倍數關系。(下圖單位是km²，平方公里)

level 0 就是正方體的六個面之一。地球表面積約等于510,100,000 km²。level 0 的面積就是地球表面積的六分之一。level 30 能表示的最小的面積0.48cm²，最大也就0.93cm² 。

4. 與Geohash對比
   

Geohash有12級，從5000km 到3.7cm。中間每一級的變化比較大。有時候可能選擇上一級會大很多，選擇下一級又會小一些。比如選擇字符串長度為4，它對應的 cell 寬度是39.1km，需求可能是50km，那么選擇字符串長度為5，對應的 cell 寬度就變成了156km，瞬間又大了3倍了。這種情況選擇多長的 Geohash 字符串就比較難選。選擇不好，每次判斷可能就還需要取出周圍的8個格子再次進行判斷。Geohash 需要 12 bytes 存儲。S2有30級，從 0.7cm2 到 85,000,000km2 。中間每一級的變化都比較平緩，接近于4次方的曲線。所以選擇精度不會出現Geohash 選擇困難的問題。S2 的存儲只需要一個uint64 即可存下。S2庫里面不僅僅有地理編碼，還有其他很多幾何計算相關的庫。地理編碼只是其中的一小部分。本文沒有介紹到的 S2 的實現還有很多很多，各種向量計算，面積計算，多邊形覆蓋，距離問題，球面球體上的問題，它都有實現。S2還能解決多邊形覆蓋的問題。比如給定一個城市，求一個多邊形剛剛好覆蓋住這個城市。

如上圖，生成的多邊形剛剛好覆蓋住下面藍色的區域。這里生成的多邊形可以有大有小。不管怎么樣，最終的結果也是剛剛覆蓋住目標物。

用相同的 Cell 也可以達到相同的目的，上圖就是用相同 Level 的 Cell 覆蓋了整個圣保羅城市。這些都是 Geohash 做不到的。多邊形覆蓋利用的是近似的算法，雖然不是嚴格意義上的最優解，但是實踐中效果特別好。額外值得說明的一點是，Google 文檔上強調了，這種多邊形覆蓋的算法雖然對搜索和預處理操作非常有用，但是“不可依賴”的。理由也是因為是近似算法，并不是唯一最優算法，所以得到的解會依據庫的不同版本而產生變化。

5. S2 Cell舉例
   

先來看看經緯度和 CellID 的轉換，以及矩形面積的計算。

latlng := s2.LatLngFromDegrees(31.232135, 121.41321700000003)

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\n”, cellID)

smallCell := s2.CellFromCellID(cellID.Parent(10))

fmt.Printf(“smallCell
level = %d\n”, smallCell.Level())

fmt.Printf(“smallCell
id = %b\n”, smallCell.ID())

fmt.Printf(“smallCell
ApproxArea = %v\n”,
smallCell.ApproxArea())

fmt.Printf(“smallCell
AverageArea = %v\n”,
smallCell.AverageArea())

fmt.Printf(“smallCell
ExactArea = %v\n”,
smallCell.ExactArea())

這里 Parent 方法參數可以直接指定返回改點的對應 level 的 CellID。上面那些方法打印出來的結果如下：

latlng = [31.2321350, 121.4132170]

cell level = 30

cell = 3869277663051577529

****Parent **** 10000000000000000000000000000000000000000

smallCell level = 10

smallCell id = 11010110110010011011110000000000000000000000000000000000000000

smallCell ApproxArea = 1.9611002454714756e-06

smallCell AverageArea = 1.997370817559429e-06

smallCell ExactArea = 1.9611009480261058e-06

再舉一個覆蓋多邊形的例子。我們先隨便創建一個區域。

rect = s2.RectFromLatLng(s2.LatLngFromDegrees(48.99, 1.852))

rect = rect.AddPoint(s2.LatLngFromDegrees(48.68, 2.75))

rc := &s2.RegionCoverer{MaxLevel: 20, MaxCells: 10, MinLevel: 2}

r := s2.Region(rect.CapBound())

covering := rc.Covering?

覆蓋參數設置成 level 2 - 20，最多的 Cell 的個數是10個。

接著我們把 Cell 至多改成20個。

最后再改成30個。

可以看到相同的 level 的范圍，cell 個數越多越精確目標范圍。這里是匹配矩形區域，匹配圓形區域也同理。

代碼就不貼了，與矩形類似。這種功能 Geohash 就做不到，需要自己手動實現了。最后舉一個多邊形匹配的例子。

func testLoop() {

ll1 := s2.LatLngFromDegrees(31.803269, 113.421145)

ll2 := s2.LatLngFromDegrees(31.461846, 113.695803)

ll3 := s2.LatLngFromDegrees(31.250756, 113.756228)

ll4 := s2.LatLngFromDegrees(30.902604, 113.997927)

ll5 := s2.LatLngFromDegrees(30.817726, 114.464846)

ll6 := s2.LatLngFromDegrees(30.850743, 114.76697)

ll7 := s2.LatLngFromDegrees(30.713884, 114.997683)

ll8 := s2.LatLngFromDegrees(30.430111, 115.42615)

ll9 := s2.LatLngFromDegrees(30.088491, 115.640384)

ll10 := s2.LatLngFromDegrees(29.907713, 115.656863)

ll11 := s2.LatLngFromDegrees(29.783833, 115.135012)

ll12 := s2.LatLngFromDegrees(29.712295, 114.728518)

ll13 := s2.LatLngFromDegrees(29.55473, 114.24512)

ll14 := s2.LatLngFromDegrees(29.530835, 113.717776)

ll15 := s2.LatLngFromDegrees(29.55473, 113.3772)

ll16 := s2.LatLngFromDegrees(29.678892, 112.998172)

ll17 := s2.LatLngFromDegrees(29.941039, 112.349978)

ll18 := s2.LatLngFromDegrees(30.040949, 112.025882)

ll19 := s2.LatLngFromDegrees(31.803269, 113.421145)

point1 := s2.PointFromLatLng(ll1)

point2 := s2.PointFromLatLng(ll2)

point3 := s2.PointFromLatLng(ll3)

point4 := s2.PointFromLatLng(ll4)

point5 := s2.PointFromLatLng(ll5)

point6 := s2.PointFromLatLng(ll6)

point7 := s2.PointFromLatLng(ll7)

point8 := s2.PointFromLatLng(ll8)

point9 := s2.PointFromLatLng(ll9)

point10 := s2.PointFromLatLng(ll10)

point11 := s2.PointFromLatLng(ll11)

point12 := s2.PointFromLatLng(ll12)

point13 := s2.PointFromLatLng(ll13)

point14 := s2.PointFromLatLng(ll14)

point15 := s2.PointFromLatLng(ll15)

point16 := s2.PointFromLatLng(ll16)

point17 := s2.PointFromLatLng(ll17)

point18 := s2.PointFromLatLng(ll18)

point19 := s2.PointFromLatLng(ll19)

points := []s2.Point{}

points = append(points, point19)

points = append(points, point18)

points = append(points, point17)

points = append(points, point16)

points = append(points, point15)

points = append(points, point14)

points = append(points, point13)

points = append(points, point12)

points = append(points, point11)

points = append(points, point10)

points = append(points, point9)

points = append(points, point8)

points = append(points, point7)

points = append(points, point6)

points = append(points, point5)

points = append(points, point4)

points = append(points, point3)

points = append(points, point2)

points = append(points, point1)

loop := s2.LoopFromPoints(points)

fmt.Println("---- loop search (gets too
much) -----")

// fmt.Printf(“Some loop status items: empty:%t full:%t \n”, loop.IsEmpty(),
loop.IsFull())

// ref: https://github.com/golang/geo/issues/14#issuecomment-257064823

defaultCoverer := &s2.RegionCoverer{MaxLevel: 20,
MaxCells: 1000, MinLevel: 1}

// rg := s2.Region(loop.CapBound())

// cvr := defaultCoverer.Covering(rg)

cvr := defaultCoverer.Covering(loop)

// fmt.Println(poly.CapBound())

for _, c3 :=
range cvr {

fmt.Printf("%d,\n", c3)

}

}

這里用到了 Loop 類，這個類的初始化的最小單元是 Point，Point 是由經緯度產生的。最重要的一點需要注意的是，多邊形是按照逆時針方向，左手邊區域確定的。如果一不小心點是按照順時針排列的話，那么多邊形確定的是外層更大的面，意味著球面除去畫的這個多邊形以外的都是你想要的多邊形。舉個具體的例子，假如我們想要畫的多邊形是下圖這個樣子的：

如果我們用順時針的方式依次存儲 Point 的話，并用順時針的這個數組去初始化 Loop，那么就會出現“奇怪”的現象。如下圖：

這張圖左上角的頂點和右下角的頂點在地球上是重合的。如果把這個地圖重新還原成球面，那么就是整個球面中間挖空了一個多邊形。把上圖放大，如下圖：

這樣就可以很清晰的看到了，中間被挖空了一個多邊形。造成這種現象的原因就是按照順時針的方向存儲了每個點，那么初始化一個 Loop 的時候就會選擇多邊形外圈的更大的多邊形。使用 Loop 一定要切記，順時針表示的是外圈多邊形，逆時針表示的是內圈多邊形。多邊形覆蓋的問題同之前舉的例子一樣：相同的 MaxLevel = 20，MinLevel = 1，MaxCells 不同，覆蓋的精度就不同，下圖是 MaxCells = 100 的情況：

下圖是 MaxCells = 1000 的情況：

6. S2的應用
   

S2 主要能用在以下 8 個地方：涉及到角度，間隔，緯度經度點，單位矢量等的表示，以及對這些類型的各種操作。單位球體上的幾何形狀，如球冠（“圓盤”），緯度 - 經度矩形，折線和多邊形。支持點，折線和多邊形的任意集合的強大的構造操作（例如聯合）和布爾謂詞（例如，包含）。對點，折線和多邊形的集合進行快速的內存索引。針對測量距離和查找附近物體的算法。用于捕捉和簡化幾何的穩健算法（該算法具有精度和拓撲保證）。用于測試幾何對象之間關系的有效且精確的數學謂詞的集合。支持空間索引，包括將區域近似為離散“S2單元”的集合。此功能可以輕松構建大型分布式空間索引。最后一點空間索引相信在工業生產中使用的非常廣泛。S2 目前應用比較多，用在和地圖相關業務上更多。Google Map 就直接大量使用了 S2 ，速度有多快讀者可以自己體驗體驗。Uber 在搜尋最近的出租車也是用的 S2 算法進行計算的。場景的例子就是本篇文章引語里面提到的場景。滴滴應該也有相關的應用，也許有更加優秀的解法。現在很火的共享單車也會用到這些空間索引算法。最后就是外賣行業和地圖關聯也很密切。美團和餓了么相信也在這方面有很多應用，具體哪里用到了，就請讀者自己想象吧。當然 S2 也有不適合的使用場景：平面幾何問題（有許多精細的現有平面幾何圖庫可供選擇）。轉換常見的 to/from GIS格式（要閱讀這種格式，請使用OGR等外部庫）。

五. 最后

本篇文章里面著重介紹了谷歌的 S2 算法的基礎實現。雖然 Geohash 也是空間點索引算法，但是性能方面比谷歌的 S2 略遜一籌。并且大公司的數據庫也基本上開始采用谷歌的 S2 算法進行索引。關于空間搜索其實還有一大類問題，如何搜索多維空間線，多維空間面，多維空間多邊形呢？他們都是由無數個空間點組成的。實際的例子，比如街道，高樓，鐵路，河流。要搜索這些東西，數據庫表如何設計？如何做到高效的搜索呢？還能用 B+ 樹來做么？答案當然是也可以實現高效率的搜索，那就需要用到 R 樹，或者 R 樹 和 B+樹。這部分就不在本文的范疇內了，可參考《多維空間多邊形索引算法》。

GitHub Repo：Halfrost-Field

Follow: halfrost · GitHub

Source:
https://halfrost.com/go_spatial_search/

總結

以上是生活随笔為你收集整理的空间点像素索引（三）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。