多边形面积计算及顶点顺-逆时针方向判断

格林（Green）公式告诉我们，在平面闭区域D上的二重积分可以通过沿闭区域D的边界曲线L上的曲线积分来表达。即，设闭区域 $D$ 由分段光滑的曲线 $L$ 围成，函数 $P(x,y)$ 及 $Q(x,y)$ 在 $D$ 上具有一阶连续偏导数，则有 $\iint_{D}(\frac{\partial Q}{\partial x}-\frac{\partial P}{\partial y})dxdy=\oint_{L}Pdx+Qdy$ 其中 $L$ 是 $D$ 的取正向的边界曲线。根据公式推导过程（参考高等数学同济版下），有如下关系存在： $-\iint_{D}\frac{\partial P}{\partial y}dxdy=\oint_{L}Pdx\\\iint_{D}\frac{\partial Q}{\partial x}dxdy=\oint_{L}Qdy$

对于上面两个式子，分别取 $P=-y$ ， $Q=x$ 可得到： $\begin{align}&\iint_{D}dxdy=A=-\oint_{L}ydx\\&\iint_{D}dxdy=A=\oint_{L}xdy\end{align}$

其中， $A$ 为闭区域 $D$ 的面积。根据式(1)和(2)也可得出： $A=\frac{1}{2}\oint_{L}xdy-ydx$

根据公式(2)来计算多边形面积，参考下面的示意图，曲线 $L$ 由多段光滑的直线段 $L_i$ 连接而成。由于 $L$ 是分段光滑的，则有向曲线 $L$ 上对坐标的曲线积分等于在光滑的各段上对坐标的曲线积分之和。即 $A=\oint_Lxdy=\int_{L_0}xdy+...+\int_{L_6}xdy$

直线段 $L_i$ 从点 $(x_i,y_i)$ 到点 $(x_{i+1},y_{i+1})$ 的参数方程为 $L_i(t)=\left ((x_{i+1}-x_i)t+x_i,(y_{i+1}-y_i)t+y_i\right),\quad 0\leq t\leq 1$ ，于是对线段 $L_i$ 的坐标积分为： $\begin{aligned}\int_{Li}xdy&=\int_0^1((x_{i+1}-x_i)t+x_i)(y_{i+1}-y_i)dt\\&=(y_{i+1}-y_i)(x_it+\frac{1}{2}t^2(x_{i+1}-x_i))\big|_{0}^{1}\\&=\frac{1}{2}(x_{i+1}+x_i)(y_{i+1}-y_i)\end{aligned}$

将每段结果相加得到多边形的面积为： $\boxed{A=\sum_{i=0}^{n}\frac{(x_{i+1}+x_i)(y_{i+1}-y_i)}{2}}$

其中， $n$ 为多边形顶点数-1， $(x_{n+1},y_{n+1})=(x_0,y_0)$ 。类似的，根据公式(1)可计算出多边形的面积为： $\boxed{A=\sum_{i=0}^{n}\frac{(x_i-x_{i+1})(y_i+y_{i+1})}{2}}$

多边形面积和顶点数组的顺-逆时针方向可以根据上面的公式计算，即如果是逆时针（ $L$ 为正向），则计算出的面积为正值；如果是顺时针，计算出的面积为负值。可参考下面的代码：

class Point:
    def __init__(self, x, y):
        self.x = x
        self.y = y

def isCounterclockwise(points):
    s = 0
    points_count = len(points)
    for i in range(points_count):
        point = points[i]
        point2 = points[(i + 1) % points_count]
        s += (point.x - point2.x) * (point.y + point2.y)
    return s > 0


def poly_area(points):
    s = 0
    points_count = len(points)
    for i in range(points_count):
        point = points[i]
        point2 = points[(i + 1) % points_count]
        s += (point.x - point2.x) * (point.y + point2.y)
    return abs(s/2)


if __name__ == "__main__":
    polygon1 = [Point(0,0), Point(0.5,0), Point(0.5,0.5), Point(0,1)] # counter-clockwise
    polygon2 = polygon1[::-1]  # clockwise
    print(isCounterclockwise(polygon1), isCounterclockwise(polygon2)) # True, False
    print(poly_area(polygon1)) # 0.375