首先界定一下本算法中多边形的范围,不存在自交和内环。如果多边形存在自交,需要在自交位置分割再应用本算法;如果包含内环,本文算法能单独计算每个环的面积。
更多几何算法,欢迎关注公众号:几何算法
算法一:三角形面积累加法
先看一个最简单的情况,如下图所示。凸多边形,我们选择的点恰巧在多边形内部,这种情况很直观的就能知道,多边形的面积就等于这些三角形面积的累加。
三角形的面积可以通过向量叉积来计算,三角形ABC面积计算公式如下:
公式 1 : S = ∣ ( P C − P A ) × ( P B − P A ) ∣ 2 \space \space \space \space \space \space {\color{Red}公式1:} \space \space \space S=\frac{|(P_C - P_A)×(P_B - P_A)|}{2} 公式1: S=2∣(PC−PA)×(PB−PA)∣
关于向量叉积,可参考如下文章
下面探讨点在逆时针凸多边形外部的情况,如下图所示。 多边形的面积不再是所有三角形面积的和,三角形OBC、OCD、ODE、OEA面积的和再减去三角形OBA的面积。
这时候如果我们让三角形OBC、OCD、ODE、OEA的面积为正,三角形OAB的面积为负,那么就和前面简单的例子统一了起来,多边形的面积同样是它们面积的累加了。
如上图所示,假设多边形顶点按逆时针旋转(顶点顺序为A->B->C->D->E),面积为正的三角形,O在边的左侧;面积为负的三角形,O在边的右侧。
以三角形OCD为例,沿着C->D的方向看,O在CD的左侧,三角形面积为正,此时向量OC与向量OD的叉积也为正;同理,在三角形OAB中,沿A->B方向看,O在AB的右侧,这时候三角形面积为负,此时向量OA与向量OB的叉积也为负。
我们规定面积的符号和叉积的符号相同。
这时,当多边形是逆时针时,去掉公式1中的绝对值,我们就得到了带符号的三角形ABC的面积公式,如下所示:
公式 2 : S = ( P C − P A ) × ( P B − P A ) 2 \space \space \space \space \space \space {\color{Red}公式2:} \space \space \space S=\frac{(P_C - P_A)×(P_B - P_A)}{2} 公式2: S=2(PC−PA)×(PB−PA)
这些三角形带符号面积的累加就得到了逆时针多边形的面积,公式如下:
公式 3 : S ′ = ∑ i = 0 n − 1 ( P i − 1 − P O ) × ( P i − P O ) 2 \space \space \space \space \space \space {\color{Red}公式3:} \space \space \space S^{'}=\sum_{i=0}^{n-1}\frac{(P_{i-1} - P_O)×(P_{i} - P_O)}{2} 公式3: S′=i=0∑n−12(Pi−1−PO)×(Pi−PO)
如果多边形是顺时针的,和前面同理,所有三角形面积的符号正好取反,这时候累加出来的多边形面积是实际面积的相反数。
也就是说,无论多边形顺逆时针,对公式3取绝对值就能得到最终的面积,公式如下:
公式 4 : S = ∣ ∑ i = 0 n − 1 ( P i − 1 − P O ) × ( P i − P O ) 2 ∣ \space \space \space \space \space \space {\color{Red}公式4:} \space \space \space S=|\sum_{i=0}^{n-1}\frac{(P_{i-1} - P_O)×(P_{i} - P_O)}{2}| 公式4: S=∣i=0∑n−12(Pi−1−PO)×(Pi−PO)∣
反过来,我们可以用公式3来判断多边形的顺逆时针。公式3为正,多边形是逆时针;为负,多边形是顺时针。
同样可验证,对于凹多边形,公式3和公式4仍然是成立的。凹多边形和凸多边形的区别是,凸多边形外被三角形覆盖的区域加了1次又减了1次,最终是不加不减;凹多边形是加了n次又减了n次,最终同样是不加不减。也就是说对任意多边形,公式3和公式4成立。
到这里我们就得到了一个通用多边形面积计算算法(公式4)和一个多边形顺逆时针判断的算法(公式3的正负)。
如果选择的O点为(0,0)点,公式4可简化如下:
公式 5 : S = ∣ ∑ i = 0 n − 1 P i − 1 × P i 2 ∣ \space \space \space \space \space \space {\color{Red}公式5:} \space \space \space S=|\sum_{i=0}^{n-1}\frac{P_{i-1}×P_{i}}{2}| 公式5: S=∣i=0∑n−12Pi−1×Pi∣
这个公式有个专门的名字叫鞋带公式。
算法二:梯形累加
首先还是先看个简单情况,假设多边形是逆时针的凸多边形,整个多边形在x轴上方。我们可以用多边形的每条边和x轴形成一个如下图所示的梯形。可以看出,多边形的面积就是位于上面的边形成的所以梯形面积减去位于下面的边形成的所有梯形面积。
当我们让上面边形成梯形的面积为正,下面边形成梯形的面积为负时,多边形的面积就是所有梯形面积的累加。
下面重点是如何让上面边梯形面积为正,下面边梯形面积为负。
观察上图可发现,假设多边形顶点按逆时针旋转(顶点顺序为A->B->C->D->E),上面的边和下面的边的顶点y坐标都是大于0的,上面的边朝向x减少的方向,下面的边朝向x增大的方向。所以我们用下面公式计算的梯形面积恰巧是满足上面边形成梯形的面积为正,下面边形成梯形的面积为负的。
公式 6 : S = ( y i − 1 + y i ) × ( x i − 1 − x i ) 2 \space \space \space \space \space \space {\color{Red}公式6:} \space \space \space S=\frac{(y_{i-1} + y_{i})×(x_{i-1} - x_{i})}{2} 公式6: S=2(yi−1+yi)×(xi−1−xi)
多边形面积公式为
公式 7 : S = ∑ i = 0 n − 1 ( y i − 1 + y i ) × ( x i − 1 − x i ) 2 \space \space \space \space \space \space {\color{Red}公式7:} \space \space \space S=\sum_{i=0}^{n-1}\frac{(y_{i-1} + y_{i})×(x_{i-1} - x_{i})}{2} 公式7: S=i=0∑n−12(yi−1+yi)×(xi−1−xi)
对于整个多边形都在x轴下方的情况,上述公式因为y变成负值了,公式6会使下面的边形成梯形面积变成正值,上面边变成负值,计算结果仍然正确。
同理可验证,x轴穿过多边形,公式7计算出的面积也是正确的,自己简单画一下就清楚了。
对于顺时针多边形,也很容易验证,公式7得到的结果是多边形面积的相反数,取绝对值就是面积了。
同样,这个公式计算面积的同时,也可以得到多边形的顺逆时针方向。
对于凹多边形公式7仍然成立。和凸多边形的区别是,凸多边形被梯形覆盖区域的多边形外面积被加了1次又减了1次;凹多边形是加了n次又减了n次。
梯形累加这个算法要比分割三角形算法有更广的应用范围。因为它本质上就是对多边形求积分。当多边形的边包含曲线时(如下图所示),这个方法依然可用,我们只需要将用公式6计算面积改为数值积分计算面积即可。当然,如果非要说割三角形那个也是积分我也认可,只能说在积分这个角度理解梯形更直观。