浏览量:262

4.2. 点与三角形

本篇文章禁止用于任何商业目的,版权申明、版本说明等见《前言》。

PDF文档和源码下载地址:https://github.com/twinklingstar20/Programmers_Computational_Geometry

4.2. 点与三角形

4.2.1. 点与三角形的关系

1. 二维空间

问题描述:在二维空间上,判断点\(P\)与三角形\({V_0}{V_1}{V_2}\)的关系,其中,\({V_i} = \left( {{x_i},{y_i}} \right),i \in \{ 0,1,2\} \)

2015-3-23 1-17-48

图4.3 点与三角形的三种关系

二维空间上点与三角形共有三种关系:三角形内,三角形边上和三角形外。如图4.3所示,\({P_0}\)在三角形内部,\({P_1}\)在三角形边上,\({P_2}\)在三角形外。

在二维空间上的叉积遵守右手定则,设有三个点\({B_0}\),\({B_1}\)\({B_2}\),令\(d = \left( {{B_1} - {B_0}} \right) \times \left( {{B_2} - {B_0}} \right)\),如果\(d \succ 0\),则\({B_0}\),\({B_1}\)\({B_2}\)三个点是逆时针方向分布;如果\(d \prec 0\),则是顺时针方向分布;无否则,有\(d = 0\),则三点是共线的。

如果点\(P\)\(\Delta {V_0}{V_1}{V_2}\)内,则\(\left( {{V_0},{V_1},P} \right),\left( {{V_1},{V_2},P} \right),\left( {{V_2},{V_0},P} \right)\)三组点都是顺时针方向或者逆时针方向,令

\[\left\{ \begin{array}{l}{d_0} = \left( {{V_1} - {V_0}} \right) \times \left( {P - {V_0}} \right)\\{d_1} = \left( {{V_2} - {V_1}} \right) \times \left( {P - {V_1}} \right)\\{d_2} = \left( {{V_0} - {V_2}} \right) \times \left( {P - {V_2}} \right)\end{array} \right. \tag{4.4}\]

如果\({d_0}\),\({d_1}\)\({d_2}\)既存在正值也存在负值,则点在三角形外;如果\({d_0}\),\({d_1}\)\({d_2}\)至少有一个值为零且其它值全是正或者都是负,则点在三角形的边上;否则,\({d_0}\),\({d_1}\)\({d_2}\)全是正或者全是负,则点在三角形内部。

算法伪码如下所示:

在二维空间上,判断点\(P\)\(\Delta {V_0}{V_1}{V_2}\)的关系,其中,\({V_i} = ({x_i},{y_i}),i \in \{ 0,1,2\} \)
return ({OUSIDE, INSIDE, ON_EDGE})
1.    \({d_0} = \left( {{V_1} - {V_0}} \right) \times \left( {P - {V_0}} \right)\);
2.    \({d_1} = \left( {{V_2} - {V_1}} \right) \times \left( {P - {V_1}} \right)\);
3.    \({d_2} = \left( {{V_0} - {V_2}} \right) \times \left( {P - {V_2}} \right)\);
4.    \({s_0} = {d_0} \cdot {d_1},{s_1} = {d_1} \cdot {d_2}\);
5.    if\({s_0} \prec 0\)或者\({s_1} \prec 0\), then
6.           return OUSIDE;
7.    else\({s_0} = 0\)或者\({s_1} = 0\), then
8.           return ON_EDGE;
9.    else
10.          return INSIDE;
11.   end if

2. 三维空间

问题描述:在三维空间上,判断点\(P\)\(\Delta {V_0}{V_1}{V_2}\)的关系,其中,三角形所在的平面是\(\pi \)\({V_i} = ({x_i},{y_i},{z_i}),i \in \{ 0,1,2\} \)

2015-3-23 1-18-00

图4.4 点与三角形的四种关系

在三维空间上,点与三角形之间共有四种关系:点在平面\(\pi \)外,点在三角形内,点在三角形的边上,点在平面\(\pi \)上但在三角形外。如图4.4所示,点\({P_0}\)位于平面\(\pi \)外,点\({P_1}\)位于三角形内,点\({P_2}\)在三角形的边上,点\({P_3}\)在平面\(\pi \)上但在三角形外。

有两种算法来判定点与三维三角形的关系,第一种算法是把点\(P\)和三角形沿着坐标轴投影到二维的平面上,将问题简化为在二维空间上判定点和三角形的关系;第二种算法是基于平面法向量的方法。这里介绍第二种算法。

2015-3-23 1-18-14

图4.5 三角形上的边所在的直线把平面分割为两个半平面

如图4.5所示,平面\(\pi \)可以表示为\(\vec n \cdot X + d = 0\),由三角形\(T\)的三个顶点,可以得到平面\(\pi \)的参数表示,即\(\vec n = ({V_1} - {V_0}) \times ({V_2} - {V_0}),d = - \vec n \cdot {V_i},i \in \{ 0,1,2\} \)

首先,判断点\(P\)是否在平面\(\pi \)上。计算点\(P\)到平面的距离,即\({\mathop{\rm dist}\nolimits} (P,\pi ) = \left| {\vec n \cdot P + d} \right|/\left\| {\vec n} \right\|\)。这里只需要判断距离是否为零,因此可以避免除法操作,即若\(\left| {\vec n \cdot P + d} \right| = 0\),则点\(P\)在平面上,否则,在平面外。

接着,考虑点\(P\)在平面\(\pi \)上的情况。以边\(\overline {{V_0}{V_1}} \)为例,它所在的直线把平面\(\pi \)分为两个半平面,如图4.5所示,点\({P_0}\)\({V_2}\)不在同一个半平面上,则\({P_0}\)一定在三角形外;否则,\(P\)可能在三角形边上,也可能在三角形外或三角形内,如图中\({P_1},{P_2},{P_3}\)所示。设\(\overline {{V_i}{V_j}} \)把平面分成两个半平面,判断点\(P\)\({V_k}\)是否在同一个半平面可以采用如下等式:

\[f\left( {\overline {{V_i}{V_j}} ,Q,P} \right) = \left[ {({V_j} - {V_i}) \times (Q - {V_i})} \right] \cdot \left[ {({V_j} - {V_i}) \times (P - {V_i})} \right] \tag{4.5}\]

如果\(f\left( {\overline {{V_i}{V_j}} ,P,Q} \right) \succ 0\),表示点\(P\)\(Q\)在同一个半平面上;如果\(f\left( {\overline {{V_i}{V_j}} ,P,Q} \right) \prec 0\),表示不在同一个半平面上;否则,点\(P\)在边\(\overline {{V_i}{V_j}} \)所在的直线上。

\[\left\{ \begin{array}{l}{d_0} = f\left( {\overline {{V_0}{V_1}} ,{V_2},P} \right)\\{d_1} = f\left( {\overline {{V_1}{V_2}} ,{V_0},P} \right)\\{d_2} = f\left( {\overline {{V_2}{V_0}} ,{V_1},P} \right)\end{array} \right. \tag{4.6}\]

那么,点\(P\)在三角形内当且仅当\({d_0} \succ 0,{d_1} \succ 0,{d_2} \succ 0\)。如果对于所有的\({d_i},i \in \{ 0,1,2\} \)满足\({d_i} \ge 0\),且至少有一个\({d_i} = 0\),则点\(P\)在三角形的边上;否则,说明至少存在一个\({d_i} \prec 0\),则点\(P\)在三角形外。

又由叉积的性质,易知\(({V_1} - {V_0}) \times ({V_2} - {V_0}) = ({V_2} - {V_1}) \times ({V_0} - {V_1}) = ({V_0} - {V_2}) \times ({V_1} - {V_2}) = \vec n\),所以等式(4.6)可以简化:

\[\left\{ \begin{array}{l}{d_0} = \vec n \cdot \left[ {\left( {{V_1} - {V_0}} \right) \times \left( {P - {V_0}} \right)} \right]\\{d_1} = \vec n \cdot \left[ {\left( {{V_2} - {V_1}} \right) \times \left( {P - {V_1}} \right)} \right]\\{d_2} = \vec n \cdot \left[ {\left( {{V_0} - {V_2}} \right) \times \left( {P - {V_2}} \right)} \right]\end{array} \right. \tag{4.7}\]

如图4.5所示,显示了平面\(\pi \)上不同的点相对于三角形\({V_0}{V_1}{V_2}\)\({d_0},{d_1},{d_2}\)的值与它们和三角形的关系。

算法的伪码如下所示:

在三维空间上,判断点\(P\)与三角形\(\Delta {V_0}{V_1}{V_2}\)的关系,其中,三角形所在的平面是\(\pi \)\({V_i} = ({x_i},{y_i},{z_i}),i \in \{ 0,1,2\} \)
return ({NOT_ON_PLANE, OUSIDE, INSIDE, ON_EDGE})
1.    \(\vec n = ({V_1} - {V_0}) \times ({V_2} - {V_0}),d = - \vec n \cdot {V_0}\)
2.    if\(P \cdot \vec n + d \ne 0\), then
3.           return NOT_ON_PLANE;                   //点\(P\)在平面\(\pi \)外;
4.    end if
5.    \({d_0} = \vec n \cdot \left[ {\left( {{V_1} - {V_0}} \right) \times \left( {P - {V_0}} \right)} \right]\);
6.    \({d_1} = \vec n \cdot \left[ {\left( {{V_2} - {V_1}} \right) \times \left( {P - {V_1}} \right)} \right]\);
7.    \({d_2} = \vec n \cdot \left[ {\left( {{V_0} - {V_2}} \right) \times \left( {P - {V_2}} \right)} \right]\);
8.    if \({d_0} \prec 0\)|| \({d_1} \prec 0\)|| \({d_2} \prec 0\), then
9.           return OUSIDE;
10.   else if \({d_0} = = 0\)||\({d_1} = = 0\)|||\({d_2} = = 0\)|, then
11.          return ON_EDGE;
12.   else
13.          return INSIDE;
14.   end if;

4.2.2. 点与三角形的距离

问题描述:在三维空间上,计算点\(P\)到三角形\(\Delta {V_0}{V_1}{V_2}\)的距离,其中,三角形所在的平面是\(\pi \)\({V_i} = ({x_i},{y_i},{z_i}),i \in \{ 0,1,2\} \)

2015-3-23 1-18-27

图4.6 点到三角形的距离

首先阐述一种朴素算法,求出点\(P\)到平面\(\pi \)上的投影点\(P'\),若点\(P'\)在三角形内,则点\(P\)到投影点\(P'\)的距离就是最短的距离;否则,求出点\(P'\)到三角形三条边的距离,取其中距离最小的值,就是点\(P\)到三角形的距离。如图4.6所示,点\({P_1}\)到平面的投影点\({P_1}'\)在三角形内,所以线段\(\overline {{P_1}{P_1}'} \)的长度就是点\({P_1}\)到三角形的距离;点\({P_0}\)到平面的投影点\({P_0}'\)不在三角形内,而线段\(\overline {{P_0}'{P_0}''} \)是投影点\({P_0}'\)到三角形边的最小距离。这种算法主要涉及两个问题:1)点到平面的投影;2)点到线段的距离计算。

这里介绍另外一种效率更高的算法,通过求解方程来求解最小距离,Schneider et al[2]对该算法作了详细地描述,三角形内任意一个点,可以用参数方程的形式表示为:

\[T(s,t) = B + s{\vec e_0} + t{\vec e_1} \tag{4.8}\]

其中,\((s,t) \in D = \left\{ {(s,t):s \in [0,1],t \in [0,1],s + t \le 1} \right\},B = {V_0},{\vec e_0} = {V_1} - {V_0},{\vec e_1} = {V_2} - {V_0}\)

三角形上任意一点与点\(P\)之间的距离是:

\[{\mathop{\rm dist}\nolimits} (P,\Delta ) = \left\| {T(s,t) - P} \right\| \tag{4.9}\]

使用它的平方函数形式,可以表示为:

\[f(s,t) = {\left\| {T(s,t) - P} \right\|^2} = a{s^2} + 2bst + c{t^2} + 2ds + 2et + f \tag{4.10}\]

其中,\(a = {\vec e_0} \cdot {\vec e_0},b = {\vec e_0} \cdot {\vec e_1},c = {\vec e_1} \cdot {\vec e_1},d = {\vec e_0} \cdot (B - P),e = - {\vec e_1} \cdot (B - P),f = (B - P) \cdot (B - P)\)

2015-3-23 1-18-37

图4.7 \(f(s,t)\)的二次函数的三维空间形状

显然,\(a \succ 0,c \succ 0\),函数\(f(s,t)\)是一个关于\(s\)\(t\)的二次函数,它的形状是一个抛物面,如图4.7所示。从上往下看,函数示意图相当于由中心,向外部扩展的一个个椭圆,称每个椭圆是二次函数的一个阶层曲线,同一个椭圆有相同的阶层值,所谓的阶层值就是指\(f(s,t)\)的值,即距离的平方,越靠外的椭圆的阶层值越大,图4.8表示出了抛物面的阶层曲线。该抛物面中心点,就是极值点,极值点处的阶层值最小。

2015-3-23 1-18-50

图4.8 二元二次函数\(f(s,t)\)的阶层曲线

根据微积分的知识,可以对\(f(s,t)\)\(s\)\(t\)的偏导,令偏导值为0,求出的解就是极值点,设极值点是\(({s_c},{t_c})\),如下面的等式所示:

\[\left\{ \begin{array}{l}{f_s}'(s,t) = \frac{{\partial f(s,t)}}{{\partial s}} = 2(as + bt + d) = 0\\{f_t}'(s,t) = \frac{{\partial f(s,t)}}{{\partial t}} = 2(ct + bs + e) = 0\end{array} \right. \tag{4.11}\]

三角形的三个顶点不在同一条直线上,因此,\(ac - {b^2} \succ 0\)一定成立,解等式(4.11),得到极值点\(({s_c},{t_c})\)

\[\left\{ \begin{array}{l}{s_c} = \frac{{be - cd}}{{ac - {b^2}}}\\{t_c} = \frac{{bd - ae}}{{ac - {b^2}}}\end{array} \right. \tag{4.12}\]

2015-3-23 1-19-02

图4.9 \((s,t)\)取值范围把平面划分划分为7个区域

如果求解的\(({s_c},{t_c}) \in D\),那么把它代入(4.10)式求解出来的距离就是最小值,否则,最小值必定出现在三角形的边上。现在分析\(({s_c},{t_c})\)所处的位置,条件\(D\)\(s\)\(t\)的限制,用函数示意图表示,是一个\(s\)轴方向和\(t\)轴方向上\([0,1]\)之间的三角形,如图4.9所示。它把平面划分成了7个区域,若\(({s_c},{t_c}) \in D\),即\(({s_c},{t_c})\)处在区域0中。

2015-3-23 1-19-14

图4.10 极值点\(({s_c},{t_c})\)在区域1,由极值点向外部扩展的阶层曲线

如果\(({s_c},{t_c})\)处在区域1中,函数\(f(s,t)\)是一个抛物面的形状,它的阶层曲线是由极值点中心\(({s_c},{t_c})\)往外部扩展的一个个椭圆,离中心越远的椭圆具有越大的阶层值,因此可以认定最先与直线\(s + t = 1\)接触的椭圆,具有最小的阶层值,如图4.10所示,第一个与直线\(s + t = 1\)接触的椭圆,必定与直线\(s + t = 1\)相切,把\(s = 1 - t\)代入函数\(f(s,t)\)中,得到等式

\[f(1 - t,t) = (a + c - 2b){t^2} + 2(b + e - d - a)t + (a + f) \tag{4.13}\]

\(f(1 - t,t)\)是一个一元二次方程式,因为\(a + c - 2b \succ 0\)\(f(1 - t,t)\)是一个开口向上的抛物线,抛物线的极值点,就是最小值点,对抛物线求导数,得

\[f'(1 - t,t) = 2\left[ {(a + c - 2b)t - (a + d - b - e)} \right] = 0 \tag{4.14}\]

解得

\[t = \frac{{a + d - b - e}}{{a + c - 2b}} \tag{4.16}\]

现在在一维的情况下,分析一元二次方程式\(f(1 - t,t)\)的最小值,\(t\)的取值范围是\([0,1]\),区间\([0,1]\)把整条直线分成三段,分别是\(( - \infty ,0)\)\([0,1]\)\((1, + \infty )\)。若\({t_c} \prec 0\),则\(f(1 - t,t)\)\([0,1]\)上是递增的,当\(t = 0\)时值最小,所以\(f(1,0)\)取得最小值;若\({t_c} \succ 0\),则\(f(1 - t,t)\)\([0,1]\)上是递减的,当\(t = 1\)时值最小,所以\(f(0,1)\)取得最小值;若\(0 \le {t_c} \le 1\),极值点就是最小值点,当\(t = {t_c}\)时值最小,所以\(f(1 - {t_c},{t_c})\)取得最小值点。

同样,也可以把\(t = 1 - s\)代入函数\(f(s,t)\)中,得到一元二次方程式\(f(s,1 - s)\),通过求导数,解得

\[s = \frac{{c + e - b - d}}{{a + c - 2b}} \tag{4.17}\]

接着对一元二次方程式\(f(s,1 - s)\)做分析,也能取得相同的结果。

对这个问题,以另外一种更直观的方式来理解,就是函数\(f(s,t)\)是如图4.8所示的抛物面,若极值点在区域1中,那么最小值点一定落在\(s + t = 1\)的平面上,平面\(s + t = 1\)与抛物面相交的横截面是一条抛物线,相当于把\(s = 1 - t\)代入函数\(f(s,t)\)中,得到抛物线的方程式(4.14),接下来的分析与上面介绍的一样。

区域3和区域5的处理方法与区域1类似,这里不再赘述。

如果\(({s_c},{t_c})\)在区域2内,则函数\(f(s,t)\)与三角形第一个接触的阶层曲线可能边\(s = 0\)或者在\(s + t = 1\)上。\(\partial f/ds = {f_s}'(s,t)\)\(\partial f/dt = {f_t}'(s,t)\)分别表示二次函数\(f(s,t)\)\(s\)\(t\)方向上的偏导数,如果\(({d_x},{d_y}) \cdot \left( {{f_s}'(s,t),{f_t}'(s,t)} \right) \succ 0\),则函数\(f(s,t)\)\((s,t)\)处,沿着\(({d_x},{d_y})\)方向上是递增的,如果\(({d_x},{d_y}) \cdot \left( {{f_s}'(s,t),{f_t}'(s,t)} \right) \prec 0\),则函数\(f(s,t)\)在位置\((s,t)\),沿着\(({d_x},{d_y})\)方向上是递减的。由于极值点落在区域2上,第一个与长方形相切的阶层曲线,只能在\(s = 0\)或者\(s + t = 1\)上,所以函数\(f(s,t)\)在点\((0,1)\)处,方向是\((1,-1)\)和方向\((0,-1)\)上的偏导数不可能都是负的,即\((1, - 1) \cdot \left( {{f_s}'(0,1),{f_t}'(0,1)} \right) = {f_s}'(1,1) - {f_t}'(1,1)\)\((0, - 1) \cdot \left( {{f_s}'(0,1),{f_t}'(0,1)} \right) = - {f_t}'(1,1)\)两个值不可能都是负的。若\({f_s}'(1,1) - {f_t}'(1,1) \prec 0\),则最小值出现在\(s + t = 1\)上;否则,最小值出现在\(s = 0\)上。在确定了最小值出现的直线,就可以把该直线方程,代入函数\(f(s,t)\)中,消去一元,再对一维的一元二次方程式进行分析,该方程式的形状是一个抛物线。例如,确定了最小值出现在直线\(s = 0\)上,代入函数\(f(s,t)\),解得\(f(0,t) = c{t^2} + 2et + f\),该函数的极值点是\(\hat t = - e/c\),若\(\hat t \prec 0\)\(f(0,0)\)取得最小值;若\(\hat t \succ 1\)\(f(0,1)\)取得最小值;若\(0 \le \hat t \le 1\)\(f(0,\hat t)\)取得最小值。

区域4和区域6的处理方法与区域2类似,这里不再赘述。

算法伪码如下所示:

在三维空间上,计算点\(P\)到三角形\(\Delta {V_0}{V_1}{V_2}\)的距离,其中,三角形所在的平面是\(\pi \)\({V_i} = ({x_i},{y_i},{z_i}),i \in \{ 0,1,2\} \))。
1.    \(a = {\vec e_0} \cdot {\vec e_0},b = {\vec e_0} \cdot {\vec e_1},c = {\vec e_1} \cdot {\vec e_1},d = {\vec e_0} \cdot (B - P),e = - {\vec e_1} \cdot (B - P),f = (B - P) \cdot (B - P)\);
2.    \(s = be - cd,t = bd - ae,det = ac - b*b\);
3.    if \(s + t \le det\), then
4.           if\(s \prec 0\), then
5.                  if \(t \prec 0\), then
6.                         //区域4
7.                         if\(d \prec 0\), then
8.                                \(t = 0\);
9.                                if \( - d \ge a\), then
10.                                      \(s = 1\);
11.                               else
12.                                      \(s = - d/a\);
13.                               end if;
14.                        else
15.                               \(s = 0\);
16.                               if\( - e \le 0\), then
17.                                      \(t = 0\);
18.                               else if\( - e \ge c\), then
19.                                      \(t = 1\);
20.                               else
21.                                      \(t = - e/c\);
22.                               end if;
23.                        end if ;
24.                 else
25.                        //区域3
26.                        \(s = 0\);
27.                        if\( - e \le 0\), then
28.                               \(t = 0\);
29.                        else if\( - e \ge c\), then
30.                               \(t = 1\);
31.                        else
32.                               \(t = - e/c\);
33.                        end if;
34.                 end if;
35.          else
36.                 if \(t \prec 0\), then
37.                        //区域5
38.                        \(t = 0\);
39.                        if\( - d \le 0\), then
40.                               \(s = 0\);
41.                        else if \( - d \ge a\), then
42.                               \(s = 1\);
43.                        else
44.                               \(s = - d/a\)
45.                        end if;
46.                 else
47.                        //区域0
48.                        \(t = t/det\),\(s = s/det\);
49.                 else if;
50.          end if;
51.   else
52.          if \(s \prec 0\), then
53.                 //区域2
54.                 if \((b + d) - (c + e) \prec 0\), then
55.                        \(numer = c + e - b - d,denom = a + c - 2b\);
56.                        if \(numer \ge {\rm{de}}nom\), then
57.                               \(s = 1\);
58.                        else
59.                               \(s = numer/denom\);
60.                        end if;
61.                        \(t = 1 - s\);
62.                 else
63.                        \(t = 1 - s\);
64.                        if \( - e \le 0\), then
65.                               \(t = 0\);
66.                        else if \( - e \ge c\), then
67.                               \(t = 1\);
68.                        else
69.                               \(t = - e/c\);
70.                        end if;
71.                 end if;
72.          else if \(t \prec 0\), then
73.                 //区域6
74.                 if \((b + e) - (a + d) \prec 0\), then
75.                        \(numer = c + e - b - d,denom = a + c - 2b\);
76.                        if \(numer \le 0\), then
77.                               \(s = 0\)
78.                        else if \(numer \ge {\rm{de}}nom\), then
79.                               \(s = 1\);
80.                        else
81.                               \(s = numer/denom\);
82.                        end if;
83.                        \(t = 1 - s\);
84.                 else
85.                        \(t = 0\);
86.                        if \( - d \le 0\), then
87.                               \(s = 0\);
88.                        else if \( - d \ge a\), then
89.                               \(s = 1\)
90.                        else
91.                               \(s = - d/a\)
92.                        end if;
93.                 end if;
94.          else
95.                 //区域1
96.                 \(numer = c + e - b - d,denom = a + c - 2b\);
97.                 if \(numer \le 0\), then
98.                        \(s = 0\);
99.                 else if \(numer \ge {\rm{de}}nom\), then
100.                      \(s = 1\);
101.               else
102.                      \(s = numer/denom\);
103.               end if;
104.               \(t = 1 - s\);
105.        end if;
106. end if;
107. return \(\sqrt {a{s^2} + 2bst + c{t^2} + 2ds + 2et + f} \);

spacer

Leave a reply