PKUWWT

消遣时光

两个平面三角形的相交测试

平面三角形求交

下面给出两个平面三角形的求交测试。测试过程分两步,第一步检测两个三角形边与边之间是否相交(9次),第二步检测一个三角形是否完全在一个三角形内部。

线段相交

检测线段是否相交是一个简单的解析几何问题。设第一条线段为AB, 起始坐标分别为,第二线段为CD,起始坐标分别为,交点为$E$。令。则有

因此,可以求得两个系数

只需要求即可。具体计算时,并不需要真正地把系数求出来,只需要计算出分子分母三个数,然后,比较符号和大小即可。如果分母为0,说明两条线段平行,如果分子为0,表明两线段交于某一端点。此算法中主要涉及到6次乘法和若干次加法。python代码如下

def line_intersect2(v1,v2,v3,v4):
    '''
    judge if line (v1,v2) intersects with line(v3,v4)
    '''
    d = (v4[1]-v3[1])*(v2[0]-v1[0])-(v4[0]-v3[0])*(v2[1]-v1[1])
    u = (v4[0]-v3[0])*(v1[1]-v3[1])-(v4[1]-v3[1])*(v1[0]-v3[0])
    v = (v2[0]-v1[0])*(v1[1]-v3[1])-(v2[1]-v1[1])*(v1[0]-v3[0])
    if d<0:
        u,v,d= -u,-v,-d
    return (0<=u<=d) and (0<=v<=d)

三角形内部判定

判别点P是否在三角形ABC内部,最直观的想法是绕三角形逆时针旋转时,判断此点是否在三角形所有边的左侧即可。但这个方法计算量较大,更好的办法是使用重心坐标系(实际上新坐标系原点不在三角形重心上)。具体而言,以A为坐标原点,AC和AB分别坐标系的两个轴(这不是直角坐标系),只要它们不重合,就可以为空间中任意一点指定一个坐标。在AC轴上的0坐标对应着A点,1坐标对应着C点。同理,AB轴上0对应A,1对应B。另一方面,两个坐标值之和如果都大于等于0,且和小于1,则此点位于三角形内部(或边上)。

那么,现在的问题就变成了,如何求重心坐标。它实际上是将矢量AP分解为平行于AB和AC的两个矢量。下面是用D3写的一个演示图。

,则构成一个线性方程组

结果为

此算法同样只涉及到6次乘法和若干次加法。python代码如下

def point_in_triangle2(A,B,C,P):
    v0 = [C[0]-A[0], C[1]-A[1]]
    v1 = [B[0]-A[0], B[1]-A[1]]
    v2 = [P[0]-A[0], P[1]-A[1]]
    cross = lambda u,v: u[0]*v[1]-u[1]*v[0]
    u = cross(v2,v0)
    v = cross(v1,v2)
    d = cross(v1,v0)
    if d<0:
        u,v,d = -u,-v,-d
    return u>=0 and v>=0 and (u+v) <= d

最后总的三角形相交测试代码为

def tri_intersect2(t1, t2):
    '''
    judge if two triangles in a plane intersect 
    '''
    if line_intersect2(t1[0],t1[1],t2[0],t2[1]): return True
    if line_intersect2(t1[0],t1[1],t2[0],t2[2]): return True
    if line_intersect2(t1[0],t1[1],t2[1],t2[2]): return True
    if line_intersect2(t1[0],t1[2],t2[0],t2[1]): return True
    if line_intersect2(t1[0],t1[2],t2[0],t2[2]): return True
    if line_intersect2(t1[0],t1[2],t2[1],t2[2]): return True
    if line_intersect2(t1[1],t1[2],t2[0],t2[1]): return True
    if line_intersect2(t1[1],t1[2],t2[0],t2[2]): return True
    if line_intersect2(t1[1],t1[2],t2[1],t2[2]): return True
    inTri = False
    inTri = inTri or point_in_triangle2(t1[0],t1[1],t1[2], t2[0])
    inTri = inTri or point_in_triangle2(t1[0],t1[1],t1[2], t2[1])
    inTri = inTri or point_in_triangle2(t1[0],t1[1],t1[2], t2[2])
    if inTri == True: return True
    inTri = False
    inTri = inTri or point_in_triangle2(t2[0],t2[1],t2[2], t1[0])
    inTri = inTri or point_in_triangle2(t2[0],t2[1],t2[2], t1[1])
    inTri = inTri or point_in_triangle2(t2[0],t2[1],t2[2], t1[2])
    if inTri == True: return True
    return False

最后给出一个测试程序,程序中有两个可拖动的三角形,如果三角形相交,则会在标题栏上显示’intersect’,否则显示’non-intersect’。

#!/usr/bin/python
import numpy as np
from numpy.random import random
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon

def line_intersect2(v1,v2,v3,v4):
    '''
    judge if line (v1,v2) intersects with line(v3,v4)
    http://thirdpartyninjas.com/blog/2008/10/07/line-segment-intersection/
    '''
    v12 = [v2[0]-v1[0],(v2[1]-v1[1])]
    v34 = [v4[0]-v3[0],(v4[1]-v3[1])]
    v31 = [v1[0]-v3[0],(v1[1]-v3[1])]
    cross = lambda u,v: u[0]*v[1]-u[1]*v[0]
    d = cross(v12,v34)
    u = cross(v34,v31)
    v = cross(v12,v31)
    if d<0:
        u,v,d= -u,-v,-d
    return (0<=u<=d) and (0<=v<=d)

def point_in_triangle2(A,B,C,P):
    v0 = [C[0]-A[0], C[1]-A[1]]
    v1 = [B[0]-A[0], B[1]-A[1]]
    v2 = [P[0]-A[0], P[1]-A[1]]
    cross = lambda u,v: u[0]*v[1]-u[1]*v[0]
    u = cross(v2,v0)
    v = cross(v1,v2)
    d = cross(v1,v0)
    if d<0:
        u,v,d = -u,-v,-d
    return u>=0 and v>=0 and (u+v) <= d


def tri_intersect2(t1, t2):
    '''
    judge if two triangles in a plane intersect 

    e.g.
        tri_intersect2([(0,0),(1,0),(0,1)], [(1,0),(2,0),(1,1)])

    '''
    if line_intersect2(t1[0],t1[1],t2[0],t2[1]): return True
    if line_intersect2(t1[0],t1[1],t2[0],t2[2]): return True
    if line_intersect2(t1[0],t1[1],t2[1],t2[2]): return True
    if line_intersect2(t1[0],t1[2],t2[0],t2[1]): return True
    if line_intersect2(t1[0],t1[2],t2[0],t2[2]): return True
    if line_intersect2(t1[0],t1[2],t2[1],t2[2]): return True
    if line_intersect2(t1[1],t1[2],t2[0],t2[1]): return True
    if line_intersect2(t1[1],t1[2],t2[0],t2[2]): return True
    if line_intersect2(t1[1],t1[2],t2[1],t2[2]): return True
    inTri = True 
    inTri = inTri and point_in_triangle2(t1[0],t1[1],t1[2], t2[0])
    inTri = inTri and point_in_triangle2(t1[0],t1[1],t1[2], t2[1])
    inTri = inTri and point_in_triangle2(t1[0],t1[1],t1[2], t2[2])
    if inTri == True: return True
    inTri = True
    inTri = inTri and point_in_triangle2(t2[0],t2[1],t2[2], t1[0])
    inTri = inTri and point_in_triangle2(t2[0],t2[1],t2[2], t1[1])
    inTri = inTri and point_in_triangle2(t2[0],t2[1],t2[2], t1[2])
    if inTri == True: return True
    return False

class DraggableTriangles:
    def __init__(self,tri):
        self.tri = tri
        self.center0 = tri.get_xy().mean(0)
        self.offset = np.zeros((2,))
        self.press = None
    def connect(self):
        'connect to all events we need'
        self.cidpress = self.tri.figure.canvas.mpl_connect(
                'button_press_event', self.on_press)
        self.cidrelease = self.tri.figure.canvas.mpl_connect(
                'button_release_event', self.on_release)
        self.cidmotion = self.tri.figure.canvas.mpl_connect(
                'motion_notify_event', self.on_motion)
    def on_press(self,event):
        if event.inaxes != self.tri.axes: return
        contains, attrd = self.tri.contains(event)
        if not contains: return
        self.press = event.xdata,event.ydata
    def on_motion(self,event):
        if self.press is None: return
        if event.inaxes != self.tri.axes: return
        xpress, ypress = self.press
        dx = event.xdata - xpress + self.offset[0]
        dy = event.ydata - ypress + self.offset[1]
        trans = mpl.transforms.Affine2D().translate(dx,dy)+self.tri.axes.transData
        self.tri.set_transform(trans)
        self.tri.figure.canvas.draw()
    def on_release(self,event):
        if self.press is None: return
        xpress, ypress = self.press
        dx = event.xdata - xpress + self.offset[0]
        dy = event.ydata - ypress + self.offset[1]
        self.offset = (dx,dy)
        self.press = None
        self.tri.figure.canvas.draw()
    def disconnect(self):
        self.tri.figure.canvas.mpl_disconnect(self.cidpress)
        self.tri.figure.canvas.mpl_disconnect(self.cidrelease)
        self.tri.figure.canvas.mpl_disconnect(self.cidmotion)

fig = plt.figure()
ax = fig.add_subplot(111)
dtris = []
tris = [Polygon(random((3,2))*4) for i in range(2)]

def on_motion(event):
    t1 = tris[0].get_verts()
    t2 = tris[1].get_verts()
    if tri_intersect2(t1[0:3],t2[0:3]):
        ax.set_title('intersect')
    else:
        ax.set_title('non-intersect')

fig.canvas.mpl_connect('motion_notify_event', on_motion)


for tri in tris:
    tri.set_color((random(),random(),random()))
    ax.add_patch(tri)
    dtri = DraggableTriangles(tri)
    dtri.connect()
    dtris.append(dtri)

ax.autoscale_view()
plt.show()

参考文献

  • http://thirdpartyninjas.com/blog/2008/10/07/line-segment-intersection/
  • http://www.blackpawn.com/texts/pointinpoly/default.html

注意,本文中的三角形内部判定算法要优于参考文献中的算法,不过两者的基本原理是一致的。