Google S2 Geometry 浅解

前言

　　Google S2 Geometry（以下简称 S2）是 Google 发明的基于单位球的一种地图投影和空间索引算法，该算法可快速进行覆盖以及邻域计算。更多详见 S2Geometry，Google’s S2, geometry on the sphere, cells and Hilbert curve，halfrost 的空间索引系列文章。虽然使用 S2 已有一年的时间，但确实没有比较系统的看过其源码，这次借着这段空闲时间，将 Shaun 常用的功能系统的看看其具体实现，下文将结合 S2 的 C++，Java，Go 的版本一起看，由于 Java 和 Go 的都算是 C++ 的衍生版，所以以 C++ 为主，捎带写写这三种语言实现上的一些区别，Java 版本时隔 10 年更新了 2.0 版本，喜大普奔。

坐标篇

　　S2 的投影方式可简单想象为一个单位球外接一个立方体，从球心发出一条射线得到球面上的点到立方体上 6 个面的投影，即将球面投影为立方体，当然中间为了使面积分布更为均匀，还做了些其他坐标变换。

S2LatLng 坐标

　　首先是经纬度坐标，默认用弧度（Radians）构造，取值范围为经度 [-π，+π]，纬度 [-π/2，+π/2]，当然也可使用 S1Angle 将角度（Degrees）转成弧度来构造。

S2Point 坐标

　　然后球面笛卡尔坐标，这是个三维坐标，由 S2LatLng 到 S2Point 相当于将单位球的极坐标表示法转换为笛卡尔坐标表示法，具体公式为 \(x=\cos(lat)cos(lng);　y=cos(lat)sin(lng);　z=sin(lat)\)。

FaceUV 坐标

　　这个坐标并没实际的类与其对应，face 指的是立方体的面，值域为 [0,5]，而 uv 坐标是指面上的点，值域为 [-1,1]。首先需要知道 S2Point 会投影到哪个面上，可以知道 S2 的笛卡尔坐标 X 轴正向指向 0 面，Y 轴正向指向 1 面，Z 轴正向指向 2 面，X 轴负向指向 3 面，Y 轴负向指向 4 面，Z 轴负向指向 5 面，所以 S2Point xyz 哪个分量的绝对值最大，就会投影到哪个轴指向的面，若该分量为正值，则取正向指的面，若该分量为负值，则取负向指的面。至于 uv 的计算方式就是直线与平面的交点了，之前的一篇「计算几何基础」中写过，但这里的平面和直线都比较特殊，所以有快速算法，就直接贴 Go 的代码吧：

// validFaceXYZToUV given a valid face for the given point r (meaning that
// dot product of r with the face normal is positive), returns
// the corresponding u and v values, which may lie outside the range [-1,1].
func validFaceXYZToUV(face int, r r3.Vector) (float64, float64) {
	switch face {
	case 0:
		return r.Y / r.X, r.Z / r.X
	case 1:
		return -r.X / r.Y, r.Z / r.Y
	case 2:
		return -r.X / r.Z, -r.Y / r.Z
	case 3:
		return r.Z / r.X, r.Y / r.X
	case 4:
		return r.Z / r.Y, -r.X / r.Y
	}
	return -r.Y / r.Z, -r.X / r.Z
}

　　这里需要注意的是 S2Point xyz 三分量构成的向量与平面法向量的点积必须是正数时 uv 才算正确有效，Go 在计算时没做校验，C++ 和 Java 都有校验，使用时需要注意。

FaceST 坐标

　　之所以引入 ST 坐标是因为同样的球面面积映射到 UV 坐标面积大小不一，大小差距比较大（离坐标轴越近越小，越远越大），所以再做一次 ST 变换，将面积大的变小，小的变大，使面积更均匀，利于后面在立方体面上取均匀格网（cell）时，每个 cell 对应球面面积差距不大。S2 的 ST 变换有三种：1、线性变换，基本没做任何变形，只是简单将 ST 坐标的值域变换为 [0, 1]，cell 对应面积最大与最小比大约为 5.2；2、二次变换，一种非线性变换，能起到使 ST 空间面积更均匀的作用，cell 对应面积最大与最小比大约为 2.1；3、正切变换，同样能使 ST 空间面积更均匀，且 cell 对应面积最大与最小比大约为 1.4，不过其计算速度相较于二次变换要慢 3 倍，所以 S2 权衡考虑，最终采用了二次变换作为默认的 UV 到 ST 之间的变换。二次变换公式为：

public double stToUV(double s) {
  if (s >= 0.5) {
    return (1 / 3.) * (4 * s * s - 1);
  } else {
    return (1 / 3.) * (1 - 4 * (1 - s) * (1 - s));
  }
}

public double uvToST(double u) {
  if (u >= 0) {
    return 0.5 * Math.sqrt(1 + 3 * u);
  } else {
    return 1 - 0.5 * Math.sqrt(1 - 3 * u);
  }
}

FaceIJ 坐标

　　IJ 坐标是离散化后的 ST 坐标，将 ST 空间的平面划分为 \(2^{30}×2^{30}\) 个网格，取网格所在的横纵坐标得到 IJ 坐标，所以由 ST 到 IJ 坐标的变换就比较简单了：

public static int stToIj(double s) {
  return Math.max(
    0, Math.min(1073741824 - 1, (int) Math.round(1073741824 * s - 0.5))
  );
}

S2CellId

　　这个 id 其实是个一维坐标，而是利用希尔伯特空间填充曲线将 IJ 坐标从二维变换为一维，该 id 用一个 64 位整型表示，高 3 位用来表示 face（0~5），后面 61 位来保存不同的 level（0~30）对应的希尔伯特曲线位置，每增加一个 level 增加两位，后面紧跟一个 1，最后的位数都补 0。注：Java 版本的 id 是有符号 64 位整型，而 C++ 和 Go 的是无符号 64 位整型，所以在跨语言传递 id 的时候，在南极洲所属的最后一个面（即 face = 5）需要小心处理。

HilbertCurve

　　上面两张图很明了的展示了希尔伯特曲线的构造过程，该曲线的构造基本元素由 ABCD 4 种“U”形构成，而 BCD 又可由 A 依次逆时针旋转 90 度得到，所以也可以认为只有一种“U”形，每个 U 占 4 个格子，以特定方式进行 1 分 4 得到下一阶曲线形状。

每个 U 坐标与希尔伯特位置（用二进制表示）对应关系如下：

A：00 -> (0,0); 01 -> (0,1); 10 -> (1,1); 11 -> (1,0);
B：00 -> (1,1); 01 -> (0,1); 10 -> (0,0); 11 -> (1,0);
C：00 -> (1,1); 01 -> (1,0); 10 -> (0,0); 11 -> (0,1);
D：00 -> (0,0); 01 -> (1,0); 10 -> (1,1); 11 -> (0,1);

每个 U 一分四对应关系如下：

A：D -> A -> A -> B
B：C -> B -> B -> A
C：B -> C -> C -> D
D：A -> D -> D -> C

　　根据以上两个对应关系就能找到右手坐标系任意阶数的希尔伯特位置及坐标对应关系。以初始 1 阶曲线 A 为例，占据四个格子，然后进行一分四操作，四个格子分成 16 个格子，A 分为 DAAB 四个“U”形，连接起来即为 2 阶曲线，位置与坐标对应关系为（都用二进制表示）：

0000 -> (00, 00); 0001 -> (01, 00); 0010 -> (01, 01); 0011 -> (00, 01)；

0100 -> (00, 10); 0101 -> (00, 11); 0110 -> (01, 11); 0111 -> (01, 10)；

1000 -> (10, 10); 1001 -> (10, 11); 1010 -> (11, 11); 1011 -> (11, 10)；

1100 -> (11, 01); 1101 -> (10, 01); 1110 -> (10, 00); 1111 -> (11, 00)；

　　从二进制中很容易看出随着阶数的增加，位置与坐标的对应关系：每增加一阶，位置往后增加两位，坐标分量各增加一位，位置增加的两位根据一分四对应关系拼接，坐标各分量增加的一位需先找到一分四对应关系，再找对应位置与坐标对应关系，将得到的坐标分量对应拼接。以一阶的 01 -> (0,1) 到二阶的 0110 -> (01, 11) 为例，首先根据 01 得到当前所属一阶第二块，查找一分四对应关系知道，下一阶这块还是 A，根据 0110 后两位 10 可知这块属于 A 的第三个位置，查找坐标得到是 (1,1)，结合一阶的 (0,1)，对应分量拼接得到坐标 (01,11)，即 (1, 3)，同理可根据第二阶的坐标反查第二阶的位置。有了这些关系，就能生成希尔伯特曲线了，下面就看看 S2 是怎么生成 id 的。

S2Id

　　首先 S2 中用了两个二维数组分别保存位置到坐标以及坐标到位置的对应的关系：

// kIJtoPos[orientation][ij] -> pos
const int kIJtoPos[4][4] = {
  // (0,0) (0,1) (1,0) (1,1)
  {     0,    1,    3,    2  },  // canonical order
  {     0,    3,    1,    2  },  // axes swapped
  {     2,    3,    1,    0  },  // bits inverted
  {     2,    1,    3,    0  },  // swapped & inverted
};

// kPosToIJ[orientation][pos] -> ij
const int kPosToIJ[4][4] = {
  // 0  1  2  3
  {  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)
};

// kPosToOrientation[pos] -> orientation_modifier
const int kPosToOrientation[4] = {1, 0, 0, 3};

　　方向 0（canonical order）相当于上文中 A，方向 1（axes swapped）相当于上文中 D，方向 2（bits inverted）相当于上文中 C，方向 3（swapped & inverted）相当于上文中 B，kPosToOrientation 代表 S2 中方向 0 一分四的对应关系，而方向 1，2，3 的对应关系可由该值推出，计算公式为 orientation ^ kPosToOrientation，eg：1 -> 1^kPosToOrientation=[0, 1, 1, 2]; 3 -> 3^kPosToOrientation=[2, 3, 3, 0]，与上文中一分四对应关系一致。

　　随后 S2 初始化了一个 4 阶希尔伯特曲线位置与坐标的对应关系查找表，见 C++ 版的 MaybeInit() 方法，

int ij = (i << 4) + j;
lookup_pos[(ij << 2) + orig_orientation] = (pos << 2) + orientation;
lookup_ij[(pos << 2) + orig_orientation] = (ij << 2) + orientation;

　　orig_orientation 代表 4 个初始方向，orientation 代表该位置或坐标下一阶一分四的方向，数组中每个元素是 16 位数，2 个字节，一个四阶希尔伯特曲线是 \(2^4×2^4=256\) 个位置，一个初始方向对应一个四阶希尔伯特曲线，所以一个查找表共占内存 \(2×256×4=2048=2KB\)，正好一级缓存能放下，再大的话，一级缓存可能放不下，反而会降低查找速度。这两个查找表就相当于 4 个超“U”形的位置与坐标对应关系，同时一分四对应关系保持不变，以超“U”作为基本元素做下一阶希尔伯特曲线，每增加一阶位置往后增加 8 位，IJ 坐标各往后增加 4 位，如此，以更快的速度迭代到 S2 想要的 30 阶希尔伯特曲线。C++ 的这份代码就很精妙了：

S2CellId S2CellId::FromFaceIJ(int face, int i, int j) {
  // 初始化超“U”形查找表
  MaybeInit();

  // face 向左移 60 位
  uint64 n = absl::implicit_cast<uint64>(face) << (kPosBits - 1);

  // 确定每个面的初始“U”形方向，使每个面都保持相同的右手坐标系，6 个面生成的希尔伯特曲线可以依次相连
  uint64 bits = (face & kSwapMask);

  // 基于超“U”形得到 30 阶希尔伯特曲线 IJ 坐标对应位置
#define GET_BITS(k) do { \
    const int mask = (1 << kLookupBits) - 1; \
    bits += ((i >> (k * kLookupBits)) & mask) << (kLookupBits + 2); \
    bits += ((j >> (k * kLookupBits)) & mask) << 2; \
    bits = lookup_pos[bits]; \
    n |= (bits >> 2) << (k * 2 * kLookupBits); \
    bits &= (kSwapMask | kInvertMask); \
  } while (0)

  // IJ 只有 30 位，7 这个调用只会导致位置移 4 位，后续调用都移 8 位，得到 4 + 8 * 7 = 60 位
  GET_BITS(7); 
  GET_BITS(6);
  GET_BITS(5);
  GET_BITS(4);
  GET_BITS(3);
  GET_BITS(2);
  GET_BITS(1);
  GET_BITS(0);
#undef GET_BITS

  // 整个 n 向右移一位，再以 1 结尾
  return S2CellId(n * 2 + 1);
}

再来看看根据 id 反算 IJ 坐标：

int S2CellId::ToFaceIJOrientation(int* pi, int* pj, int* orientation) const {
  // 与上面一样
  MaybeInit();

  int i = 0, j = 0;
  int face = this->face();
  int bits = (face & kSwapMask);

  // 反算 IJ 坐标，k == 7 时，取希尔伯特曲线位置高 4 位，IJ 各前 2 位，其余依次取位置 8 位， IJ 各 4 位
#define GET_BITS(k) do { \
    const int nbits = (k == 7) ? (kMaxLevel - 7 * kLookupBits) : kLookupBits; \
    bits += (static_cast<int>(id_ >> (k * 2 * kLookupBits + 1)) \
             & ((1 << (2 * nbits)) - 1)) << 2; \
    bits = lookup_ij[bits]; \
    i += (bits >> (kLookupBits + 2)) << (k * kLookupBits); \
    j += ((bits >> 2) & ((1 << kLookupBits) - 1)) << (k * kLookupBits); \
    bits &= (kSwapMask | kInvertMask); \
  } while (0)

  GET_BITS(7);
  GET_BITS(6);
  GET_BITS(5);
  GET_BITS(4);
  GET_BITS(3);
  GET_BITS(2);
  GET_BITS(1);
  GET_BITS(0);
#undef GET_BITS

  *pi = i;
  *pj = j;

  if (orientation != nullptr) {
    S2_DCHECK_EQ(0, kPosToOrientation[2]);
    S2_DCHECK_EQ(kSwapMask, kPosToOrientation[0]);
    // 0x1111111111111111ULL may be better?
    if (lsb() & 0x1111111111111110ULL) {
      bits ^= kSwapMask;
    }
    *orientation = bits;
  }
  return face;
}

　　这里的 orientation 实际是指当前位置的方向，即其周围必有 3 个位置与其方向相同，最后一行注释 Shaun 之所以认为应该是 0x1111111111111111ULL，是因为第 30 阶希尔伯特曲线位置（leaf cell）按理说同样需要做异或操作得到方向，不过整个 S2 库都没有需要用到 leaf cell 的方向，所以这就倒无关紧要了。之所以需要做异或操作，是因为 bits 是该位置下一阶一分四的方向，而对于同一个希尔伯特曲线位置，奇数阶与奇数阶下一阶一分四方向相同，偶数阶与偶数阶下一阶一分四方向相同，lsb() 表示二进制 id 从右往左数第一个 1 所代表的数，所以有 0x1111111111111110ULL 这一魔术数，而异或操作正好能将下一阶一分四方向调整为当前阶方向。

　　如此 S2 的坐标以及 id 的生成以及反算就很明了了，下面就是 S2 如何使用 id 做计算了。

FaceSiTi 坐标

　　这个是 S2 内部计算使用的坐标，一般用来计算 cell 的中心坐标，以及根据当前 s 和 t 坐标的精度（小数点后几位）判断对应的级别（level）。由于 S2 本身并不显式存储 ST 坐标（有存 UV 坐标），所以 ST 坐标只能计算出来，每个 cell 的中心点同样如此。计算公式为 \(Si=s*2^{31};Ti=t*2^{31}\)。至于为啥是 \(2^{31}\)，是因为该坐标是用来描述从 0~ 31 阶希尔伯特曲线网格的中心坐标，0 阶中心以 \(1/2^1\) 递增，1 阶中心以 \(1/2^2\) 递增，2 阶中心以 \(1/2^3\) 递增，……，30 阶中心以 \(1/2^{31}\) 递增。S2 计算 id 对应的格子中心坐标，首先就会计算 SiTi 坐标，再将 SiTi 转成 ST 坐标。

算法篇

邻域算法

　　S2 计算邻域，最关键的是计算不同面相邻的 leaf cell id，即：

S2CellId S2CellId::FromFaceIJWrap(int face, int i, int j) {
  // 限制 IJ 最大最小取值为 -1~2^30, 刚好能超出 IJ 正常表示范围 0~2^30-1
  i = max(-1, min(kMaxSize, i));
  j = max(-1, min(kMaxSize, j));

  static const double kScale = 1.0 / kMaxSize;
  static const double kLimit = 1.0 + DBL_EPSILON;
  S2_DCHECK_EQ(0, kMaxSize % 2);
  // IJ -> SiTi -> ST -> UV
  double u = max(-kLimit, min(kLimit, kScale * (2 * (i - kMaxSize / 2) + 1)));
  double v = max(-kLimit, min(kLimit, kScale * (2 * (j - kMaxSize / 2) + 1)));

  face = S2::XYZtoFaceUV(S2::FaceUVtoXYZ(face, u, v), &u, &v);
  return FromFaceIJ(face, S2::STtoIJ(0.5*(u+1)), S2::STtoIJ(0.5*(v+1)));
}

　　这个算法主要用来计算超出范围（0~2^30-1）的 IJ 对应的 id，核心思想是先将 FaceIJ 转为 XYZ，再使用 XYZ 反算得到正常的 FaceIJ，进而得到正常的 id。中间 IJ -> UV 中坐标实际经过了 3 步，对于 leaf cell，IJ -> SiTi 的公式为 \(Si=2×I+1\)，而对于 ST -> UV，这里没有采用二次变换，就是线性变换 \(u=2*s-1\)，官方注释上说明用哪个变换效果都一样，所以采用最简单的就行。

边邻域

　　边邻域代码很简单，也很好理解：

void S2CellId::GetEdgeNeighbors(S2CellId neighbors[4]) const {
  int i, j;
  int level = this->level();
  // 计算当前 level 一行或一列对应多少个 30 级的 cell（leaf cell） 2^(30-level)
  int size = GetSizeIJ(level);
  int face = ToFaceIJOrientation(&i, &j, nullptr);

  // Edges 0, 1, 2, 3 are in the down, right, up, left directions.
  neighbors[0] = FromFaceIJSame(face, i, j - size, j - size >= 0)
                 .parent(level);
  neighbors[1] = FromFaceIJSame(face, i + size, j, i + size < kMaxSize)
                 .parent(level);
  neighbors[2] = FromFaceIJSame(face, i, j + size, j + size < kMaxSize)
                 .parent(level);
  neighbors[3] = FromFaceIJSame(face, i - size, j, i - size >= 0)
                 .parent(level);
}

　　分别计算当前 IJ 坐标下右上左坐标对应 id，FromFaceIJSame 表示若邻域在相同面，则走 FromFaceIJ，否则走 FromFaceIJWrap，由于这两个函数得到都是 leaf cell，要上升到指定 level，需要用到 parent 方法，即将希尔伯特曲线位置去掉右 \(2*(30-level)\) 位，再组合成新的 id，位运算也很有意思：

static uint64 lsb_for_level(int level) {
  return uint64{1} << (2 * (kMaxLevel - level));
}

inline S2CellId S2CellId::parent(int level) const {
  uint64 new_lsb = lsb_for_level(level);
  // 取反加一实际是取负数
  return S2CellId((id_ & (~new_lsb + 1)) | new_lsb);
}

点邻域

　　S2 的点邻域并不是指常规意义上 4 个顶点相邻左上右上右下左下的 id，而是一种比较特殊的相邻关系，以直角坐标系 (0,0),(0,1),(1,1),(1,0) 为例，(0,0) 的点邻域为 (0,0),(0,-1),(-1,-1),(-1,0)，(0,1) 的点邻域为 (0,1),(0,2),(-1,2),(-1,1)，(1,1) 的点邻域为 (1,1),(1,2),(2,2),(2,1)，(1,0) 的点邻域为 (1,0),(1,-1),(2,-1),(2,0)。具体代码如下：

void S2CellId::AppendVertexNeighbors(int level,
                                     vector<S2CellId>* output) const {
  // level < this->level()
  S2_DCHECK_LT(level, this->level());
  int i, j;
  int face = ToFaceIJOrientation(&i, &j, nullptr);

  // 判断 IJ 落在 level 对应 cell 的哪个方位？（左下左上右上右下，对应上文的(0,0),(0,1),(1,1),(1,0)坐标）
  int halfsize = GetSizeIJ(level + 1);
  int size = halfsize << 1;
  bool isame, jsame;
  int ioffset, joffset;
  if (i & halfsize) {
    ioffset = size;
    isame = (i + size) < kMaxSize;
  } else {
    ioffset = -size;
    isame = (i - size) >= 0;
  }
  if (j & halfsize) {
    joffset = size;
    jsame = (j + size) < kMaxSize;
  } else {
    joffset = -size;
    jsame = (j - size) >= 0;
  }

  output->push_back(parent(level));
  output->push_back(FromFaceIJSame(face, i + ioffset, j, isame).parent(level));
  output->push_back(FromFaceIJSame(face, i, j + joffset, jsame).parent(level));
  // 则邻域的 IJ 与当前 cell 都不在同一个面，则说明只有三个点邻域
  if (isame || jsame) {
    output->push_back(FromFaceIJSame(face, i + ioffset, j + joffset,
                                     isame && jsame).parent(level));
  }
}

　　上面的代码算是比较清晰了，3 个点邻域的情况一般出现在当前 id 位于立方体 6 个面的角落，该方法的参数 level 必须比当前 id 的 level 要小。

全邻域

　　所谓全邻域，即为当前 id 对应 cell 周围一圈 cell 对应的 id，若周围一圈 cell 的 level 与当前 id 的 level 一样，则所求即为正常的 9 邻域。具体代码如下：

void S2CellId::AppendAllNeighbors(int nbr_level,
                                  vector<S2CellId>* output) const {
  // nbr_level >= level
  S2_DCHECK_GE(nbr_level, level());
  int i, j;
  int face = ToFaceIJOrientation(&i, &j, nullptr);

  // 先归一 IJ 坐标，将 IJ 坐标调整为当前 cell 左下角 leaf cell 的坐标
  int size = GetSizeIJ();
  i &= -size;
  j &= -size;

  int nbr_size = GetSizeIJ(nbr_level);
  S2_DCHECK_LE(nbr_size, size);

  for (int k = -nbr_size; ; k += nbr_size) {
    bool same_face;
    if (k < 0) {
      same_face = (j + k >= 0);
    } else if (k >= size) {
      same_face = (j + k < kMaxSize);
    } else {
      same_face = true;
      // 生成外包围圈下上两边的 id, 顺序为从左往右
      output->push_back(FromFaceIJSame(face, i + k, j - nbr_size,
                                       j - size >= 0).parent(nbr_level));
      output->push_back(FromFaceIJSame(face, i + k, j + size,
                                       j + size < kMaxSize).parent(nbr_level));
    }
    // 生成外包围圈左右两边以及四个边角的 id, 顺序为从下往上
    output->push_back(FromFaceIJSame(face, i - nbr_size, j + k,
                                     same_face && i - size >= 0)
                      .parent(nbr_level));
    output->push_back(FromFaceIJSame(face, i + size, j + k,
                                     same_face && i + size < kMaxSize)
                      .parent(nbr_level));
    if (k >= size) break;
  }
}

　　知道这个函数的作用，再看代码就很明了了，这个方法的参数 nbr_level 必须大于或等于当前 id 的 level，因为一旦外包围圈的 cell 面积比当前 cell 还大，就无法得到正确的外包围圈。

覆盖算法

　　S2 的覆盖，是指给定一块区域，能用多少 id 对应的 cell 完全覆盖该区域（GetCovering），当然也有尽量覆盖的算法（GetInteriorCovering），下面主要解析 GetCovering，因为 GetInteriorCovering 也差不多，就是覆盖策略略有不同。

GetCovering 的区域入参是 S2Region，比较典型的 S2Region 有以下几种：

S2Cell：S2 id 对应的网格，会保存左下右上两个 UV 坐标，也是覆盖算法使用的基本元素；
S2CellUnion：多个 S2Cell 集合体，GetCovering 的返回值；
S2LatLngRect：经纬度矩形区域；
S2Cap：球帽区域，类比于二维圆的圆弧，球帽的构造比较奇怪，球帽的中心 S2Point 是需要，但另一个变量不是球帽的圆弧角，而是半个圆弧角（S2 代码库对应的 S1Angle 弧度，90 度代表半球，180 度代表全球）所对应弦长的平方，最大值为 4，之所以采用弦长的平方作为默认构造，是因为这就是 3 维中距离，在进行距离比较的场景时会更方便，比如测试是否包含一个 S2Point，计算覆盖多边形时，就不用再比较角度，毕竟角度计算代价比较大；
S2Loop：多边形的基本组成元素，第一个点与最后一个点隐式连接，逆时针代表封闭，顺时针代表开孔取外围区域，不允许自相交；
S2Polygon：非常正常的复杂多边形，由多个 S2Loop 构成，S2Loop 之间不能相交；
S2Polyline：一条折线，同样不能自相交；
还有些其它不常用的：S2R2Rect（S2Point 矩形区域），S2RegionIntersection（集合相交区域），S2RegionUnion（集合合并区域），……等。

　　S2 覆盖算法的本质是一种启发式算法，先取满足当前条件最基本的元素，再依照条件进行迭代优化，所以该算法得到的只是一个近似最优解。GetCovering 需要依次满足以下条件：

生成的 S2Cell level 不能比指定的 minLevel 小；（必须满足）
生成的 S2Cell 的个数不能比指定的 maxCells 多；（可以满足，当满足 1 时，数目已经 maxCells 多，迭代停止）
生成的 S2Cell level 不能比指定的 maxLevel 大；（必须满足）

　　以上 3 个条件对应 GetCovering 的其他三个参数，当然还有一个参数是 levelModel，表示从 minLevel 向下分到 maxLevel 时，是 1 分 4，还是 1 分 16，还是 1 分 64，对应一次升 1 阶曲线，还是一次升 2 阶，或是一次升 3 阶。下面就来具体看看 GetCovering 的算法流程（代码就不贴了，太多了）：

首先获取候选种子 S2Cell。先构造一个临时覆盖器，设置 maxCells 为 4，minLevel 为 0，以快速得到初始覆盖结果，做法为：先得到覆盖输入区域的 S2Cap，再用 S2CellUnion 覆盖该 S2Cap，根据 S2Cap 圆弧度计算 S2Cell 的 level，若最终 level < 0，则说明 S2Cap 非常大，需要取 6 个面对应的 S2Cell，否则只需要取 S2Cap 中心点对应 S2Cell 的 level 级的点邻域 4 个 S2Cell 作为初始候选 S2Cell。
然后标准化候选种子。第一步，如果候选 S2Cell level 比 maxLevel 大或者候选 S2Cell 的 level 不符合 levelModel，则调整候选 S2Cell 的 level，用指定父级 S2Cell 来代替；第二步，归一化候选 S2Cell，先对 S2Cell 按 id 排序，去除被包含的 id，以及对 id 剪枝（若连续 4 个 S2Cell 共有同一个 parent，则用 parent 代替这 4 个 S2Cell）；第三步，反归一化候选 S2Cell，若候选 S2Cell level 比 minLevel 小或不满足 levelModel，则需要将 S2Cell 分裂，用指定级别的孩子来取代该 S2Cell；第四步，检查是否满足全部条件，若满足，则标准化完成，若不满足，则看候选 S2Cell 的数目是否足够多，若足够多，则需要迭代进行 GetCovering，这样会极大降低算法性能，若不是很多，则迭代合并相同祖先的两个 S2Cell（当然祖先的 level 不能比 minLevel 小），最后再次检查所有候选 S2Cell 是否达到标准化要求，并调整 S2Cell level。
构造优先级队列。将符合条件（与入参区域相交）的候选 S2Cell 放进一个优先级队列中，优先级会依次根据三个参数进行判断，1、S2Cell 的大小（level 越大，S2Cell 越小），越大的优先级越高；2、入参区域与候选 S2Cell 孩子相交（这里的相交是指相交但不完全包含）的个数，越少优先级越高；3、入参区域完全包含候选 S2Cell 孩子和与无法再细分的孩子的个数，同样是越少优先级越高。在构造这个优先级队列的同时，会输出一些候选 S2Cell 作为覆盖算法的正式结果，这些 S2Cell 满足任意以下条件：1、被入参区域完全覆盖；2、与入参区域相交但不可再细分；3、入参区域包含或相交全部孩子。如此留在优先级队列中的，就都是些与入参区域边界相交的 S2Cell，这些就是真正的候选 S2Cell。
最后，处理优先级队列中的 S2Cell。处理方式也比较简单粗暴，继续细分并入队，满足上面3个出队条件的任意一个，即可出队作为正式结果，当然，若分到后面可能正式的 S2Cell 太多，甚至超过 maxCells，这时不再细分强行出队作为正式结果。最后，再对正式结果做一次标准化处理，即进行第 2 步，得到最终的覆盖结果。

　　以上就是 S2 覆盖算法的大致流程，更加细节的东西，还是得看代码，文字有些不是很好描述，代码里面计算候选 S2Cell 的优先级就很有意思。

　　当然 S2 中还有很多其他算法（凸包，相交，距离），这里就不做太多介绍了，Shaun 平常用的最多的就是覆盖算法，之前一直没有细看，就简单用用 api，同时为了对一块大的 S2Cell 做多线程处理，需要了解 S2Cell 一分四的方向，经过这次对 S2 的了解，发现之前的用法存在一些问题，可见调包侠同样需要对包有一定的了解才能调好包 ╮(╯▽╰)╭。

后记

　　正如许多经典的算法一样，看完之后总有种我上我也行的感觉，但实际完全不行，S2 全程看下来有些地方确实比较晦涩，而且这一连串的想法也很精妙（单位球立方体投影，ST 空间面积优化，64 位 id 生成等），Shaun 或许能有部分想法，但这么多奇思妙想组合起来，就完全不行。

附录

HilbertCurve 绘制

　　在网上随便找了三种实现方式，并用 threejs 简单绘制了一下：

import * as THREE from "three";

export default class HilbertCurve {
    order = 3; // 阶数
    size = 1 << this.order; // 行列数
    totalSize = this.size * this.size; // 总网格数，希尔伯特长度

    // https://www.youtube.com/watch?v=dSK-MW-zuAc
    getPath_V1() {
        let path = [];
        let origOrientation = [
            [0, 0],
            [0, 1],
            [1, 1],
            [1, 0],
        ]; // 倒 U 形

        for (let i = 0; i < this.totalSize; i++) {
            path.push(hilbertToXY(i, this.order));
        }

        return path;

        function hilbertToXY(i: number, order: number) {
            let index = i & 3;
            let curCoord = origOrientation[index].slice();

            for (let ord = 1; ord < order; ord++) {
                i = i >>> 2;
                index = i & 3;
                let delta = 1 << ord;
                if (index === 0) {
                    // 顺时针旋转 90°
                    let tmp = curCoord[0];
                    curCoord[0] = curCoord[1];
                    curCoord[1] = tmp;
                } else if (index === 1) {
                    curCoord[1] += delta;
                } else if (index === 2) {
                    curCoord[0] += delta;
                    curCoord[1] += delta;
                } else if (index === 3) {
                    // 逆时针旋转 90°
                    let tmp = delta - 1 - curCoord[0];
                    curCoord[0] = delta - 1 - curCoord[1];
                    curCoord[1] = tmp;
                    curCoord[0] += delta;
                }
            }

            return curCoord;
        }
    }

    // Hacker's Delight
    getPath_V2() {
        let path: number[][] = [];
        let x = -1,
            y = 0;
        let s = 0; // along the curve

        step(0);
        hilbert(0, 1, this.order);

        return path;

        function step(dir: number) {
            switch (dir & 3) {
                case 0:
                    x += 1;
                    break;
                case 1:
                    y += 1;
                    break;
                case 2:
                    x -= 1;
                    break;
                case 3:
                    y -= 1;
                    break;
            }

            path[s] = [x, y];

            s += 1;
        }

        function hilbert(dir: number, rot: number, order: number) {
            if (order === 0) return;

            dir += rot;
            hilbert(dir, -rot, order - 1);
            step(dir);

            dir -= rot;
            hilbert(dir, rot, order - 1);
            step(dir);

            hilbert(dir, rot, order - 1);

            dir -= rot;
            step(dir);
            hilbert(dir, -rot, order - 1);
        }
    }

    // https://en.wikipedia.org/wiki/Hilbert_curve
    getPath_V3() {
        let path: number[][] = [];

        // for (let i = 0; i < this.totalSize; i++) {
        //     path.push(hilbertToXY(this.size, i));
        // }

        for (let y = 0; y < this.size; y++) {
            for (let x = 0; x < this.size; x++) {
                path[xyToHilbert(this.size, x, y)] = [x, y];
            }
        }

        return path;

        function rot(N: number, rx: number, ry: number, xy: number[]) {
            if (ry === 0) {
                if (rx === 1) {
                    xy[0] = N - 1 - xy[0];
                    xy[1] = N - 1 - xy[1];
                }

                let t = xy[0];
                xy[0] = xy[1];
                xy[1] = t;
            }
        }

        function hilbertToXY(N: number, h: number) {
            let t = h;
            let xy = [0, 0];
            for (let s = 1; s < N; s *= 2) {
                let rx = 1 & (t / 2);
                let ry = 1 & (t ^ rx);
                rot(s, rx, ry, xy);
                xy[0] += s * rx;
                xy[1] += s * ry;
                t /= 4;
            }

            return xy;
        }

        function xyToHilbert(N: number, x: number, y: number) {
            let h = 0;
            let xy = [x, y];
            for (let s = N / 2; s > 0; s /= 2) {
                let rx = (xy[0] & s) > 0 ? 1 : 0;
                let ry = (xy[1] & s) > 0 ? 1 : 0;
                h += s * s * ((3 * rx) ^ ry);
                rot(N, rx, ry, xy);
            }

            return h;
        }
    }

    draw() {
        let lineGeometry = new THREE.Geometry();
        this.getPath_V3().forEach((vertice) => {
            let vecot = new THREE.Vector3().fromArray(vertice);
            vecot.setZ(0);
            lineGeometry.vertices.push(vecot);
        });
        let lineMaterial = new THREE.LineBasicMaterial({ color: 0x00ffff, linewidth: 1 });
        let line = new THREE.Line(lineGeometry, lineMaterial);

        return line;
    }
}