计算几何常用算法及numpy仿真

在很久之前的一篇的文章点乘和叉乘及其物理意义（C++STL实现），我们用C++（STL）实现了对向量内积和叉积的定义与简单计算，最后演示了如何用几何的方法计算点到直线的距离，计算任意三角形的面积等问题。把这些放在更大的范围内，其实就是今天的主角，计算几何（Computational geometry）。

本文所有的仿真程序均在numpy（python第三方的科学计算库）下进行：

import numpy as np

1. 矢量减法

设二维矢量 $P =(x_1,y_1)$ ，$Q =(x_2,y_2)$
则矢量减法定义（对应位相减）为： $P - Q = ( x_1 - x_2 , y_1 - y_2 )$
显然有性质 $P - Q = - ( Q - P )$
如不加说明，下面所有的点都看作矢量，两点的减法就是矢量相减；

def sub(p, q):
    return p-q
p, q = np.array([1, 2]), np.array([3, 5])
print(sub(p, q))
print(-sub(q, p))

2. 矢量叉积

设矢量$P = (x_1,y_1), Q = (x_2,y_2)$（本文以二维为例，可自然推广到三维）
则矢量叉积定义为： $P × Q = x_1*y_2 - x_2*y_1$ 得到的是一个标量
显然有性质 $P × Q = - ( Q × P ), P × ( - Q ) = - ( P × Q )$
如不加说明，下面所有的点都看作矢量，点的乘法看作矢量叉积；

叉乘的重要性质：

若 P × Q > 0 , 则P 在Q的顺时针方向
若 P × Q < 0 , 则P 在Q的逆时针方向（与性质1等价）
若 P × Q = 0 , 则P 与Q共线，但可能同向也可能反向

def cross_prod(p, q):
    assert len(p)==2 and len(q)==2
    return p[0]*q[1]-p[1]*q[0]
p, q = np.array([1, 0]), np.array([0, 1])
                        # p在q的顺时针方向
print(cross_prod(p, q))
                        # 1>0
p, q = np.array([1, 2], [2, 4])         
                        # p, q 同向
print(cross_prod(p, q))
                        # 0

3. 判断点在线段上

设点为$Q$，线段为$P_1P_2$ ，判断点$Q$在该线段上的依据是：
$(Q - P_1) × ( P_2 - P_1 ) = 0$ 且 $Q$ 在以 $P_1,P_2$ 为对角顶点的矩形内；
这是叉积性质的直接应用，也即$(Q-P_1)\times (P_2-P_1)=0$，则 $Q-P_1$与$P_2-P_1$共线；

def is_point_on_line(q, p1, p2):
    return True if cross_prod(q-p1, p2-p1)==0 else False
p1, p2 = np.array([1, 0]), np.array([0, 1])
q = np.array([1/2, 1/2])
print(is_point_on_line(q, p1, p2))
                    # True
q = np.array([1/2, 1/3])
print(is_point_on_line(q, p1, p2))
                    # False

4. 判断两线段是否相交

我们分两步确定两条线段是否相交：

1. 快速排斥试验

   设以线段 $P_1P_2$ 为对角线的矩形为$R$， 设以线段 $Q_1Q_2$ 为对角线的矩形为$T$，如果$R$和$T$不相交，显然两线段不会相交；

def rect_overlap(r1, r2):
    return not (((r1[1][0] < r2[0][0])|(r1[0][1] > r2[1][1]))
               |((r2[1][0] < r1[0][0])|(r2[0][1] > r1[1][1]))

1. 跨立试验

   如果两线段相交，则两线段必然相互跨立对方，如图1所示。在图1中，$P_1P_2$跨立$Q_1Q_2$ ，则 矢量 $( P1 - Q1 )$ 和$(P_2 - Q_1)$位于矢量$(Q_2 - Q_1)$ 的两侧，即$(P1-Q1)×(Q2-Q1) * (P2-Q1 )×(Q2-Q1)<0$ 上式可改写成 ( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) > 0，
   当 ( P1 - Q1 ) × ( Q2 - Q1 ) = 0 时，说明 ( P1 - Q1 ) 和 ( Q2 - Q1 )共线，但是因为已经通过快速排斥试验，所以 P1 一定在线段 Q1Q2上；同理，( Q2 - Q1 ) ×(P2 - Q1 ) = 0 说明 P2 一定在线段 Q1Q2上。 所以判断P1P2跨立Q1Q2的依据是： ( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) ≥ 0 同理判断Q1Q2跨立P1P2的依据是： ( Q1 - P1 ) × ( P2 - P1 ) * ( P2 - P1 ) × ( Q2 - P1 ) ≥ 0至此已经完全解决判断线段是否相交的问题。

def line_intersection(p1, p2, q1, q2):
    return False if not rect_overlap([p1, p2], [q1, q2]) else \             # 不相交直接返回
            (True if cross_prod(p1-q1, q2-q1)*cross_prod(q2-q1, p2-q1) >= 0 else False)

5. 判断线段和直线是否相交

从标题即可看出，这里的直线和线段是不相同的；
如果线段 $P_1P_2$和直线 $Q_1Q_2$相交，则 $P_1P_2$跨立 $Q_1Q_2$，根据4的结论：

$$(P_1-Q_1)\times(Q_2-Q_1) \;*\; (Q_2-Q_1)\times(P_2,Q_1)\beq 0$$

6. 判断点是否在矩形内

只需判断该点的横纵坐标是否夹在矩形的左右边，上下边之间。
判断线段、折线、多边形是否在矩形中，因为矩形是个凸集，所以只要判断所有端点是否都在矩形中即可；

7. 判断矩形是否在矩形中

只要把比较左右边界和上下边界即可；

8. 判断圆是否在矩形中：

圆在矩形内的充要条件：
圆心在矩形中且圆的半径小于等于圆心到矩形四条边的距离的最小值。

def circle_in_rect(c, r, lb, rt):
                # c: center
                # r: radius
                # lb: left bottom
                # rt: right top
    if ( ((center[0]<lb[0]) | (center[0]>rt[0])) | ((center[1]<lb[1])|(center[1]>rt[1])) ):
                # 圆心是否在矩形内
        return False
    return True r <= min(center[0]-lb[0], rt[0]-center[0], center[1]-lb[1], rt[1]-center[1]) if else False

9. 判断点是否在多边形中

References

[1] 计算几何常用算法