空间点像素索引（三）

点与坐标轴点相互转换

在 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² 。

```
与Geohash对比
```

Geohash有12级，从5000km 到3.7cm。中间每一级的变化比较大。有时候可能选择上一级会大很多，选择下一级又会小一些。比如选择字符串长度为4，它对应的 cell 宽度是39.1km，需求可能是50km，那么选择字符串长度为5，对应的 cell 宽度就变成了156km，瞬间又大了3倍了。这种情况选择多长的 Geohash 字符串就比较难选。选择不好，每次判断可能就还需要取出周围的8个格子再次进行判断。Geohash 需要 12 bytes 存储。S2有30级，从 0.7cm² 到 85,000,000km² 。中间每一级的变化都比较平缓，接近于4次方的曲线。所以选择精度不会出现Geohash 选择困难的问题。S2 的存储只需要一个uint64 即可存下。S2库里面不仅仅有地理编码，还有其他很多几何计算相关的库。地理编码只是其中的一小部分。本文没有介绍到的 S2 的实现还有很多很多，各种向量计算，面积计算，多边形覆盖，距离问题，球面球体上的问题，它都有实现。S2还能解决多边形覆盖的问题。比如给定一个城市，求一个多边形刚刚好覆盖住这个城市。

如上图，生成的多边形刚刚好覆盖住下面蓝色的区域。这里生成的多边形可以有大有小。不管怎么样，最终的结果也是刚刚覆盖住目标物。

用相同的 Cell 也可以达到相同的目的，上图就是用相同 Level 的 Cell 覆盖了整个圣保罗城市。这些都是 Geohash 做不到的。多边形覆盖利用的是近似的算法，虽然不是严格意义上的最优解，但是实践中效果特别好。额外值得说明的一点是，Google 文档上强调了，这种多边形覆盖的算法虽然对搜索和预处理操作非常有用，但是“不可依赖”的。理由也是因为是近似算法，并不是唯一最优算法，所以得到的解会依据库的不同版本而产生变化。

```
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 的情况：

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

空间点像素索引（三）相关推荐

空间点像素索引（二）
空间点像素索引(二) 三. Hilbert Curve 希尔伯特曲线希尔伯特曲线的定义希尔伯特曲线一种能填充满一个平面正方形的分形曲线(空间填充曲线),由大卫·希尔伯特在1891年提出.由于它能填 ...
空间点像素索引（一）
空间点像素索引(一) 给你一个需求,查找定位点附近的一定范围内所有餐馆,你会怎么实现呢? 一一计算定位点与所有餐馆的距离,然后取出最小的距离?如果同时有很多遍布全国的请求都在查找附近的餐馆,按照上述的 ...
【opencv学习笔记】第五篇：访问图像中像素的三种方式、ROI区域图像叠加和图像混合
1. 访问图像中像素的三种方式任何图像处理算法,都是从操作每个像素开始的.在OpenCV中,提供了三种访问每个像素的方法. 方法1:指针访问:C操作符[] 方法2:迭代器iterator 方法3:动 ...
Opencv3编程入门学习笔记（三）之访问图像像素的三种方法
访问图像像素的三种方法:指针访问,迭代器访问,动态地址访问.访问最快的为指针访问,以下算法在几毫秒,但指针访问容易造成内存泄漏:其次为迭代器访问:最后为动态地址访问. 以下程序是根据<OpenC ...
OpenCV精进之路（零）：访问图像中像素的三种方法
访问像素的三种方法指针访问:最快迭代器iterator:较慢,非常安全,指针访问可能出现越界问题动态地址计算:更慢,通过at()实现.适用于访问具体某个第i行,j列的像素,而不适用遍历像素这里 ...
mysql索引三个字段查询两个字段_mysql中关于关联索引的问题——对a,b,c三个字段建立联合索引，那么查询时使用其中的2个作为查询条件，是否还会走索引？...
情况描述:在MySQL的user表中,对a,b,c三个字段建立联合索引,那么查询时使用其中的2个作为查询条件,是否还会走索引? 根据查询字段的位置不同来决定,如查询a, a,b a,b, ...
mysql索引三个字段查询两个字段_mysql中关于关联索引的问题——对a,b,c三个字段建立联合索引，那么查询时使用其中的2个作为查询条件，是否还会走索引？......
情况描述:在MySQL的user表中,对a,b,c三个字段建立联合索引,那么查询时使用其中的2个作为查询条件,是否还会走索引? 根据查询字段的位置不同来决定,如查询a, a,b a,b, ...
OpenCV（图像处理）—访问像素的三种方法
方法一:用指针访问像素 1 #include <opencv2/opencv.hpp> 2 #include <opencv2/core/core.hpp> 3 #includ ...
react如何获取option的下标和值_select选中获取索引三种写法
$('#someId').prop('selectedIndex'); $('option:selected', '#someId').index(); $('#someId option').ind ...

空间点像素索引（三）

空间点像素索引（三）相关推荐

最新文章

热门文章