空间点像素索引（二）
三. Hilbert Curve 希尔伯特曲线
1. 希尔伯特曲线的定义
希尔伯特曲线一种能填充满一个平面正方形的分形曲线（空间填充曲线），由大卫·希尔伯特在1891年提出。由于它能填满平面，它的豪斯多夫维是2。取它填充的正方形的边长为1，第n步的希尔伯特曲线的长度是2^n - 2^(-n)。
2. 希尔伯特曲线的构造方法
一阶的希尔伯特曲线，生成方法就是把正方形四等分，从其中一个子正方形的中心开始，依次穿线，穿过其余3个正方形的中心。二阶的希尔伯特曲线，生成方法就是把之前每个子正方形继续四等分，每4个小的正方形先生成一阶希尔伯特曲线。然后把4个一阶的希尔伯特曲线首尾相连。三阶的希尔伯特曲线，生成方法就是与二阶类似，先生成二阶希尔伯特曲线。然后把4个二阶的希尔伯特曲线首尾相连。n阶的希尔伯特曲线的生成方法也是递归的，先生成n-1阶的希尔伯特曲线，然后把4个n-1阶的希尔伯特曲线首尾相连。
3. 为何要选希尔伯特曲线？
看到这里可能就有读者有疑问了，这么多空间填充曲线，为何要选希尔伯特曲线？因为希尔伯特曲线有非常好的特性。
(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其实是来自几何数学中的一个数学符号 S²，它表示的是单位球。S2 这个库其实是被设计用来解决球面上各种几何问题的。值得提的一点是，除去 golang 官方 repo 里面的 geo/s2 完成度目前只有40%，其他语言，Java，C++，Python 的 S2 实现都完成100%了。本篇文章讲解以 Go 的这个版本为主。接下来就看看怎么用 S2 来解决多维空间点索引的问题的。
1. 球面坐标转换
按照之前我们处理多维空间的思路，先考虑如何降维，再考虑如何分形。众所周知，地球是近似一个球体。球体是一个三维的，如何把三维降成一维呢？球面上的一个点，在直角坐标系中，可以这样表示：
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]这个区间之内了。
2. 球面变平面
从球心向外切正方体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坐标。
3. 球面矩形投影修正
上一步我们把球面上的球面矩形投影到正方形的某个面上，形成的形状类似于矩形，但是由于球面上角度的不同，最终会导致即使是投影到同一个面上，每个矩形的面积也不大相同。
上图就表示出了球面上个一个球面矩形投影到正方形一个面上的情况。
经过实际计算发现，最大的面积和最小的面积相差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.) * (4*s*s - 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 - 3*u);
}
#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 - 3*u)
注意上面虽然变换公式只写了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.) * (4*s*s - 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.5*math.Sqrt(1-3*u)
}
经过修正变换以后，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()，以致于不会影响计算性能？