空间点像素索引（二）
空間點像素索引（二）
三. Hilbert Curve
希爾伯特曲線
- 希爾伯特曲線的定義
希爾伯特曲線一種能填充滿一個平面正方形的分形曲線（空間填充曲線），由大衛·希爾伯特在1891年提出。由于它能填滿平面，它的豪斯多夫維是2。取它填充的正方形的邊長為1，第n步的希爾伯特曲線的長度是2^n - 2^(-n)。
- 希爾伯特曲線的構造方法
一階的希爾伯特曲線，生成方法就是把正方形四等分，從其中一個子正方形的中心開始，依次穿線，穿過其余3個正方形的中心。二階的希爾伯特曲線，生成方法就是把之前每個子正方形繼續四等分，每4個小的正方形先生成一階希爾伯特曲線。然后把4個一階的希爾伯特曲線首尾相連。三階的希爾伯特曲線，生成方法就是與二階類似，先生成二階希爾伯特曲線。然后把4個二階的希爾伯特曲線首尾相連。n階的希爾伯特曲線的生成方法也是遞歸的，先生成n-1階的希爾伯特曲線，然后把4個n-1階的希爾伯特曲線首尾相連。
- 為何要選希爾伯特曲線？
看到這里可能就有讀者有疑問了，這么多空間填充曲線，為何要選希爾伯特曲線？因為希爾伯特曲線有非常好的特性。
(1) 降維
首先，作為空間填充曲線，希爾伯特曲線可以對多維空間有效的降維。上圖就是希爾伯特曲線在填滿一個平面以后，把平面上的點都展開成一維的線了。可能有人會有疑問，上圖里面的希爾伯特曲線只穿了16個點，怎么能代表一個平面呢？
當然，當n趨近于無窮大的時候，n階希爾伯特曲線就可以近似填滿整個平面了。
(2) 穩定
當n階希爾伯特曲線，n趨于無窮大的時候，曲線上的點的位置基本上趨于穩定。舉個例子：上圖左邊是希爾伯特曲線，右邊是蛇形的曲線。當n趨于無窮大的時候，兩者理論上都可以填滿平面。但是為何希爾伯特曲線更加優秀呢？在蛇形曲線上給定一個點，當n趨于無窮大的過程中，這個點在蛇形曲線上的位置是時刻變化的。這就造成了點的相對位置始終不定。再看看希爾伯特曲線，同樣是一個點，在n趨于無窮大的情況下：從上圖可以看到，點的位置幾乎沒有怎么變化。所以希爾伯特曲線更加優秀。
(3) 連續
希爾伯特曲線是連續的，所以能保證一定可以填滿空間。連續性是需要數學證明的。具體證明方法這里就不細說了，感興趣的可以點文章末尾一篇關于希爾伯特曲線的論文，那里有連續性的證明。接下來要介紹的谷歌的 S2 算法就是基于希爾伯特曲線的。現在讀者應該明白選擇希爾伯特曲線的原因了吧。
四. 算法
Google’s S2 library is a real treasure, not only due to its
capabilities for spatial indexing but also because it is a library that was
released more than 4 years ago and it didn’t get the attention it deserved.
上面這段話來自2015年一位谷歌工程師的博文。他由衷的感嘆 S2 算法發布4年沒有得到它應有的贊賞。不過現在 S2 已經被各大公司使用了。在介紹這個重量級算法之前，先解釋一些這個算法的名字由來。S2其實是來自幾何數學中的一個數學符號 S2，它表示的是單位球。S2 這個庫其實是被設計用來解決球面上各種幾何問題的。值得提的一點是，除去 golang 官方 repo 里面的 geo/s2 完成度目前只有40%，其他語言，Java，C++，Python 的 S2 實現都完成100%了。本篇文章講解以 Go 的這個版本為主。接下來就看看怎么用 S2 來解決多維空間點索引的問題的。
- 球面坐標轉換
按照之前我們處理多維空間的思路，先考慮如何降維，再考慮如何分形。眾所周知，地球是近似一個球體。球體是一個三維的，如何把三維降成一維呢？球面上的一個點，在直角坐標系中，可以這樣表示：
x = r * sin θ * cos φ
y = r * sin θ * sin φ
z = r * cos θ
通常地球上的點我們會用經緯度來表示。
再進一步，我們可以和球面上的經緯度聯系起來。不過這里需要注意的是，緯度的角度 α 和直角坐標系下的球面坐標 θ 加起來等于 90°。所以三角函數要注意轉換。于是地球上任意的一個經緯度的點，就可以轉換成 f(x,y,z)。在 S2 中，地球半徑被當成單位 1 了。所以半徑不用考慮。x，y，z的值域都被限定在了[-1,1] x [-1,1] x [-1,1]這個區間之內了。
- 球面變平面
從球心向外切正方體6個面分別投影。S2 是把球面上所有的點都投影到外切正方體的6個面上。
這里簡單的畫了一個投影圖，上圖左邊的是投影到正方體一個面的示意圖，實際上影響到的球面是右邊那張圖。
從側面看，其中一個球面投影到正方體其中一個面上，邊緣與圓心的連線相互之間的夾角為90°，但是和x，y，z軸的角度是45°。我們可以在球的6個方向上，把45°的輔助圓畫出來，見下圖左邊。
上圖左邊的圖畫了6個輔助線，藍線是前后一對，紅線是左右一對，綠線是上下一對。分別都是45°的地方和圓心連線與球面相交的點的軌跡。這樣我們就可以把投影到外切正方體6個面上的球面畫出來，見上圖右邊。投影到正方體以后，我們就可以把這個正方體展開了。一個正方體的展開方式有很多種。不管怎么展開，最小單元都是一個正方形。以上就是 S2 的投影方案。接下來講講其他的投影方案。首先有下面一種方式，三角形和正方形組合。這種方式展開圖如下圖。這種方式其實很復雜，構成子圖形由兩種圖形構成。坐標轉換稍微復雜一點。再還有一種方式是全部用三角形組成，這種方式三角形個數越多，就能越切近于球體。
上圖最左邊的圖，由20個三角形構成，可以看的出來，菱角非常多，與球體相差比較大，當三角形個數越來越多，就越來越貼近球體。
20個三角形展開以后就可能變成這樣。最后一種方式可能是目前最好的方式，不過也可能是最復雜的方式。按照六邊形來投影。
六邊形的菱角比較少，六個邊也能相互銜接其他的六邊形。看上圖最后邊的圖可以看出來，六邊形足夠多以后，非常近似球體。
六邊形展開以后就是上面這樣。當然這里只有12個六邊形。六邊形個數越多越好，粒度越細，就越貼近球體。Uber 在一個公開分享上提到了他們用的是六邊形的網格，把城市劃分為很多六邊形。這塊應該是他們自己開發的。也許滴滴也是劃分六邊形，也許滴滴有更好的劃分方案也說不定。在 Google S2 中，它是把地球展開成如下的樣子：
如果上面展開的6個面，假設都用5階的希爾伯特曲線表示出來的話，6個面會是如下的樣子：
回到 S2 上面來，S2是用的正方形。這樣第一步的球面坐標進一步的被轉換成 f(x,y,z) -> g(face,u,v)，face是正方形的六個面，u，v對應的是六個面中的一個面上的x，y坐標。
- 球面矩形投影修正
上一步我們把球面上的球面矩形投影到正方形的某個面上，形成的形狀類似于矩形，但是由于球面上角度的不同，最終會導致即使是投影到同一個面上，每個矩形的面積也不大相同。
上圖就表示出了球面上個一個球面矩形投影到正方形一個面上的情況。
經過實際計算發現，最大的面積和最小的面積相差5.2倍。見上圖左邊。相同的弧度區間，在不同的緯度上投影到正方形上的面積不同。現在就需要修正各個投影出來形狀的面積。如何選取合適的映射修正函數就成了關鍵。目標是能達到上圖右邊的樣子，讓各個矩形的面積盡量相同。這塊轉換的代碼在 C++ 的版本里面才有詳細的解釋，在 Go 的版本里面只一筆帶過了。害筆者懵逼了好久。
線性變換是最快的變換，但是變換比最小。tan() 變換可以使每個投影以后的矩形的面積更加一致，最大和最小的矩形比例僅僅只差0.414。可以說非常接近了。但是 tan() 函數的調用時間非常長。如果把所有點都按照這種方式計算的話，性能將會降低3倍。最后谷歌選擇的是二次變換，這是一個近似切線的投影曲線。它的計算速度遠遠快于 tan() ，大概是 tan() 計算的3倍速度。生成的投影以后的矩形大小也類似。不過最大的矩形和最小的矩形相比依舊有2.082的比率。上表中，ToPoint 和 FromPoint 分別是把單位向量轉換到 Cell ID 所需要的毫秒數、把 Cell ID 轉換回單位向量所需的毫秒數（Cell ID 就是投影到正方體六個面，某個面上矩形的 ID，矩形稱為 Cell，它對應的 ID 稱為 Cell ID）。ToPointRaw 是某種目的下，把 Cell ID 轉換為非單位向量所需的毫秒數。在 S2 中默認的轉換是二次轉換。
#define S2_PROJECTION S2_QUADRATIC_PROJECTION
詳細看看這三種轉換到底是怎么轉換的。
#if S2_PROJECTION == S2_LINEAR_PROJECTION
inline double S2::STtoUV(double s) {
return
2 * s - 1;
}
inline double S2::UVtoST(double u) {
return
0.5 * (u + 1);
}
#elif S2_PROJECTION == S2_TAN_PROJECTION
inline double S2::STtoUV(double s) {
//
Unfortunately, tan(M_PI_4) is slightly less than 1.0. This isn’t due to
// a
flaw in the implementation of tan(), it’s because the derivative of
//
tan(x) at x=pi/4 is 2, and it happens that the two adjacent floating
//
point numbers on either side of the infinite-precision value of pi/4 have
//
tangents that are slightly below and slightly above 1.0 when rounded to
// the
nearest double-precision result.
s =tan(M_PI_2 * s - M_PI_4);
return
s + (1.0 / (GG_LONGLONG(1) << 53)) * s;
}
inline double S2::UVtoST(double u) {
volatile double a = atan(u);
return
(2 * M_1_PI) * (a + M_PI_4);
}
#elif S2_PROJECTION == S2_QUADRATIC_PROJECTION
inline double S2::STtoUV(double s) {
if (s
= 0.5) return (1/3.) * (4ss - 1);
else return (1/3.) * (1 -4*(1-s)*(1-s));
}
inline double S2::UVtoST(double u) {
if (u>= 0) return 0.5 * sqrt(1 + 3*u);
else return 1 - 0.5 sqrt(1 - 3u);
}
#else
#error Unknown value for S2_PROJECTION
#endif
上面有一處對 tan(M_PI_4) 的處理，是因為精度的原因，導致略小于1.0 。所以投影之后的修正函數三種變換應該如下：
// 線性轉換
u = 0.5 * ( u + 1)
// tan() 變換
u = 2 / pi * (atan(u) + pi / 4) = 2 * atan(u)
/ pi + 0.5
// 二次變換
u >= 0，u = 0.5 * sqrt(1+ 3*u)
u < 0， u = 1 - 0.5 sqrt(1 - 3u)
注意上面雖然變換公式只寫了u，不代表只變換u。實際使用過程中，u，v都分別當做入參，都會進行變換。這塊修正函數在 Go 的版本里面就直接只實現了二次變換，其他兩種變換方式找遍整個庫，根本沒有提及。
// stToUV converts an s or t value to the corresponding u
or v value.
// This is a non-linear transformation from
[-1,1] to [-1,1] that
// attempts to make the cell sizes more
uniform.
// This uses what the C++ version calls ‘the
quadratic transform’.
func stToUV(s float64) float64 {
if s >= 0.5 {
return (1 / 3.) * (4ss - 1)
}
return (1 / 3.) * (1 - 4*(1-s)*(1-s))
}
// uvToST is the inverse of the stToUV
transformation. Note that it
// is not always true that uvToST(stToUV(x))
== x due to numerical
// errors.
func uvToST(u float64) float64 {
if u >= 0 {
return 0.5 * math.Sqrt(1+3*u)
}
return 1 - 0.5math.Sqrt(1-3u)
}
經過修正變換以后，u，v都變換成了s，t。值域也發生了變化。u，v的值域是[-1,1]，變換以后，是s，t的值域是[0,1]。至此，小結一下，球面上的點S(lat,lng) -> f(x,y,z) -> g(face,u,v)
-> h(face,s,t)。目前總共轉換了4步，球面經緯度坐標轉換成球面xyz坐標，再轉換成外切正方體投影面上的坐標，最后變換成修正后的坐標。到目前為止，S2 可以優化的點有兩處，一是投影的形狀能否換成六邊形？二是修正的變換函數能否找到一個效果和 tan() 類似的函數，但是計算速度遠遠高于 tan()，以致于不會影響計算性能？
總結
以上是生活随笔為你收集整理的空间点像素索引（二）的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 空间点像素索引（一）
- 下一篇： 空间点像素索引（三）