本文介绍了本三角剖分算法数学运算的核心内容——计算射线与线段的交点,通过引入新的数学工具大大简化了该运算,并利用代码和图示进行了详细阐述。
上篇文章介绍了步骤显示功能的具体实现方式,并且给出了我对面向对象设计的一些建议。本篇将会介绍在本三角剖分算法中涉及到的数学知识,并向您展示在引入新的数学工具后实现的便捷。
新的算法由寻找辅助线、寻找相交点、寻找三角形三个大块组成,这几大块的大致思路已经在《第二回:漫谈新思路,是我们自己干的时候了》中有所交代。其中,寻找相交点需要进行大量的数学运算,而关键的数学运算是计算射线与线段的交点。
一、老工具
计算射线与线段的交点可以有很多方法,比如使用解析几何。我大概的描述一下使用解析几何的做法:先找到射线与线段所在的两条直线,得到两条直线的函数,联立为二元一次方程组,解出结果得到两条直线的交点,再判断交点是否在射线上以及交点是否在线段上。
这种做法原不原始尚不讨论,单就其运算的繁复就已经大打折扣,当然使用这种做法有一个好处:一个耐心的高中生就可搞定,省却了看这篇文章的时间
二、新工具
隆重推出新的数学工具:向量。我好似看到一堆鸡蛋、西红柿朝我飞来,向量还算新?高中就学过!没错,但是在这里我会引入一些新的用法,从而使这种数学工具焕发新的青春。
在这里我假设大家已经了解向量的概念,基本的运算法则,在以后的系列中我会详细讲解向量。
由于是解决平面三角剖分问题,所以使用的向量是二维的,因而这里不会涉及到向量叉乘的问题,但是点乘会经常涉及,垂直向量点乘为零的性质也会经常用到。
在本篇文章中粗体小写字母代表的是向量,斜体字母代表的是标量,大写字母代表的是点。下面是二维向量类和二维点类的代码:
Vector2
1Public Class Vector2Class Vector2
2 Private xx As Decimal '由于需要精度较高,这里使用Decimal类
3 Private yy As Decimal
4
5 Public Sub New()Sub New(ByVal a As Decimal, ByVal b As Decimal)
6 Me.xx = a
7 Me.yy = b
8 End Sub
9
10 Public Property X()Property X() As Decimal
11 Get
12 Return Me.xx
13 End Get
14 Set(ByVal value As Decimal)
15 Me.xx = value
16 End Set
17 End Property
18
19 Public Property Y()Property Y() As Decimal
20 Get
21 Return Me.yy
22 End Get
23 Set(ByVal value As Decimal)
24 Me.yy = value
25 End Set
26 End Property
27
28 Public ReadOnly Property K()Property K() As Double '斜率
29 Get
30 If (Me.xx = 0) Then '如果是竖直方向上的向量
31 If (Me.yy > 0) Then '如果指向上
32 Return Double.PositiveInfinity '返回正无穷
33 ElseIf (Me.yy < 0) Then '如果指向下
34 Return Double.NegativeInfinity '返回负无穷
35 Else '如果是零向量
36 Return 0 '返回零,代表水平向量
37 End If
38 End If
39 Return Me.yy / Me.xx
40 End Get
41 End Property
42
43 Public ReadOnly Property Angel()Property Angel() As Decimal '与水平线夹角
44 Get
45 Return Math.Atan(Me.K)
46 End Get
47 End Property
48
49 Public ReadOnly Property Length()Property Length() As Decimal '长度
50 Get
51 Return Math.Sqrt(X * X + Y * Y)
52 End Get
53 End Property
54
55 Public Function perp()Function perp() As Vector2 '返回与该矢量长度相同,逆时针旋转至垂直的向量
56 Return New Vector2(-Me.Y, Me.X)
57 End Function
58
59 Public Overrides Function ToString()Function ToString() As String
60 Return xx & "," & yy
61 End Function
62
63 Public Shared Operator +(ByVal vector As Vector2, ByVal point As Point2) As Point2
64 Return New Point2(point.X + vector.X, point.Y + vector.Y)
65 End Operator
66
67 Public Shared Operator +(ByVal vector1 As Vector2, ByVal vector2 As Vector2) As Vector2
68 Return New Vector2(vector1.X + vector2.X, vector1.Y + vector2.Y)
69 End Operator
70
71 Public Shared Operator *(ByVal scaler As Decimal, ByVal vector As Vector2) As Vector2
72 Return New Vector2(scaler * vector.X, scaler * vector.Y)
73 End Operator
74
75 Public Shared Operator *(ByVal vector As Vector2, ByVal scaler As Decimal) As Vector2
76 Return New Vector2(scaler * vector.X, scaler * vector.Y)
77 End Operator
78
79 Public Shared Operator *(ByVal vector1 As Vector2, ByVal vector2 As Vector2) As Decimal
80 Return vector1.X * vector2.X + vector1.Y * vector2.Y
81 End Operator
82
83
84End Class
85
Point2
1Public Class Point2Class Point2
2 Private xx As Decimal
3 Private yy As Decimal
4
5 Public Sub New()Sub New(ByVal x As Decimal, ByVal y As Decimal)
6 Me.xx = x
7 Me.yy = y
8 End Sub
9
10 Public Property X()Property X() As Decimal
11 Get
12 Return Me.xx
13 End Get
14 Set(ByVal value As Decimal)
15 Me.xx = value
16 End Set
17 End Property
18
19 Public Property Y()Property Y() As Decimal
20 Get
21 Return Me.yy
22 End Get
23 Set(ByVal value As Decimal)
24 Me.yy = value
25 End Set
26 End Property
27
28 Public Overrides Function Equals()Function Equals(ByVal obj As Object) As Boolean
29 If (Not (TypeOf obj Is Point2)) Then
30 Return False
31 End If
32 Dim point As Point2 = CType(obj, Point2)
33 Dim vector As Vector2 = point - Me
34 If (vector.Length > 0.000001D) Then
35 Return False
36 End If
37 Return True
38 End Function
39
40 Public Overrides Function ToString()Function ToString() As String
41 Return xx & "," & yy
42 End Function
43
44 Public Shared Operator -(ByVal point1 As Point2, ByVal point2 As Point2) As Vector2
45 Return New Vector2(point1.X - point2.X, point1.Y - point2.Y)
46 End Operator
47
48 Public Shared Operator +(ByVal point As Point2, ByVal vector As Vector2) As Point2
49 Return New Point2(point.X + vector.X, point.Y + vector.Y)
50 End Operator
51
52End Class
53
在上面的代码中,重载了一系列的运算符,各自功用不难理解,不过是为了方便编程罢了。Point2的Equals方法不是简单的比较两个点的坐标,而是将待比较的两个点连成一个向量,根据向量的长短来进行比较,当向量足够短时就认定两个点相同。Vector2的perp方法返回一个与该向量长度相同,逆时针旋转至垂直的新向量,说着挺玄,看看代码其实就是简单的将横纵坐标交换,然后将原来的y反向,别看它的定义如此简单,后面可有大用途,在我看来该方法的引入使向量大大进化,运算方法更加灵活,向引入该方法的数学家(可惜不知道是谁)致敬!在本篇文章中a的perp以符号~a表示。
三、用向量描述直线
大家知道,向量是没有位置的,所以要表现一条直线需要一个点和一个向量,而点和向量所决定的直线为A + ct ,其中t的取值为(-∞,+∞),当取其中某个值时确定一个点,特别的,当t为零时,为向量的起始点A,当t为1时,为向量的终止点A + c,当t小于零时为向量反方向上的某个点,详情如下图
你是否有什么惊奇的发现,看到这个t的取值图有没有眉头微释,没错,得到了t你就找到了点,而t落在的区间就决定了它与射线、线段的关系:t>0时点落在射线上,t∈[0,1]时点落在线段上,哈哈,就是这么简单,当你找到了点,它的位置就尽在掌握,岂一个爽字了得。
四、射线与线段的交点
上面给出了向量描述直线的方法,求射线与线段的交点实质上就是求它们所在直线的交点,显然的,当A + ct = M + eu时正是交点,如下图:
整理一下上面的等式,得到ct = (M–A) + eu(一),这里的M–A就是一个从M开始到A结束的向量,当然由于上面的代码中已经给出了二维点的减号重载,程序里直接将两个点相减就可以了。
下面的任务就是求出t和u,继而得到交点了,这里需要一些技巧,还记得上文中提到的perp吗?
在上述等式的两边同时乘以e的perp,也即~e,由于(~e)•e = 0,所以上面的等式就恒等变形为ct•(~e) = (M–A)•(~e),继而t = (M–A)• (~e) / c•(~e),这些都是已知量了,t也就求出来了,将t代入(一)式,u也即求出,当然也可如法炮制在(一)式两边同时乘以个~c。t和u都求出后,发现t>1,故而两直线的交点落在了蓝色射线上,但是没有落在AB线段上,而0<u<1,故而该交点落在了线段MN上,怎么样,简单!
另外,当考虑某条射线和一系列线段相交时将求出的各个t记录下来,经过将t排序可以按照射线的方向将各个交点排序,这在《第七回:寻找三角形,夺取红旗》中会用到。
射线与线段的交点通过引入向量这一数学工具用一种较为简便的方式解决了,在程序中只需三行就可以得到t和u,继而找到交点的位置,配之操作符的种种重载,代码犹如数学公式,一目了然,优雅之至!在《第六回:寻找交点,离胜利就剩一步 之 开找》中你将见到这段代码。向量知识博大而精深,这篇小文只是九牛之一毛,但对于该算法已经够用了。下一篇将详细讲述数据结构的设计,用一种较为合适的方式来存储运算过程中得到的各种数据以及最后的运算结果。