計算一個多邊形的重心點座標 (Swift 程式碼實現)
在之前的 《如何判斷一個多邊形是否合法》 一文中有提到，用無人機規劃飛行路線前，往往需要框選一個多邊形的區域。
而在地圖控制元件上顯示這個多邊形區域時，往往會遇到這樣一個需求：需要把所要測繪的多邊形區域移動到地圖中心。
實現這個需求的基本思路就是：獲取到多邊形區域的重心點座標，然後利用地圖控制元件的 setCenter 方法，就可以把地圖的顯示中心移動到多邊形區域重心了。那麼問題來了，如何求出一個多邊形的重心點座標呢？
這裡所說的重心，也常常叫幾何中心
這裡首先給出一個公式:
平面多邊形 可以被剖分為 n個有限的簡單圖形 ，這些簡單圖形的重心點為 ，面積為 ，那麼這個平面多邊形的重心點座標為
公式參考:維基百科
一般來說我們可以給多邊形進行三角剖分，而 即為多邊形的總面積，那麼這個公式可以理解為：
多邊形重心橫座標 = 多邊形剖分的每一個三角形重心的橫座標 * 該三角形的面積之和 / 多邊形總面積
多邊形重心縱座標 = 多邊形剖分的每一個三角形重心的縱座標 * 該三角形的面積之和 / 多邊形總面積
所以這裡就把問題拆分成了三個小問題：
- 求每個剖分出來的三角形的重心。
- 求每個剖分出來的三角形的面積。
- 求多邊形的面積。
演算法解析
1. 求三角形的重心
三角形的重心：三條中線的交點。其中重心到其中一個頂點的距離是重心到該頂點對邊中點的距離的2倍。
即：GC = 2 * GP，也就是說重心座標在 CP 線段上距離 AB 的中點 P 的 1/3 處。 假設 A,B,C 三點的座標為：
那麼通過簡單座標計算，可以得出其重心座標為
2. 求三角形面積
計算三角形的面積，我們這裡利用 向量積 來計算，我們知道平面中的兩個向量的叉乘的模等於以這兩個向量為邊的平行四邊形的面積，那麼以這個兩個向量為邊的三角形，則是這個平行四邊形的面積的一半。
參考：向量叉積
如上圖，已知平面上兩點 ，以 A，B和座標原點 構成的三角形的面積 S 為：
這裡給出運算草稿：
為什麼這裡我們會以原點作為第三個點構成三角形呢？其實是跟接下來求多邊形面積是有關聯的。
3. 求多邊形的面積
我們在上面給出的求平面多邊形重心的公式中有說到，一般我們會把多邊形剖分為多個三角形。 那麼這個剖分點 P 我們可以設在哪裡呢？這裡先給出結論：這個剖分點可以設定在多邊形的內部，也可以設定到外部。
為什麼這個剖分點可以設定到外部呢？我們可以通過簡單的三角形情況來推廣到多邊形的情況。 對於三角形ABC，我們把剖分點設定在其外部 P 的一點上
如果大家還記得 《如何判斷一個多邊形是否合法》 一文中有講過向量叉積是有正負之分的，並且根據上面所說的計算三角形面積，那麼以 P 為剖分點，通過向量積可以得出這個三角形的面積 A 為：
因為 向量PB 在 向量PA 的順時針方向，所以 的結果是負數的。那麼上面的面積計算公式其實就可以理解為：
三角形ABC的面積 = 三角形PBC面積 + 三角形PCA面積 - 三角形PAB面積
假設這四個點的座標為： ，通過上面的公式進行計算，具體的演算過程我就不給出了，這裡直接給出計算結果：
我們可以發現，計算結果中沒有 的項，因為它們在計算過程中給消去了，數學就是這麼奇妙！所以我們可以得出一個結論，多邊形的面積結果與這個剖分點的位置是無關的。那麼為了計算方便，我們當然選擇把這個 P 點設定到原點上啦。
那麼只要我們知道多邊形的每一個頂點，通過原點進行剖分成多個三角形，然後通過向量的叉乘求出每個三角的面積，最後相加，就可以求出多邊形的面積了。
示例程式碼及解析
好了，說到這裡，我們已經找到所有滿足最開始的計算多邊形重心點座標的所有計算元素了。是時候上程式碼了，這裡構建一個函式 calculatePolygonGravityCenter(coordinates: [CLLocationCoordinate2D]) ，這個函式傳入的引數是多邊形在地圖上的座標點陣列。
func calculatePolygonGravityCenter(coordinates: [CLLocationCoordinate2D]) -> CLLocationCoordinate2D {
var area = 0.0 // 多邊形面積
var gravityLat = 0.0 // 重心點 latitude
var gravityLng = 0.0 // 重心點 longitude
for (index, coordinate) in coordinates.enumerated() {
// 1
let lat = coordinate.latitude
let lng = coordinate.longitude
let nextLat = coordinates[(index + 1) % coordinates.count].latitude
let nextLng = coordinates[(index + 1) % coordinates.count].longitude
// 2
let tempArea = (nextLat * lng - nextLng * lat) / 2.0
// 3
area += tempArea
// 4
gravityLat += tempArea * (lat + nextLat) / 3
gravityLng += tempArea * (lng + nextLng) / 3
}
// 5
gravityLat = gravityLat / area
gravityLng = gravityLng / area
return CLLocationCoordinate2D(latitude: gravityLat, longitude: gravityLng)
}
複製程式碼
對應上面程式碼的註釋：
- 拿到多邊形上連續兩個點的座標，我們可以把 latitude 看做橫座標，longitude 是縱座標。
- 利用向量叉乘計算這兩個點與原點組成的三角形的面積。
- 所有面積之和得出多邊形的面積，就是求公式 中的 。
-
(lat + nextLat) / 3是以這兩個點和原點組成的三角形的重心橫座標，這樣的累加gravityLat += tempArea * (lat + nextLat) / 3其實是求公式 中的 的值。 - 到這一步就簡單了，直接套用公式 。