【算法】推步聚顶计算多边形的面积

计算多边形的面积

“数性至朴，算学是天下最诚实的东西，一加一永远是二，五乘四永远是二十，而十二自实永远是一百四十四。”

这世上，人心比算学更复杂。人际关系让人疲惫，还是数学比较单纯。这是我在电视剧《显微镜下的大明：丝绢案》中看到的，“人心叵测，难以琢磨，而数性至朴，从不虚饰，数与数之间的关联，就好比天上的星宿一样，万年难易，只要你掌握了其中的运转之妙，计算之理，便可以俯仰天地，上下求索”，帅家默说，“这便是我所求的道，任凭数字无穷无尽，任凭难题变化莫测，吾以一道而御之。也许在你看来，这只是一间装满了断烂朝报的破屋子，可在我看来，这是直通仙府的大道，我在其中神游天外，上穷碧落下黄泉，无数的景致只在一念之间，何其的快活。”

我留意到剧中有一个丈量妖田的算法————推步聚顶之术，于是心血来潮研究了一下，总结几种计算多边形面积的方法。

算法一：三角剖分法

这个算法使用了多边形三角剖分的思想来计算多边形的面积。具体来说，它将多边形划分为若干个三角形，然后计算每个三角形的面积，最后将所有三角形的面积相加得到多边形的面积。

具体实现步骤如下：

首先，将多边形的顶点按照顺序排列。
从任意一个顶点开始，顺次连接所有相邻的顶点，形成若干个三角形。
对于每个三角形，使用海龙公式（海伦公式）计算其面积。海龙公式的计算公式为：

S = √[p(p-a)(p-b)(p-c)]

其中，S 表示三角形的面积，a、b、c 分别为三角形的三条边的长度，p 表示半周长，即 (a+b+c)/2。
将所有三角形的面积相加，即可得到多边形的面积。

在实现时，我们可以用一个循环遍历多边形的每个顶点，计算每个三角形的面积，并将其加入总面积中。具体地，我们可以计算以当前顶点和相邻的两个顶点构成的三角形面积，再累加到总面积中。

需要注意的是，这个算法要求多边形的顶点按照顺序排列。如果顺序不正确，计算结果将不正确。

// 三角剖分法的时间复杂度通常为 O(n log n) 或 O(n^2)
func polygonArea(x, y []float64) float64 {n := len(x)area := 0.0for i := 0; i < n; i++ {j := (i + 1) % narea += x[i] * y[j] - x[j] * y[i]}area = math.Abs(area) / 2.0// Triangulate the polygontriArea := 0.0for i := 1; i < n-1; i++ {triArea += triangleArea(x[0], y[0], x[i], y[i], x[i+1], y[i+1])}return triArea
}func triangleArea(x1, y1, x2, y2, x3, y3 float64) float64 {a := math.Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2))b := math.Sqrt((x2-x3)*(x2-x3) + (y2-y3)*(y2-y3))c := math.Sqrt((x3-x1)*(x3-x1) + (y3-y1)*(y3-y1))s := (a + b + c) / 2.0return math.Sqrt(s * (s - a) * (s - b) * (s - c))
}

计算三角形面积可以使用叉积来实现，即 S = 0.5 * (x1y2 - x2y1 + x2y3 - x3y2 + … + xny1 - x1yn)，其中 x 和 y 分别表示多边形顶点的 x 坐标和 y 坐标。

叉积是向量运算中的一种，用于计算两个向量构成的平行四边形的面积、方向和法线。叉积也称为叉积积、矢量积或向量积。

两个向量 a 和 b 的叉积记为 a × b，其结果是一个新的向量，其大小等于 a 和 b 所在平行四边形的面积，方向垂直于 a 和 b 所在平面，且满足右手定则。右手定则是一个基本规律，它规定：将右手伸开，将大拇指指向向量 a 的方向，食指指向向量 b 的方向，那么中指所指的方向就是 a × b 的方向。

计算 a 和 b 叉积的公式为：

a × b = ∣ a ∣ ∣ b ∣ s i n θ n a × b =|a||b|sinθn a×b=∣a∣∣b∣sinθn

其中，|a| 和 |b| 分别表示向量 a 和 b 的长度，θ 表示 a 和 b 之间的夹角，n 表示垂直于 a 和 b 所在平面的单位向量。

在计算多边形面积时，我们可以使用叉积来计算每个三角形的面积。具体地，假设三角形的两个边向量分别为 a 和 b，则三角形的面积为 0.5 * |a × b|。这是因为 a 和 b 所在平行四边形的面积就是 a × b 的大小，而三角形的面积等于平行四边形的面积的一半。

叉积也常用于计算向量之间的夹角、判断向量之间的关系、计算平面方程等。

算法二：高斯公式法（又称原点内切圆法）

高斯公式法是一种基于环绕定理的方法。它将多边形分解为若干个三角形，然后分别计算每个三角形的面积，并将其相加得到多边形的面积。具体来说，假设多边形的顶点按照顺序排列，则其面积为：

S = 0.5 * |Σ(xiyi+1 - xi+1yi)|

其中，xi 和 yi 分别表示第 i 个顶点的 x 坐标和 y 坐标，n 表示多边形的顶点数。

// 高斯公式计算多边形面积的时间复杂度为 O(n)，其中 n 表示多边形的顶点数。
// 具体地，这个算法的时间复杂度主要来自于计算多边形顶点坐标的行列式，其计算复杂度为 O(n^3)。但由于这个行列式是一个三角矩阵，因此实际的计算复杂度可以通过高斯消元算法等方法降低到 O(n^2)。
func polygonArea(x, y []float64) float64 {n := len(x)area := 0.0for i := 0; i < n; i++ {j := (i + 1) % narea += x[i] * y[j] - x[j] * y[i]}return 0.5 * area
}

通过计算多边形顶点坐标的行列式来计算多边形面积。它要求多边形顶点必须按照顺序排列，并且不能有内部空洞。实现起来也比较简单，时间复杂度为 O(n)，但是只适用于简单多边形。

另外，三角剖分法和高斯公式法的计算精度也不同。三角剖分法的精度一般较高，但在计算复杂多边形时可能会有一些误差。高斯公式法的精度也较高，但在顶点数较多时可能会导致计算溢出或精度损失。

因此，选择哪种方法应该根据实际需要和情况来确定。如果多边形较为简单或需要高精度计算，可以选择高斯公式法。如果多边形较为复杂或需要计算带有内部空洞的多边形，可以选择三角剖分法。

算法三：扫描线法

将多边形投影到 x 轴上，然后从下往上扫描每个像素，计算其所属的多边形。具体地，该算法通过扫描线的方式逐行扫描多边形的每个像素，并检查该像素是否在多边形内部。为了检查像素是否在多边形内部，我们需要先对多边形进行排序，然后根据多边形的每一条边来判断像素是否在多边形内部。

在实际的应用中，扫描线法通常用于计算二维简单多边形的面积、边界和重心等。该算法的时间复杂度为 O(n log n)，其中 n 表示多边形的顶点数。需要注意的是，对于带有内部空洞的多边形，扫描线法需要进行额外的处理，否则可能无法正确计算多边形的面积。

//  时间复杂度 O(nlogn)
type point struct {x float64y float64
}func polygonArea(vertices []point) float64 {n := len(vertices)if n < 3 {return 0}// Sort vertices by y-coordinatefor i := 0; i < n-1; i++ {for j := i + 1; j < n; j++ {if vertices[i].y > vertices[j].y {vertices[i], vertices[j] = vertices[j], vertices[i]}}}// Compute area using scanline algorithmarea := 0.0for y := vertices[0].y + 1; y <= vertices[n-1].y; y++ {var intersections []float64for i := 0; i < n; i++ {j := (i + 1) % nif (vertices[i].y <= y && y < vertices[j].y) || (vertices[j].y <= y && y < vertices[i].y) {x := (y-vertices[i].y)*(vertices[j].x-vertices[i].x)/(vertices[j].y-vertices[i].y) + vertices[i].xintersections = append(intersections, x)}}intersections = append(intersections, math.Inf(1))// Sort intersections by x-coordinatefor i := 0; i < len(intersections)-1; i++ {for j := i + 1; j < len(intersections); j++ {if intersections[i] > intersections[j] {intersections[i], intersections[j] = intersections[j], intersections[i]}}}// Compute area contribution of current scanlinefor i := 0; i < len(intersections)-1; i++ {x1 := intersections[i]x2 := intersections[i+1]area += 0.5 * (y - vertices[0].y) * (x2 - x1)}}return area
}

使用了一个包含多边形顶点的点结构体来表示多边形，即 point。polygonArea 函数将多边形顶点按照 y 坐标从小到大进行排序，然后使用扫描线算法计算多边形的面积。

算法四：分治法

分治法将多边形递归地分成若干个简单多边形，然后计算每个简单多边形的面积，最后将其相加得到多边形的面积。这个算法时间复杂度为 O(n log n)。

// 时间复杂度为 O(n log n),每次递归将多边形分割成两个子多边形，因此递归树的高度为 log n。对于每个子多边形，计算其面积的时间复杂度为 O(1)，因此总时间复杂度为 O(n log n)。
_type Point struct {X float64Y float64
}func PolygonArea(vertices []Point) float64 {if len(vertices) < 3 {return 0.0}return dividePolygonArea(vertices, 0, len(vertices)-1)
}func dividePolygonArea(vertices []Point, start int, end int) float64 {if end-start == 2 {// 三角形return triangleArea(vertices[start], vertices[start+1], vertices[end])} else if end-start == 3 {// 梯形和三角形return trapezoidArea(vertices[start], vertices[start+1], vertices[start+2], vertices[end])}// 将多边形分割成两部分mid := (start + end) / 2leftArea := dividePolygonArea(vertices, start, mid)rightArea := dividePolygonArea(vertices, mid, end)return leftArea + rightArea
}func triangleArea(p1 Point, p2 Point, p3 Point) float64 {return 0.5 * (p1.X*p2.Y + p2.X*p3.Y + p3.X*p1.Y - p1.Y*p2.X - p2.Y*p3.X - p3.Y*p1.X)
}func trapezoidArea(p1 Point, p2 Point, p3 Point, p4 Point) float64 {return 0.5 * (p2.Y-p1.Y + p3.Y-p4.Y) * (p4.X-p1.X)
}func main() {vertices := []Point{{0, 0}, {1, 0}, {2, 1}, {1, 2}}area := PolygonArea(vertices)fmt.Printf("The area of the polygon is %f\n", area)
}

算法五：推步聚顶

先牵经纬以衡量，再点原初标步长。田型取顶分别数，再算推步知地方。

先牵经纬以衡量：先按照经纬度顺序排列顶点，以便进行下一步计算；
再点原初标步长：然后确定原初标准的步长；
田型取顶分别数：接下来，取出所有顶点的最高点和最低点，并进行编号；
再算推步知地方：然后根据步长推算出每个顶点所在的地方，并将其分别编号；
按照编号的顺序进行推步聚顶计算，得到多边形的三角剖分。

该算法的基本思想是将多边形的顶点按照某种规则重新排列，使得在新的顶点顺序下，多边形的三角剖分更加容易进行。具体步骤如下：

初始化顶点序列 V 和边序列 E，将顶点按顺序插入 V 中，将边插入 E 中。
对 V 中的顶点进行重新排列，使得相邻两个顶点之间的边越短越好。
对顶点序列 V 中的每个顶点，判断其与相邻顶点连成的边是否为多边形的一条对角线，如果是，则将其剪切掉，并将多边形分成两个子多边形。
对每个子多边形递归进行三角剖分，最终将所有子多边形的三角形合并为整个多边形的三角形剖分。

type Point struct {X float64Y float64
}
func VertexReordering(vertices []point) []point {// 按照经纬度排序vertices = sortByLongitude(vertices)// 确定原初标准的步长step := computeStep(vertices)// 取出所有顶点的最高点和最低点，并进行编号top, bottom := getTopAndBottompoints(vertices)_, bottomIndex := getTopAndBottomIndices(vertices, top, bottom)// 推算出每个顶点所在的地方，并将其分别编号columnCount := int(math.Ceil((top.y - bottom.y) / step))vertexColumns := make([][]int, columnCount)for i := 0; i < len(vertices); i++ {column := int(math.Floor((vertices[i].y - bottom.y) / step))vertexColumns[column] = append(vertexColumns[column], i)}// 按照编号的顺序进行推步聚顶计算triangles := make([][3]int, 0)for i := 0; i < columnCount-1; i++ {column1 := vertexColumns[i]column2 := vertexColumns[i+1]for j := 0; j < len(column1); j++ {index1 := column1[j]for k := 0; k < len(column2); k++ {index2 := column2[k]if isDiagonal(vertices, index1, index2) {triangles = append(triangles, [3]int{index1, index2, bottomIndex + i})triangles = append(triangles, [3]int{index2, index2, bottomIndex + i})triangles = append(triangles, [3]int{index2, index1, bottomIndex + i + 1})break}}}}return vertices
}func sortByLongitude(vertices []point) []point {// TODO: 根据经度排序return vertices
}func computeStep(vertices []point) float64 {// TODO: 确定原初标准的步长return 0.1
}func getTopAndBottompoints(vertices []point) (point, point) {// TODO: 取出最高点和最低点return point{}, point{}
}func getTopAndBottomIndices(vertices []point, top point, bottom point) (int, int) {// TODO: 取出最高点和最低点的编号return 0, 0
}func isDiagonal(vertices []point, i int, j int) bool {// TODO: 判断是否为多边形的对角线return false
}
tices []point) (point, point) {// TODO: 取出最高点和最低点return point{}, point{}
}func getTopAndBottomIndices(vertices []point, top point, bottom point) (int, int) {// TODO: 取出最高点和最低点的编号return 0, 0
}func isDiagonal(vertices []point, i int, j int) bool {// TODO: 判断是否为多边形的对角线return false
}

代码地址：计算多边形面积