Google S2 是如何解決空間覆蓋最優解問題的?

前言

這篇不出意外就是 Google S2 整個系列的最終篇了。這篇裡面會把 regionCoverer 演算法都講解清楚。至於 Google S2 庫裡面還有很多其他的小演算法，程式碼同樣也很值得閱讀和學習，這裡也就不一一展開了，有興趣的讀者可以把整個庫都讀一遍。

在寫這篇文章的同時，發現了庫的作者的一些新的 commit ，比如 2017年12月4號的 commit f9610db2b871b54b17d36d4da6a4d6a2aab6018d，這次提交的改動更改了 README，別看只改了文件，其實裡面內容很多。

 -For an analogous library in C++, see
 -https://code.google.com/archive/p/s2-geometry-library/, and in Java, see
 -https://github.com/google/s2-geometry-library-java
------------------------------------------------------------------------------------
 +For an analogous library in C++, see https://github.com/google/s2geometry, in
 +Java, see https://github.com/google/s2-geometry-library-java, and Python, see
 +https://github.com/google/s2geometry/tree/master/src/python

複製程式碼

可以看到他們把原來存在 Google 官方私有程式碼倉庫裡面的程式碼放到了 github。之前都只能在程式碼歸檔裡面檢視 C++ 程式碼，現在直接可以在 github 上檢視了。方便了很多。

+More details about S2 in general are available on the S2 Geometry Website
 +[s2geometry.io](https://s2geometry.io/).

複製程式碼

在這次 commit 裡面還提到了一個新的網站，這個網站我發現也是最近才釋出出來的。因為筆者我連續關注 S2 每個 commit 了快半年了。網上每個關於 S2 的資源都有關注，這個網站非常新。關於這個網站文章最後會提到，這裡就不詳細說了。

從這個提交開始，筆者認為 Google S2 可能被重視起來了，也有可能打算大力推廣了。

好了，正式進入正題。

一. 空間型別

在 Google S2 中能進行 RegionCovering 的有以下幾種型別。基本上都是必須滿足 Region interface 的。

1. Cap 球帽

Cap 代表由中心和半徑限定的盤形區域。從技術上講，這種形狀被稱為“球形帽”（而不是圓盤），因為它不是平面的。帽子代表被飛機切斷的球體的一部分。帽的邊界是由球面與平面的交點所定義的圓。帽子是一個封閉的組合，即它包含了它的邊界。 大多數情況下，無論在平面幾何中使用光碟，都可以使用球冠。帽的半徑是沿著球體表面測量的（而不是通過內部的直線距離）。因此，一個半徑為 π/ 2 的帽是一個半球，半徑為 π 的帽覆蓋整個球。 中心是單位球面上的一個點。（因此需要它是單位長度）帽子也可以由其中心點和高度來定義。高度是從中心點到截斷平面的距離。還有支援“空”和“全”的上限，分別不包含任何點數和所有點數。 下面是帽高（h），帽半徑（r），帽中心的最大弦長（d）和帽底部半徑（a）之間的一些有用關係。

h = 1 - cos(r)
	= 2 * sin^2(r/2)
d^2 = 2 * h
	= a^2 + h^2

複製程式碼

2. Loop 迴圈

Loop 代表一個簡單的球面多邊形。它由一系列頂點組成，其中第一個頂點隱含地被認為是連線到最後一個頂點的。所有的 loop 被定義為具有 CCW 方向，即 loop 的內部在邊的左側。這意味著包圍一個小區域的順時針 loop 被解釋為包圍非常大的區域的 CCW 的 loop。 loop 不允許有任何重複的頂點（不管是否相鄰）。不允許相鄰的邊相交，而且不允許長度為180度的邊（即，相鄰的頂點不能是相反的）。loop 必須至少有3個頂點（下面討論的“空”和“全” loop 除外）。 有兩個特殊的 loop：EmptyLoop 不包含點，FullLoop 包含所有點。這些 loop 沒有任何邊，但為了保持每一個 loop 都可以表示為頂點鏈的不變數，它們被定義為每個只有一個頂點。

3. Polygon 多邊形

多邊形表示一個零或多個 loop 的序列;同樣，一個 loop 的左手邊方向定義為它的內部。 當多邊形初始化時，給定的 loop 自動轉換為“孔”的組成的規範形式。loop 被重新排序以對應於巢狀層次的預定義遍歷方式。 多邊形可以表示具有多邊形邊界的球體的任何區域，包括整個球體（稱為“完整”多邊形）。完整的多邊形由一個完整的 loop 組成，而空的多邊形完全沒有 loop。 使用 FullPolygon() 來構造一個完整的多邊形。 Polygon 的零值被視為空的多邊形。

想要 多個 loop 構成一個 Polygon 多邊形，必須滿足以下4個條件：

1. loop 不能交叉，即 loop 的邊界可能不與任何其他 loop 的內部和外部相交。
2. loop 不共享邊緣，即如果 loop 包含邊緣 AB，則其他 loop 可能不包含 AB 或 BA。
3. loop 可以共享頂點，但是在單個 loop 中不會出現兩次頂點（參見S2Loop）。
4. 不能有空的 loop。full loop 可能只出現在完整的 full polygon 中。

4. Rect 矩形

Rect 代表一個封閉的經緯度矩形。它也是 Region 型別。它能夠表示空的和完整的矩形以及單個點。它有一個 AddPoint 方法，可以方便地為一組點構造邊界矩形，包括跨越180度子午線的點集。

5. Region 區域

區域表示單位球體上的二維區域。 這個介面的目的是讓複雜的區域近似為更簡單的區域。該介面僅限於計算近似值的方法。

S2 區域表示單位球體上的二維區域。它是一個具有各種具體子型別的抽象介面，如盤形，矩形，多段線，多邊形，幾何集合，緩衝形狀等。 這個介面的主要目的是使複雜區域近似為更簡單的區域。因此，介面只能用於計算近似值的方法，而不是具有各種各樣的由所有子型別實現的虛擬方法。

6. Shape 形狀

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）。

除去上面說的這幾種常用的型別以外，還有以下這些中間型別或者更加底層的型別可以供開發者使用。

• S2LatLngRect - 經緯度座標系中的矩形。
• S2Polyline - 折線。
• S2CellUnion - 一個近似為S2CellIds集合的區域。RegionCoverer 轉換以後都是這種型別。
• S2ShapeIndexRegion - 點，多義線和多邊形的任意集合。
• S2ShapeIndexBufferedRegion - 定義和 S2ShapeIndexRegion 一樣，只是擴大了給定的半徑。
• S2RegionUnion - 任意區域的集合。
• S2RegionIntersection - 任意其他區域的交集部分。

最後，額外說一句，S2RegionTermIndexer 這個型別是支援索引和查詢任何型別的S2Region，也就是上述說的所有的型別。可以使用 S2RegionTermIndexer 來索引一組多段線，然後查詢哪些多段線與給定的多邊形相交。

二. RegionCoverer 舉例

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 個。現在中間的孔就比較完整的顯示出來了。

在舉幾個我們實際專案中用到的例子。下面是上海的一個網格的邊緣。我們先用

defaultCoverer := &s2.RegionCoverer{MaxLevel: 16, MaxCells: 100, MinLevel: 13}

複製程式碼

去轉換，得到的結果如下:

defaultCoverer := &s2.RegionCoverer{MaxLevel: 30, MaxCells: 1000, MinLevel: 1}

複製程式碼

精確度提高到1000，結果就會如下：

還有一些區域更大的情況，比如一個省，湖北省：

或者一個湖，太湖：

最後在舉一個 polygon 的例子。我們知道 polygon 是由多個 loop 組成的：

	loops := []*s2.Loop{}
	loops = append(loops, loop1)
	loops = append(loops, loop2)

	polygon := s2.PolygonFromLoops(loops)

	defaultCoverer := &s2.RegionCoverer{MaxLevel: 16, MaxCells: 100, MinLevel: 8}
	fmt.Println("----  polygon-----")

	cvr = defaultCoverer.Covering(polygon)

	for i := 0; i < len(cvr); i++ {
		fmt.Printf("%d,\n", cvr[i])
	}
	fmt.Printf("------------\n")


複製程式碼

下面一次是兩個 loop 在地圖上的樣子。

最後是 polygon,它包含了以上2個 loop。

三. RegionCoverer 核心演算法 Covering 的實現

這一章節詳細分析一下 Covering 是如何實現的。

最常見的用法就是下面幾行：

	rc := &s2.RegionCoverer{MaxLevel: 30, MaxCells: 5}
	r := s2.Region(CapFromCenterArea(center, area))
	covering := rc.Covering(r)


複製程式碼

上面例子展示的是最多轉換以後 CellUnion 裡面有5個 CellID。上面這三行覆蓋的區域是一個 Cap。

type RegionCoverer struct {
	MinLevel int // the minimum cell level to be used.
	MaxLevel int // the maximum cell level to be used.
	LevelMod int // the LevelMod to be used.
	MaxCells int // the maximum desired number of cells in the approximation.
}


複製程式碼

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，這樣就又不符合要求了。那如何能在 <= MaxCells 的情況下還能最精確的覆蓋給定的區域呢？

有幾點需要提前說明的是：

• 1. MinLevel 優先順序高於 MaxCells（注意這裡不是 MaxLevel），即低於給定 Level 的 Cell 永遠都不會被使用，即使用一個它能替代掉很多面積小（Level 大）的 Cell。
• 2. 對於 MaxCells 的最小取值範圍，如果某一種情況要求的是所需的最小單元數量（例如，如果該區域與所有六個面單元相交），則可以返回多達6個單元。 如果碰巧位於三個立方體面的交點處，即使對於非常小的凸起區域也可能返回多達3個單元格。
• 3. 如果 MinLevel 對於近似的區域來說都太大了，那麼 MaxCells 是會失去約束的限制，可以返回任意數量的單元格。
• 4. 如果 MaxCells 小於4，即使該區域是凸的，比如 cap 或者 rect ，最終覆蓋的面積也要比原生區域大。所以這種情況開發者心裡要清楚。

好了，接下來從原始碼開始看起。RegionCoverer 轉換的核心函式就是這個了。

func (rc *RegionCoverer) Covering(region Region) CellUnion {
	covering := rc.CellUnion(region)
	covering.Denormalize(maxInt(0, minInt(maxLevel, rc.MinLevel)), maxInt(1, minInt(3, rc.LevelMod)))
	return covering
}

複製程式碼

從這個函式實現我們可以看到，轉換實際上就分為2步，一步是 Normalize Cell + 轉換，另外一步是 Denormalize Cell。

(一). CellUnion

CellUnion 方法的具體實現：

func (rc *RegionCoverer) CellUnion(region Region) CellUnion {
	c := rc.newCoverer()
	c.coveringInternal(region)
	cu := c.result
	cu.Normalize()
	return cu
}


複製程式碼

這個方法主要也可以分解成三個部分，新建 newCoverer、coveringInternal、Normalize。

1. newCoverer()

func (rc *RegionCoverer) newCoverer() *coverer {
	return &coverer{
		minLevel: maxInt(0, minInt(maxLevel, rc.MinLevel)),
		maxLevel: maxInt(0, minInt(maxLevel, rc.MaxLevel)),
		levelMod: maxInt(1, minInt(3, rc.LevelMod)),
		maxCells: rc.MaxCells,
	}
}


複製程式碼

newCoverer() 方法是初始化一個 coverer 的結構體。maxLevel 是一個之前定義過的常量，maxLevel = 30。coverer 的初始化引數全部都來自於 RegionCoverer 的引數。我們在外部初始化了一個 RegionCoverer ，它主要包含的4個引數，MinLevel，MaxLevel，LevelMod，MaxCells，都會傳到這裡。上面這段初始化函式裡面用到的 maxInt、minInt 主要用來進行非法值的處理。

其實 coverer 結構體裡面包含了8個元素項。

type coverer struct {
	minLevel         int // the minimum cell level to be used.
	maxLevel         int // the maximum cell level to be used.
	levelMod         int // the LevelMod to be used.
	maxCells         int // the maximum desired number of cells in the approximation.
	region           Region
	result           CellUnion
	pq               priorityQueue
	interiorCovering bool
}


複製程式碼

除去上面初始化的4項，其實它還包含其他重要的4項，這4項會在下面用到。region 要覆蓋的區域。result 就是最終轉換的結果，結果是一個 CellUnion 的陣列，pq 是優先佇列 priorityQueue，interiorCovering 是一個 bool 變數，標誌的當前轉換是否是內部轉換。

2. coveringInternal()

接下來看看 coveringInternal 方法。

func (c *coverer) coveringInternal(region Region) {
	c.region = region

	c.initialCandidates()
	for c.pq.Len() > 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 <= c.maxCells {
			for _, child := range cand.children {
				if !c.interiorCovering || len(c.result) < c.maxCells {
					c.addCandidate(child)
				}
			}
		} else {
			cand.terminal = true
			c.addCandidate(cand)
		}
	}
	c.pq.Reset()
	c.region = nil
}


複製程式碼

coveringInternal 方法會生成覆蓋的方案，並把結果儲存在 result 內。覆蓋轉換的大體策略是：

從立方體的6個面開始。丟棄任何與該區域不相交的形狀。然後重複選擇與形狀相交的最大單元格並將其細分。

coverer 結構體裡面的8個元素，前4個是外部傳進來初始化的，後4個元素就是在這裡被用到的。首先

c.region = region

複製程式碼

初始化 coverer 的區域。另外三個元素是 result、pq、interiorCovering 在下面都會被用到。

result 中只是包含將成為最終輸出的一部分的滿足條件的 Cell，而 pq 優先佇列裡面包含可能仍然需要繼續再細分的 Cell。

如果一個 Cell 100% 完全被包含在覆蓋區域內，就會被立即新增到輸出中，而完全不和該區域有任何相交的部分的 Cell 會立即丟棄。所以 pq 優先佇列中只會包含部分與該區域相交的 Cell。

pq 優先佇列出隊的策略是：

1. 首先根據 Cell 大小（首先開始是大的 Cell）優先考慮候選人
2. 然後根據相交的孩子的數量（最少的孩子優先順序高，先出列）
3. 最後按完全容納的孩子的數量（最少的孩子優先順序高，先出列）

經過 pq 優先佇列的篩選以後，最終留下來的 Cell 必定是優先順序最低的，即 Cell 面積是比較小的，並且和區域相交的部分較大且和完全容納孩子數量最多。也就是說和區域邊緣上最貼近的 Cell（Cell 覆蓋的區域比要覆蓋轉換的區域多餘的部分最少）是會被最終留下來的。

if c.interiorCovering || int(cand.cell.level) < c.minLevel || cand.numChildren == 1 || len(c.result)+c.pq.Len()+cand.numChildren <= c.maxCells {

}


複製程式碼

coveringInternal 函式實現中的這個判斷條件的意圖是：

對於內部覆蓋轉換，無論候選人有多少孩子，我們都會將其繼續不斷細分。如果我們在擴大所有孩子之前到達 MaxCells，我們將只使用其中的一些。對於外部覆蓋我們不能這樣做，因為結果必須覆蓋整個地區，所以所有的孩子都必須使用。

candidate.numChildren == 1

複製程式碼

在上述的情況下，我們已經有更多的 MaxCells 結果（minLevel太高）的情況下照顧情況。有一個孩子的候選人就算繼續細分在這種情況下對最終結果也沒有什麼影響。

3. initialCandidates()

接下來看看如何初始化候選人的：

func (c *coverer) initialCandidates() {
	// Optimization: start with a small (usually 4 cell) covering of the region's bounding cap.
	temp := &RegionCoverer{MaxLevel: c.maxLevel, LevelMod: 1, MaxCells: min(4, c.maxCells)}

	cells := temp.FastCovering(c.region)
	c.adjustCellLevels(&cells)
	for _, ci := range cells {
		c.addCandidate(c.newCandidate(CellFromCellID(ci)))
	}
}


複製程式碼

這個函式裡面有一個小優化，把區域第一步覆蓋轉換成4個 Cell 覆蓋區域邊緣的 Cap。initialCandidates 方法實現中比較重要的兩個方法是 FastCovering 和 adjustCellLevels。

4. FastCovering()

func (rc *RegionCoverer) FastCovering(region Region) CellUnion {
	c := rc.newCoverer()
	cu := CellUnion(region.CellUnionBound())
	c.normalizeCovering(&cu)
	return cu
}


複製程式碼

FastCovering 函式會返回一個 CellUnion 集合，這個集合裡面的 Cell 覆蓋了給定的區域，但是這個方法的不同之處在於，這個方法速度很快，得到的結果也比較粗糙。當然得到的 CellUnion 集合也是滿足 MaxCells, MinLevel, MaxLevel, 和 LevelMod 要求的。只不過結果不嘗試去使用 MaxCells 的大值。一般會返回少量的 Cell，所以結果比較粗糙。

所以把 FastCovering 這個函式作為遞迴細分 Cell 的起點，非常管用。

在這個方法中，會呼叫 region.CellUnionBound() 方法。這個要看各個 region 區域是如何實現這個 interface 的。

這裡以 loop 為例，loop 對 CellUnionBound() 的實現如下：

func (l *Loop) CellUnionBound() []CellID {
	return l.CapBound().CellUnionBound()
}



複製程式碼

上面方法就是快速計算邊界轉換的具體實現。也是實現空間覆蓋的核心部分。

CellUnionBound 返回覆蓋區域的 CellID 陣列。 Cell 沒有排序，可能有冗餘（例如包含其他單元格的單元格），可能覆蓋的區域多於必要的區域。

也由於上面這個理由，導致此方法不適用於客戶端程式碼直接使用。 客戶通常應該使用 Region.Covering 方法，Covering 方法可以用來控制覆蓋物的大小和準確性。 另外，如果你想快速覆蓋，不關心準確性，可以考慮呼叫 FastCovering（它會返回一個由此方法計算出來的被覆蓋的清理版本）。

CellUnionBound 實現應該嘗試返回覆蓋區域的小覆蓋（理想情況下為4個或更少），並且可以快速計算。所以 CellUnionBound 方法被 RegionCoverer 用作進一步改進的起點。

func (l *Loop) CapBound() Cap {
	return l.bound.CapBound()
}

複製程式碼

CapBound 返回一個邊界上限，它可能比相應的 RectBound 會有更多的填充。它的邊界是保守的，如果 loop 包含點P，那麼邊界也一定包含這個點。

func (c Cap) CellUnionBound() []CellID {

	level := MinWidthMetric.MaxLevel(c.Radius().Radians()) - 1

	// If level < 0, more than three face cells are required.
	if level < 0 {
		cellIDs := make([]CellID, 6)
		for face := 0; face < 6; face++ {
			cellIDs[face] = CellIDFromFace(face)
		}
		return cellIDs
	}

	return cellIDFromPoint(c.center).VertexNeighbors(level)
}


複製程式碼

上面這段程式碼就是轉換的最粗糙版本的核心演算法了。在這段演算法裡面核心一步就是算出要找的 level。

level := MinWidthMetric.MaxLevel(c.Radius().Radians()) - 1

複製程式碼

這裡找到的 Level 就是該 Cap 所能包含的最大的 Cell。

return cellIDFromPoint(c.center).VertexNeighbors(level)

複製程式碼

上面這句就是返回了 4 個 Cell，距離 Cap 中心點最近的。當然，如果 Cap 非常大，有可能會返回 6 個 Cell。當然，返回的這些 Cell 是沒有經過任何排序的。

5. normalizeCovering()

FastCovering 的最後一步就是 normalizeCovering。

func (c *coverer) normalizeCovering(covering *CellUnion) {
	// 1
	// 
	if c.maxLevel < maxLevel || c.levelMod > 1 {
		for i, ci := range *covering {
			level := ci.Level()
			newLevel := c.adjustLevel(minInt(level, c.maxLevel))
			if newLevel != level {
				(*covering)[i] = ci.Parent(newLevel)
			}
		}
	}
	// 2
	// 
	covering.Normalize()

	// 3
	// 
	for len(*covering) > c.maxCells {
		bestIndex := -1
		bestLevel := -1
		for i := 0; i+1 < len(*covering); i++ {
			level, ok := (*covering)[i].CommonAncestorLevel((*covering)[i+1])
			if !ok {
				continue
			}
			level = c.adjustLevel(level)
			if level > bestLevel {
				bestLevel = level
				bestIndex = i
			}
		}

		if bestLevel < c.minLevel {
			break
		}
		(*covering)[bestIndex] = (*covering)[bestIndex].Parent(bestLevel)
		covering.Normalize()
	}
	// 4
	// 
	if c.minLevel > 0 || c.levelMod > 1 {
		covering.Denormalize(c.minLevel, c.levelMod)
	}
}



複製程式碼

normalizeCovering 會對前一步的覆蓋轉換結果進行進一步的規範化，使其符合當前覆蓋引數（MaxCells，minLevel，maxLevel和levelMod）。 這種方法不會嘗試最佳結果。 特別是，如果minLevel> 0或者levelMod> 1，那麼即使這不是必需的，它也可能返回比期望的 Cell 更多的值。

在上面的程式碼實現中標註了4處需要注意的地方。

第一處，判斷的是，如果 Cell 太小了，或者不滿足 levelMod 的條件，就用他們的祖先來替換掉他們。

第二處，是對前一步的結果排序並且進一步簡化。

第三處，如果 Cell 數量上還是太多，仍然有太多的 Cell，則用 for 迴圈去查詢兩個相鄰的 Cell 的最近公共祖先 LCA 並替換掉它們倆，for 迴圈的順序就是以 CellID 的順序進行迴圈的。

這裡用到的 LCA 的具體實現在前面一篇文章裡面有詳解，文章連結，這裡就不再贅述了。

第四處，最後確保得到的結果能滿足 minLevel 和 levelMod，最好也能滿足 MaxCells。

接下來還需要進一步分析的就是 Normalize() 和 Denormalize() 這兩個函式實現了。

6. Normalize()

先看 Normalize()。

func (cu *CellUnion) Normalize() {
	sortCellIDs(*cu)

	output := make([]CellID, 0, len(*cu)) // the list of accepted cells
	
	for _, ci := range *cu {
	
		// 1
		if len(output) > 0 && output[len(output)-1].Contains(ci) {
			continue
		}
		
		// 2
		j := len(output) - 1 // last index to keep
		for j >= 0 {
			if !ci.Contains(output[j]) {
				break
			}
			j--
		}
		output = output[:j+1]

		// 3
		for len(output) >= 3 && areSiblings(output[len(output)-3], output[len(output)-2], output[len(output)-1], ci) {
			output = output[:len(output)-3]
			ci = ci.immediateParent() // checked !ci.isFace above
		}
		output = append(output, ci)
	}
	*cu = output
}


複製程式碼

Normalize 這個方法的主要意圖是想整理 CellUnion 中各個 CellID，經過排序輸出沒有冗餘的 CellUnion。

排序是第一步。

接下來整理冗餘的 Cell 是第二步，也是這個函式實現裡面關鍵的一步。冗餘分為2種，一種是完全包含，另外一種是4個小的 Cell 可以合併成一個大的。

先處理是否有一個 Cell 完全包含另外一個 Cell 的情況，這種情況下，被包含的那個 Cell 就是冗餘的 Cell，就該被丟棄。對應的是上述程式碼中標1的地方。

在實現上，我們只需要檢查最後接受的單元格。 Cell 首先經過排序以後，如果當前這個候選的 Cell 不被最後接受的 Cell 所包含，那麼它就不能被任何先前接受的 Cell 所包含。

同理，如果當前候選的 Cell 包含了之前已經接受過檢查的 Cell，那麼之前已經在 output 裡面的 Cell 也需要被丟棄掉。因為 output 維護的是一個連續的尾序列，前面也提到了 S2 Cell 是被排序了，所以這裡就不能破壞它的連續性。這裡對應的是上述程式碼中標2的地方。

最後程式碼中標3的地方是看是否最後三個單元格加上這個可以合併。 如果連續的3個 Cell 加上當前的 Cell 可以級聯到最近的一個父節點上，我們就把它們三個用大的 Cell 替換掉。

經過 Normalize 這一步的“格式化”以後，輸出的 Cell 都是有序且無冗餘的 Cell。

7. Denormalize()

func (cu *CellUnion) Denormalize(minLevel, levelMod int) {
	var denorm CellUnion
	for _, id := range *cu {
		level := id.Level()
		newLevel := level
		if newLevel < minLevel {
			newLevel = minLevel
		}
		if levelMod > 1 {
			newLevel += (maxLevel - (newLevel - minLevel)) % levelMod
			if newLevel > maxLevel {
				newLevel = maxLevel
			}
		}
		if newLevel == level {
			denorm = append(denorm, id)
		} else {
			end := id.ChildEndAtLevel(newLevel)
			for ci := id.ChildBeginAtLevel(newLevel); ci != end; ci = ci.Next() {
				denorm = append(denorm, ci)
			}
		}
	}
	*cu = denorm
}


複製程式碼

Denormalize 這個函式實現很簡單了，它是 normalize 的另一面（名稱上雖然是反義）。這個函式是為了“規範” Cell 是否滿足覆蓋轉換之前的幾個預定條件：MinLevel、MaxLevel、LevelMOD、MaxCell。

所有 level 小於 minLevel 或者（level-minLevel）不是 levelMod 的倍數的任何 Cell 都會被它的孩子節點替換，直到滿足這兩個條件或者達到 maxLevel。

Denormalize 函式的意圖是確保得到的結果能滿足 minLevel 和 levelMod，最好也能滿足 MaxCells。

到這裡讀者也應該清楚為何函式名叫 Denormalize 了吧，為了滿足條件，它把大的 Cell 替換成了自己的孩子，小的 Cell，而 Normalize 正好相反，是把4個小的孩子 Cell 用它們直系父親節點去替換。

分析到這裡 FastCovering() 函式也就都分析完畢了。

8. adjustCellLevels()

在回到 initialCandidates() 函式中，在這個函式中 FastCovering() 之後還有一步操作，adjustCellLevels。

c.adjustCellLevels(&cells)

複製程式碼

接下來就看看 adjustCellLevels 的具體實現。

func (c *coverer) adjustCellLevels(cells *CellUnion) {
	if c.levelMod == 1 {
		return
	}

	var out int
	for _, ci := range *cells {
		level := ci.Level()
		newLevel := c.adjustLevel(level)
		if newLevel != level {
			ci = ci.Parent(newLevel)
		}
		if out > 0 && (*cells)[out-1].Contains(ci) {
			continue
		}
		for out > 0 && ci.Contains((*cells)[out-1]) {
			out--
		}
		(*cells)[out] = ci
		out++
	}
	*cells = (*cells)[:out]
}


複製程式碼

adjustCellLevels 是用來確保 level> minLevel 的所有 Cell 也能滿足 levelMod，必要時可以用祖先替換它們。 小於 minLevel 的這些 Cell ，它們的 Level 的不會被修改（請參閱調整級別）。 最終得到的結果是標準化以確保不存在冗餘單元的 CellUnion。

adjustCellLevels 和 Denormalize 有點類似，都是為了滿足條件調整 CellUnion。但是兩者調整的方向不同，Denormalize 是把 Cell 往孩子的方向替換，adjustCellLevels 是把 Cell 往父親節點的方向替換。

func (c *coverer) adjustLevel(level int) int {
	if c.levelMod > 1 && level > c.minLevel {
		level -= (level - c.minLevel) % c.levelMod
	}
	return level
}


複製程式碼

adjustLevel 意圖是為了返回更小一級的 Level，以使其滿足 levelMod 條件。 小於minLevel的 level 不受影響（因為這些 level 的單元最終會被 Denormalize 函式處理）。

（二）. Denormalize

關於 Denormalize 的實現已經在上面分析過了。這裡就不再分析了。

（三）. 小結

用一張圖來表示 Google S2 對整個空間區域覆蓋演算法實現的全流程：

上圖中每個關鍵實現都分析過了，哪個節點還不明白的同學可以回過頭往上再翻一翻。

這個近似演算法並不是最優演算法，但是在實踐中效果還不錯。輸出的結果並不總是使用的滿足條件的最多的 Cell 個數，因為這樣也不是總能產生更好的近似結果(比如上面舉例的，待覆蓋的區域整好位於三個面的交點處，得到的結果比原區域要大很多) 並且 MaxCells 對搜尋的工作量和最終輸出的 Cell 的數量是一種限制。

由於這是一個近似演算法，所以不能依賴它輸出的穩定性。特別的，覆蓋演算法的輸出結果會在不同的庫的版本上有所不同。

這個演算法還可以產生內部覆蓋的 Cell，內部覆蓋的 Cell 指的是完全被包含在區域內的 Cell。如果沒有滿足條件的 Cell ，即使對於非空區域，內部覆蓋 Cell 也可能是空的。

請注意，處於效能考慮，在計算內部覆蓋 Cell 的時候，指定 MaxLevel 是明智的做法。否則，對於小的或者零面積的區域，演算法可能會花費大量時間將 Cell 細分到葉子 level ，以嘗試找到滿足條件的內部覆蓋的 Cell。

四. 最後

關於空間搜尋系列文章到這裡也告一段落了，最後當然說一些心得體會了。

在學習和實踐空間搜尋這塊知識的時候，筆者我檢視了物理，數學，演算法這三方面的資料，從物理和數學2個層次提升了整個空間和時間的認知。雖然目前個人在這方面的認知也許還很淺顯，不過對比之前實在是進步了很多。目的也達到了。

最後的最後就是推薦2個網站，這也是微博上提問問到最多的。

第一個問題是：系列文章裡面的 S2 Cell 都是怎麼畫出來的？

這其實是一個人的開源網站，s2map.com/，筆者是在這裡填入程式算好的 CellID，然後顯示出來的。這就相當於是 S2 的視覺化研究展示工具了。

第二個問題是：為什麼有些程式碼在 Go 的版本里面沒有找到相關實現？

答案是 Go 的版本實現完成度還沒有到 100%，個別的還需要去參考 C++ 和 Java 完整版的程式碼。關於 C++ 和 Java 的原始碼，Google 已經在幾天前把程式碼從私有程式碼倉庫移到了 Github 上了。更加方便學習與檢視了。官方也把一些文件整理到了 s2geometry.io/ 這個網站上。初學者建議還是先看官方 API 說明文件。在看完文件以後，原理性的問題還有些疑惑，可以來翻翻筆者這個空間搜尋系列文章，希望能對讀者有幫助。

---

空間搜尋系列文章：

如何理解 n 維空間和 n 維時空
高效的多維空間點索引演算法 — Geohash 和 Google S2
Google S2 中的 CellID 是如何生成的 ？
Google S2 中的四叉樹求 LCA 最近公共祖先
神奇的德布魯因序列
四叉樹上如何求希爾伯特曲線的鄰居 ？
Google S2 是如何解決空間覆蓋最優解問題的?

GitHub Repo：Halfrost-Field

Follow: halfrost · GitHub

Source: halfrost.com/go_s2_regio…