Google S2 是如何解决空间覆盖最优解问题的?

这篇不出意外就是 Google S2 整个系列的最终篇了。这篇里面会把 regionCoverer 算法都讲解清楚。至于 Google S2 库里面还有很多其他的小算法，代码同样也很值得阅读和学习，这里也就不一一展开了，有兴趣的读者可以把整个库都读一遍。

在写这篇文章的同时，发现了库的作者的一些新的 commit ，比如 2017年12月4号的 commit f9610db2b871b54b17d36d4da6a4d6a2aab6018d，这次提交的改动更改了 README，别看只改了文档，其实里面内容很多。

可以看到他们把原来存在 Google 官方私有代码仓库里面的代码放到了 github。之前都只能在代码归档里面查看 C++ 代码，现在直接可以在 github 上查看了。方便了很多。

在这次 commit 里面还提到了一个新的网站，这个网站我发现也是最近才发布出来的。因为笔者我连续关注 S2 每个 commit 了快半年了。网上每个关于 S2 的资源都有关注，这个网站非常新。关于这个网站文章最后会提到，这里就不详细说了。

从这个提交开始，笔者认为 Google S2 可能被重视起来了，也有可能打算大力推广了。

好了，正式进入正题。

在 Google S2 中能进行 RegionCovering 的有以下几种类型。基本上都是必须满足 Region interface 的。

Cap 代表由中心和半径限定的盘形区域。从技术上讲，这种形状被称为“球形帽”（而不是圆盘），因为它不是平面的。帽子代表被飞机切断的球体的一部分。帽的边界是由球面与平面的交点所定义的圆。帽子是一个封闭的组合，即它包含了它的边界。 大多数情况下，无论在平面几何中使用光盘，都可以使用球冠。帽的半径是沿着球体表面测量的（而不是通过内部的直线距离）。因此，一个半径为 π/ 2 的帽是一个半球，半径为 π 的帽覆盖整个球。 中心是单位球面上的一个点。（因此需要它是单位长度）帽子也可以由其中心点和高度来定义。高度是从中心点到截断平面的距离。还有支持“空”和“全”的上限，分别不包含任何点数和所有点数。 下面是帽高（h），帽半径（r），帽中心的最大弦长（d）和帽底部半径（a）之间的一些有用关系。

Loop 代表一个简单的球面多边形。它由一系列顶点组成，其中第一个顶点隐含地被认为是连接到最后一个顶点的。所有的 loop 被定义为具有 CCW 方向，即 loop 的内部在边的左侧。这意味着包围一个小区域的顺时针 loop 被解释为包围非常大的区域的 CCW 的 loop。 loop 不允许有任何重复的顶点（不管是否相邻）。不允许相邻的边相交，而且不允许长度为180度的边（即，相邻的顶点不能是相反的）。loop 必须至少有3个顶点（下面讨论的“空”和“全” loop 除外）。 有两个特殊的 loop：EmptyLoop 不包含点，FullLoop 包含所有点。这些 loop 没有任何边，但为了保持每一个 loop 都可以表示为顶点链的不变量，它们被定义为每个只有一个顶点。

多边形表示一个零或多个 loop 的序列;同样，一个 loop 的左手边方向定义为它的内部。 当多边形初始化时，给定的 loop 自动转换为“孔”的组成的规范形式。loop 被重新排序以对应于嵌套层次的预定义遍历方式。 多边形可以表示具有多边形边界的球体的任何区域，包括整个球体（称为“完整”多边形）。完整的多边形由一个完整的 loop 组成，而空的多边形完全没有 loop。 使用 FullPolygon() 来构造一个完整的多边形。 Polygon 的零值被视为空的多边形。

想要 多个 loop 构成一个 Polygon 多边形，必须满足以下4个条件：

Rect 代表一个封闭的经纬度矩形。它也是 Region 类型。它能够表示空的和完整的矩形以及单个点。它有一个 AddPoint 方法，可以方便地为一组点构造边界矩形，包括跨越180度子午线的点集。

区域表示单位球体上的二维区域。 这个接口的目的是让复杂的区域近似为更简单的区域。该接口仅限于计算近似值的方法。

S2 区域表示单位球体上的二维区域。它是一个具有各种具体子类型的抽象接口，如盘形，矩形，多段线，多边形，几何集合，缓冲形状等。 这个接口的主要目的是使复杂区域近似为更简单的区域。因此，接口只能用于计算近似值的方法，而不是具有各种各样的由所有子类型实现的虚拟方法。

Shape 算是一切图形或者形状的“基类”了。它可以最灵活的方式表示几何多边形。它是由边缘的集合构成的，可选地定义内部。由给定 Shape 表示的所有几何图形必须具有相同的尺寸，这意味着 Shape 可以表示一组点，一组多边形或一组多边形。 Shape 被定义为一个接口，以便让客户端更加方便的控制底层的数据表示。有时一个 Shape 没有自己的数据，而是包装其他类型的数据。 Shape 操作通常在 ShapeIndex 上定义，而不是单独的形状。 ShapeIndex 只是一个 Shapes 集合，可能有不同的维度（例如10个点和3个多边形），组织成一个数据结构，以便高效的访问。 Shape 的边缘由从 0 开始的连续范围的边缘 ID 索引。边缘被进一步细分为链，其中每个链由端到端连接的一系列边（多段线）组成。例如，表示两条折线 AB 和 CDE 的形状将具有分成两条链（AB）和（CD，DE）的三条边（AB，CD，DE）。类似地，代表5个点的形状将具有由一个边缘组成的5个链。 Shape具有允许使用全局编号（边缘ID）或在特定链中访问边的方法。全局编号对于大多数情况来说是足够的，但链表示对于某些算法（如交集（请参阅BooleanOperation））非常有用。

S2 中总共定义了两个用于表示几何的可扩展接口：S2Shape 和 S2Region。

它们两者不同点是： S2Shape 的目的是灵活地表示多边形几何。 （这不仅包括多边形，还包括点和折线）。大部分的核心 S2 操作将与任何实现 S2Shape 接口的类一起工作。

S2Region 的目的是计算几何的近似值。例如，有计算边界矩形和圆盘的方法，S2RegionCoverer 可以用来逼近一个区域，以任意期望的精度作为单元的集合。与S2Shape 不同，S2Region 可以表示非多边形几何形状，例如球帽（S2Cap）。

除去上面说的这几种常用的类型以外，还有以下这些中间类型或者更加底层的类型可以供开发者使用。

最后，额外说一句，S2RegionTermIndexer 这个类型是支持索引和查询任何类型的S2Region，也就是上述说的所有的类型。可以使用 S2RegionTermIndexer 来索引一组多段线，然后查询哪些多段线与给定的多边形相交。

RegionCoverer 主要是要找到一个能覆盖当前区域的近似最优解(为何不是最优解？)

转换条件主要有3个，MaxLevel, MaxCells, MinLevel。MaxCells 决定了最大 cell 的个数，但是太接近最大值以后又会导致覆盖面积太大，不精确。所以最大个数仅仅是限制在满足最大精度的条件下最多不能超过这个个数。由于这一点导致并不是满足 MaxCells 的最优解。

举几个例子：

下面是一个半径为 10 公里的 cap，并且这个 cap 位于 3 个 face 夹角处，我们假设需要最大个数为 10 的 cell 去覆盖它。结果如下：

相同的设置，我们把个数改到 20，如下：

还是相同的配置，把个数改到 50 个：

到目前，精确度马马虎虎，边缘部门覆盖的还是多于原始的 cap 了，我们继续提高精度，把个数调整到 200 。

200 个看似比较精确了，我们再提高一下，提高到 1000，看看会出现什么结果。

代码配置上虽然设置的 1000 个，实际只有 706 个 cell。原因是代码上虽然是按照 1000 个计算的，但是实际算法处理上还会进行 cell 剪枝后的合并。所以最终个数会小于 1000 个。

再看一个例子，下面这个矩形是表示一个纬度/经度矩形从北纬 60 度延伸到北纬 80 度，从 -170 度经度延伸到 +170 度。覆盖范围限于 8 个单元。请注意，中间的洞被完全覆盖。这明显是不符合我们的意图的。

我们把 cell 的个数提高到 20 个。中间的孔依旧被填补上了。

我们把参数调节到 100 个，其他配置都完全一样。现在中间的孔有一定样子了。但是日期线附近的空白还是没有出来。

最后把参数调整到 500 个。现在中间的孔就比较完整的显示出来了。

在举几个我们实际项目中用到的例子。下面是上海的一个网格的边缘。我们先用

去转换，得到的结果如下:

精确度提高到1000，结果就会如下：

还有一些区域更大的情况，比如一个省，湖北省：

或者一个湖，太湖：

最后在举一个 polygon 的例子。我们知道 polygon 是由多个 loop 组成的：

下面一次是两个 loop 在地图上的样子。

最后是 polygon,它包含了以上2个 loop。

这一章节详细分析一下 Covering 是如何实现的。

最常见的用法就是下面几行：

上面例子展示的是最多转换以后 CellUnion 里面有5个 CellID。上面这三行覆盖的区域是一个 Cap。

RegionCoverer 是一个结构体，实际上里面是包含4个元素的。MinLevel、MaxLevel、MaxCells，这3个应该不用解释了，用的很多。关键要说明说明 LevelMod。

LevelMod 这个值一旦设置了以后，在进行 RegionCover 转换的时候，选取的 Cell Level 只能是 (level - MinLevel) % LevelMod = 0，即 (level - MinLevel) 只能是 LevelMod 的倍数，满足这个条件的 Cell Level 才会被选取。这能有效地允许S2 CellID 层级的分支因子增加。当前的参数取值只能是0，1，2，3，对应的分支因子是0，4，16，64 。

再来谈谈算法的核心思想。

RegionCover 可以被抽象成这样一种问题，给定一个区域，用尽可能精确的 Cell 去覆盖它，但是个数最多不要超过 MaxCells 的个数，问如何去找到这些 Cell ？

这个问题就是一个近视最优解的问题。如果想最精确，方案当然是边缘部分全部都用 MaxLevel 去铺（Level 越大，格子越小）这样就最精确。但是这样会导致 Cell 的个数陡增，远远超过 MaxCells，这样就又不符合要求了。那如何能在 0 && (!c.interiorCovering || len(c.result) < c.maxCells) { cand := heap.Pop(&c.pq).(*candidate) if c.interiorCovering || int(cand.cell.level) < c.minLevel || cand.numChildren == 1 || len(c.result)+c.pq.Len()+cand.numChildren