空間點(diǎn)像素索引（二）

三. Hilbert Curve
希爾伯特曲線

1. 希爾伯特曲線的定義

希爾伯特曲線一種能填充滿一個(gè)平面正方形的分形曲線（空間填充曲線），由大衛(wèi)·希爾伯特在1891年提出。由于它能填滿平面，它的豪斯多夫維是2。取它填充的正方形的邊長(zhǎng)為1，第n步的希爾伯特曲線的長(zhǎng)度是2^n - 2^(-n)。

2. 希爾伯特曲線的構(gòu)造方法

一階的希爾伯特曲線，生成方法就是把正方形四等分，從其中一個(gè)子正方形的中心開始，依次穿線，穿過其余3個(gè)正方形的中心。二階的希爾伯特曲線，生成方法就是把之前每個(gè)子正方形繼續(xù)四等分，每4個(gè)小的正方形先生成一階希爾伯特曲線。然后把4個(gè)一階的希爾伯特曲線首尾相連。三階的希爾伯特曲線，生成方法就是與二階類似，先生成二階希爾伯特曲線。然后把4個(gè)二階的希爾伯特曲線首尾相連。n階的希爾伯特曲線的生成方法也是遞歸的，先生成n-1階的希爾伯特曲線，然后把4個(gè)n-1階的希爾伯特曲線首尾相連。

3. 為何要選希爾伯特曲線？

看到這里可能就有讀者有疑問了，這么多空間填充曲線，為何要選希爾伯特曲線？因?yàn)橄柌厍€有非常好的特性。

(1) 降維

首先，作為空間填充曲線，希爾伯特曲線可以對(duì)多維空間有效的降維。上圖就是希爾伯特曲線在填滿一個(gè)平面以后，把平面上的點(diǎn)都展開成一維的線了。可能有人會(huì)有疑問，上圖里面的希爾伯特曲線只穿了16個(gè)點(diǎn)，怎么能代表一個(gè)平面呢？

當(dāng)然，當(dāng)n趨近于無窮大的時(shí)候，n階希爾伯特曲線就可以近似填滿整個(gè)平面了。

(2) 穩(wěn)定

當(dāng)n階希爾伯特曲線，n趨于無窮大的時(shí)候，曲線上的點(diǎn)的位置基本上趨于穩(wěn)定。舉個(gè)例子：上圖左邊是希爾伯特曲線，右邊是蛇形的曲線。當(dāng)n趨于無窮大的時(shí)候，兩者理論上都可以填滿平面。但是為何希爾伯特曲線更加優(yōu)秀呢？在蛇形曲線上給定一個(gè)點(diǎn)，當(dāng)n趨于無窮大的過程中，這個(gè)點(diǎn)在蛇形曲線上的位置是時(shí)刻變化的。這就造成了點(diǎn)的相對(duì)位置始終不定。再看看希爾伯特曲線，同樣是一個(gè)點(diǎn)，在n趨于無窮大的情況下：從上圖可以看到，點(diǎn)的位置幾乎沒有怎么變化。所以希爾伯特曲線更加優(yōu)秀。

(3) 連續(xù)

希爾伯特曲線是連續(xù)的，所以能保證一定可以填滿空間。連續(xù)性是需要數(shù)學(xué)證明的。具體證明方法這里就不細(xì)說了，感興趣的可以點(diǎn)文章末尾一篇關(guān)于希爾伯特曲線的論文，那里有連續(xù)性的證明。接下來要介紹的谷歌的 S2 算法就是基于希爾伯特曲線的。現(xiàn)在讀者應(yīng)該明白選擇希爾伯特曲線的原因了吧。

四. 算法

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 算法發(fā)布4年沒有得到它應(yīng)有的贊賞。不過現(xiàn)在 S2 已經(jīng)被各大公司使用了。在介紹這個(gè)重量級(jí)算法之前，先解釋一些這個(gè)算法的名字由來。S2其實(shí)是來自幾何數(shù)學(xué)中的一個(gè)數(shù)學(xué)符號(hào) S2，它表示的是單位球。S2 這個(gè)庫其實(shí)是被設(shè)計(jì)用來解決球面上各種幾何問題的。值得提的一點(diǎn)是，除去 golang 官方 repo 里面的 geo/s2 完成度目前只有40%，其他語言，Java，C++，Python 的 S2 實(shí)現(xiàn)都完成100%了。本篇文章講解以 Go 的這個(gè)版本為主。接下來就看看怎么用 S2 來解決多維空間點(diǎn)索引的問題的。

1. 球面坐標(biāo)轉(zhuǎn)換

按照之前我們處理多維空間的思路，先考慮如何降維，再考慮如何分形。眾所周知，地球是近似一個(gè)球體。球體是一個(gè)三維的，如何把三維降成一維呢？球面上的一個(gè)點(diǎn)，在直角坐標(biāo)系中，可以這樣表示：

x = r * sin θ * cos φ

y = r * sin θ * sin φ

z = r * cos θ

通常地球上的點(diǎn)我們會(huì)用經(jīng)緯度來表示。

再進(jìn)一步，我們可以和球面上的經(jīng)緯度聯(lián)系起來。不過這里需要注意的是，緯度的角度 α 和直角坐標(biāo)系下的球面坐標(biāo) θ 加起來等于 90°。所以三角函數(shù)要注意轉(zhuǎn)換。于是地球上任意的一個(gè)經(jīng)緯度的點(diǎn)，就可以轉(zhuǎn)換成 f(x,y,z)。在 S2 中，地球半徑被當(dāng)成單位 1 了。所以半徑不用考慮。x，y，z的值域都被限定在了[-1,1] x [-1,1] x [-1,1]這個(gè)區(qū)間之內(nèi)了。

2. 球面變平面

從球心向外切正方體6個(gè)面分別投影。S2 是把球面上所有的點(diǎn)都投影到外切正方體的6個(gè)面上。

這里簡(jiǎn)單的畫了一個(gè)投影圖，上圖左邊的是投影到正方體一個(gè)面的示意圖，實(shí)際上影響到的球面是右邊那張圖。

從側(cè)面看，其中一個(gè)球面投影到正方體其中一個(gè)面上，邊緣與圓心的連線相互之間的夾角為90°，但是和x，y，z軸的角度是45°。我們可以在球的6個(gè)方向上，把45°的輔助圓畫出來，見下圖左邊。

上圖左邊的圖畫了6個(gè)輔助線，藍(lán)線是前后一對(duì)，紅線是左右一對(duì)，綠線是上下一對(duì)。分別都是45°的地方和圓心連線與球面相交的點(diǎn)的軌跡。這樣我們就可以把投影到外切正方體6個(gè)面上的球面畫出來，見上圖右邊。投影到正方體以后，我們就可以把這個(gè)正方體展開了。一個(gè)正方體的展開方式有很多種。不管怎么展開，最小單元都是一個(gè)正方形。以上就是 S2 的投影方案。接下來講講其他的投影方案。首先有下面一種方式，三角形和正方形組合。這種方式展開圖如下圖。這種方式其實(shí)很復(fù)雜，構(gòu)成子圖形由兩種圖形構(gòu)成。坐標(biāo)轉(zhuǎn)換稍微復(fù)雜一點(diǎn)。再還有一種方式是全部用三角形組成，這種方式三角形個(gè)數(shù)越多，就能越切近于球體。

上圖最左邊的圖，由20個(gè)三角形構(gòu)成，可以看的出來，菱角非常多，與球體相差比較大，當(dāng)三角形個(gè)數(shù)越來越多，就越來越貼近球體。

20個(gè)三角形展開以后就可能變成這樣。最后一種方式可能是目前最好的方式，不過也可能是最復(fù)雜的方式。按照六邊形來投影。

六邊形的菱角比較少，六個(gè)邊也能相互銜接其他的六邊形。看上圖最后邊的圖可以看出來，六邊形足夠多以后，非常近似球體。

六邊形展開以后就是上面這樣。當(dāng)然這里只有12個(gè)六邊形。六邊形個(gè)數(shù)越多越好，粒度越細(xì)，就越貼近球體。Uber 在一個(gè)公開分享上提到了他們用的是六邊形的網(wǎng)格，把城市劃分為很多六邊形。這塊應(yīng)該是他們自己開發(fā)的。也許滴滴也是劃分六邊形，也許滴滴有更好的劃分方案也說不定。在 Google S2 中，它是把地球展開成如下的樣子：

如果上面展開的6個(gè)面，假設(shè)都用5階的希爾伯特曲線表示出來的話，6個(gè)面會(huì)是如下的樣子：

回到 S2 上面來，S2是用的正方形。這樣第一步的球面坐標(biāo)進(jìn)一步的被轉(zhuǎn)換成 f(x,y,z) -> g(face,u,v)，face是正方形的六個(gè)面，u，v對(duì)應(yīng)的是六個(gè)面中的一個(gè)面上的x，y坐標(biāo)。

3. 球面矩形投影修正

上一步我們把球面上的球面矩形投影到正方形的某個(gè)面上，形成的形狀類似于矩形，但是由于球面上角度的不同，最終會(huì)導(dǎo)致即使是投影到同一個(gè)面上，每個(gè)矩形的面積也不大相同。

上圖就表示出了球面上個(gè)一個(gè)球面矩形投影到正方形一個(gè)面上的情況。

經(jīng)過實(shí)際計(jì)算發(fā)現(xiàn)，最大的面積和最小的面積相差5.2倍。見上圖左邊。相同的弧度區(qū)間，在不同的緯度上投影到正方形上的面積不同。現(xiàn)在就需要修正各個(gè)投影出來形狀的面積。如何選取合適的映射修正函數(shù)就成了關(guān)鍵。目標(biāo)是能達(dá)到上圖右邊的樣子，讓各個(gè)矩形的面積盡量相同。這塊轉(zhuǎn)換的代碼在 C++ 的版本里面才有詳細(xì)的解釋，在 Go 的版本里面只一筆帶過了。害筆者懵逼了好久。

線性變換是最快的變換，但是變換比最小。tan() 變換可以使每個(gè)投影以后的矩形的面積更加一致，最大和最小的矩形比例僅僅只差0.414。可以說非常接近了。但是 tan() 函數(shù)的調(diào)用時(shí)間非常長(zhǎng)。如果把所有點(diǎn)都按照這種方式計(jì)算的話，性能將會(huì)降低3倍。最后谷歌選擇的是二次變換，這是一個(gè)近似切線的投影曲線。它的計(jì)算速度遠(yuǎn)遠(yuǎn)快于 tan() ，大概是 tan() 計(jì)算的3倍速度。生成的投影以后的矩形大小也類似。不過最大的矩形和最小的矩形相比依舊有2.082的比率。上表中，ToPoint 和 FromPoint 分別是把單位向量轉(zhuǎn)換到 Cell ID 所需要的毫秒數(shù)、把 Cell ID 轉(zhuǎn)換回單位向量所需的毫秒數(shù)（Cell ID 就是投影到正方體六個(gè)面，某個(gè)面上矩形的 ID，矩形稱為 Cell，它對(duì)應(yīng)的 ID 稱為 Cell ID）。ToPointRaw 是某種目的下，把 Cell ID 轉(zhuǎn)換為非單位向量所需的毫秒數(shù)。在 S2 中默認(rèn)的轉(zhuǎn)換是二次轉(zhuǎn)換。

#define S2_PROJECTION S2_QUADRATIC_PROJECTION

詳細(xì)看看這三種轉(zhuǎn)換到底是怎么轉(zhuǎn)換的。

#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

上面有一處對(duì) tan(M_PI_4) 的處理，是因?yàn)榫鹊脑?#xff0c;導(dǎo)致略小于1.0 。所以投影之后的修正函數(shù)三種變換應(yīng)該如下：

// 線性轉(zhuǎn)換

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。實(shí)際使用過程中，u，v都分別當(dāng)做入?yún)?#xff0c;都會(huì)進(jìn)行變換。這塊修正函數(shù)在 Go 的版本里面就直接只實(shí)現(xiàn)了二次變換，其他兩種變換方式找遍整個(gè)庫，根本沒有提及。

// 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)

}

經(jīng)過修正變換以后，u，v都變換成了s，t。值域也發(fā)生了變化。u，v的值域是[-1,1]，變換以后，是s，t的值域是[0,1]。至此，小結(jié)一下，球面上的點(diǎn)S(lat,lng) -> f(x,y,z) -> g(face,u,v)
-> h(face,s,t)。目前總共轉(zhuǎn)換了4步，球面經(jīng)緯度坐標(biāo)轉(zhuǎn)換成球面xyz坐標(biāo)，再轉(zhuǎn)換成外切正方體投影面上的坐標(biāo)，最后變換成修正后的坐標(biāo)。到目前為止，S2 可以優(yōu)化的點(diǎn)有兩處，一是投影的形狀能否換成六邊形？二是修正的變換函數(shù)能否找到一個(gè)效果和 tan() 類似的函數(shù)，但是計(jì)算速度遠(yuǎn)遠(yuǎn)高于 tan()，以致于不會(huì)影響計(jì)算性能？

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。