浏览量:334

6.5. 点与多边形的关系

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

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

6.5. 点与多边形的关系

Haines Eric[2](1994) 详细地总结各种点与简单多边形关系的判定算法,并测试了不同算法的性能。判定点与一般的简单多边形关系的算法,有交叉法(Crossings Test),角度累积法(Angle Summation),增量角度法(Incremental Angle),三角扇法(Triangle Fan),栅格法(Grid Method),基于像素的测试法(Pixel Based Test)。其中,速度最快的是栅格法,但是它需要预处理操作,并且对内存的需求特别大。

判定点与凸多边形关系,同样有很多种方法,最容易想到的一种方法就是判断点是否在凸多边形每条边的外侧,只要存在一条这样的边,就可知点就在凸多边形外,否则,在凸多边形内。这种方法的复杂度是\(O(n)\),而Schneider et al[5](1985)描述了一种基于二分查找的算法,时间复杂度是\(O(\log n)\),是最快的判定算法。

对于点与简单多边形关系的判定问题,本小节将重点介绍交叉法和增量角度法,两种算法的思路完全不同,但是性能相当。对于点与凸多边形关系的判定问题,本小节将介绍基于二分查找的算法。

6.5.1. 交叉法

问题描述:判定点\(Q\)与简单多边形\({\rm P} = \{ {V_0},{V_1},…,{V_{n – 1}}\} \)的关系,其中\({{V}_{i}}=({{x}_{i}},{{y}_{i}})\),设简单多边形是逆时针顺序。

它是基于射线与边交点个数的奇偶性来判定的,基本思想是:由待判定点\(Q\),沿任意一个方向做一条射线,求出射线与多边形的交点个数,若交点个数是奇数,则点在多边形内,若是偶数,则在多边形外。如图6.30所示,点\({{Q}_{1}}\)在多边形外,由它引出的射线与多边形边有4个交点,交点个数是偶数;点\({{Q}_{2}}\)在多边形内,由它引出的射线与多边形边有1个交点,交点个数是奇数。为了算法的简便,射线的方向通常取(1,0)或者是(0,1)。

2015-3-23 18-13-57

图6.30 根据射线和多边形边交点个数的奇偶性,判定点与多边形的关系

由待检测点引出的射线,可能经过多边形的顶点,或者与多边形的边发生重合。如图6.30所示,若\({{Q}_{2}}\)引出的射线经过顶点\({{V}_{3}}\),则计算交点时,射线与边\(\overline{{{V}_{2}}{{V}_{3}}}\)和边\(\overline{{{V}_{3}}{{V}_{4}}}\)都相交,交点个数是2,由此得出的结论“点\({{Q}_{2}}\)在多边形外”是错误的。O’Rourke [4](1998)介绍了一种解决方案,以待测试点为原点,方向是(1,0)作一条射线,求它与多边形的交点个数。对相交点的计数作如下限制:(1)若多边形一条边的一个端点严格位于射线之上,另一个端点严格位于射线之下,相交点个数加1;(2)与射线重合的边,直接忽略不计;(3)一条边位于射线之上,它的一个端点在射线上,相交点个数不变;(4)一条边位于射线之下,它的一个端点在射线上,相交点个数加1。

如图6.31所示,初始情况下,相交点计数值\(n=0\),经过点\({{V}_{15}}\),边\(\overline{{{V}_{15}}{{V}_{16}}}\)在射线上,计数值不变;边\(\overline{{{V}_{14}}{{V}_{15}}}\)与射线重合,忽略;经过点\({{V}_{14}}\),边\(\overline{{{V}_{13}}{{V}_{14}}}\)在射线下方,计数值加1,得\(n=1\);经过点\({{V}_{11}}\),边\(\overline{{{V}_{11}}{{V}_{12}}}\)在射线下方,计数值加1,得\(n=2\);依次类推,最终得到的计数值为\(n=5\),是一个奇数,判定点\(Q\)在多边形内。

设射线的原点是\(Q\)、方向是(1,0),边为\(\overline{{{V}_{i}}{{V}_{i+1}}}\),如何判断射线与边的关系呢?若边\(\overline{{{V}_{i}}{{V}_{i+1}}}\]的两个端点在\(Y\)轴上的分量大于\({{Q}_{y}}\)或者小于\({{Q}_{y}}\),则射线与边\({{V}_{i}}{{V}_{i+1}}\)一定不相交。否则,射线与边可能相交。边可以用参数方程表示为\(L(t)={{V}_{i}}+t({{V}_{i+1}}-{{V}_{i}}),0\le t\le 1\)。射线与边相交,设交点是\(T\],则交点\(Y\)轴上的分量是\({{T}_{y}}={{Q}_{y}}\),则有\({{Q}_{y}}={{y}_{i}}+t({{y}_{i+1}}-{{y}_{i}})\),因此可以解得

\[t = \frac{{{Q_y} – {y_i}}}{{{y_{i + 1}} – {y_i}}} \tag{6.8}\]

把\(t\)代入边的参数方程中,可以得到交点\(X\)轴上的分量

\[{T_x} = \frac{{({Q_y} – {y_{i,y}}) \cdot ({x_{i + 1}} – {x_i}) + {x_i}({y_{i + 1}} – {y_i})}}{{{y_{i + 1}} – {y_i}}} \tag{6.9}\]

设\(d={{T}_{x}}-{{Q}_{x}}\),则有

\[d = \frac{{({Q_y} – {y_i}) \cdot ({x_{i + 1}} – {x_i}) + ({x_i} – {Q_x}) \cdot ({y_{i + 1}} – {y_i})}}{{{y_{i + 1}} – {y_i}}} \tag{6.10}\]

现在只需要判断\(d\)的正负即可,若\(d\prec 0\),射线与边不相交,否则,相交。这里使用了除法,可能出现除数为0的情况,实际上我们并不需要计算出\(d\)的值,因此\(d\)可以改写成如下形式,避免了除法运算可能会发生的错误。可以提前计算出\(({{y}_{i+1}}-{{y}_{i}})\)的值,保存起来,这样能减少减法操作的次数。

\[d = \left[ {({Q_y} – {y_i}) \cdot ({x_{i + 1}} – {x_i}) + ({x_i} – {Q_x}) \cdot ({y_{i + 1}} – {y_i})} \right] \cdot ({y_{i + 1}} – {y_i}) \tag{6.11}\]

2015-3-23 18-14-16

图6.31 射线与多边形相交,多边形上存在点和边与射线重合的情况

判定点与多边形关系的伪码如下所示:

判定\(Q\)与多边形\({\rm P} = \{ {V_0},{V_1},…,{V_{n – 1}}\} \]的关系,其中\({{V}_{i}}=({{x}_{i}},{{y}_{i}})\),设简单多边形是逆时针顺序
return ({OUTSIDE, INSIDE, ON_EDGE})
1.    count = 0;
2.    for each edge \(\overline{{{V}_{i}}{{V}_{i+1}}}\)
3.           if \(Q=={{V}_{i}}\), then
4.                  return ON_EDGE;                                           //点Q与点\({{V}_{i}}\)重合
5.           else if \({{Q}_{y}}=={{y}_{i}}\) && \({{Q}_{y}}=={{y}_{i+1}}\), then
6.                  if \({{x}_{i}}\le {{Q}_{x}}\le {{x}_{i+1}}\) || \({{x}_{i+1}}\le {{Q}_{x}}\le {{x}_{i}}\), then
7.                         return ON_EDGE;                                    //点Q在边\(\overline{{{V}_{i}}{{V}_{i+1}}}\)上
8.                  end if;
9.           else if \({{y}_{i}}\le {{Q}_{y}}<{{y}_{i+1}}\) || \({{y}_{i+1}}\le {{Q}_{y}}<{{y}_{i}}\), then
10.                 \(d=\left[ ({{Q}_{y}}-{{y}_{i}})\cdot ({{x}_{i+1}}-{{x}_{i}})+({{V}_{i,x}}-{{Q}_{x}})\cdot ({{y}_{i+1}}-{{y}_{i}}) \right]\cdot ({{y}_{i+1}}-{{y}_{i}})\);
11.                 if d > 0, then
12.                        count = count + 1;
13.                 else if d == 0, then
14.                        return ON_EDGE;                                    //点Q在边\(\overline{{{V}_{i}}{{V}_{i+1}}}\)上
15.                 end if;
16.          end if;
17.   end for;
18.   return (count % 2 = 0) ? OUTSIDE : INSIDE;

6.5.2. 角度累积法

问题描述:判定点\(Q\)与简单多边形\({\rm P} = \{ {V_0},{V_1},…,{V_{n – 1}}\} \)的关系,其中\({V_i} = ({x_i},{y_i})\),设简单多边形是逆时针顺序。

2015-3-23 18-14-29

图6.32 角度累积法判定点与多边形的关系,(a)点\({{Q}_{1}}\)在多边形外,(b)点\({{Q}_{2}}\)在多边形内

另一种完全不同的一种算法思路,就是角度累积法,假设你站在待检测点\(Q\),待检测点\(Q\)与每条边的两个端点会有一个有向角度,沿着逆时针方向遍历多边形的每一条边,累积有向角度。若\(Q \in {\rm P}\),则最终的角度和为360,相当于绕了一圈;若\(Q \notin {\rm P}\),则最终的角度和为0。如图6.32(a)所示,点\({{Q}_{1}}\)在多边形外,边\(\overline{{{V}_{0}}{{V}_{1}}}\)与点\({{Q}_{1}}\)形成的交角,是一个正数,而边\(\overline{{{V}_{5}}{{V}_{6}}}\)与点\({{Q}_{1}}\)形成的交角是一个负数,沿着逆时针遍历多边形的每一条边,最终得到的有向角度和是0。如图(b)所示,点\({{Q}_{2}}\)在多边形内,逆时针遍历边,与所有边形成的有向角度和是360。计算夹角的方法,可以通过点积:\({{\vec{v}}_{0}}\cdot {{\vec{v}}_{1}}=\left\| {{{\vec{v}}}_{0}} \right\|\cdot \left\| {{{\vec{v}}}_{1}} \right\|\cdot \cos \theta \),或者通过叉积:\({{\vec{v}}_{0}}\times {{\vec{v}}_{1}}=\left\| {{{\vec{v}}}_{0}} \right\|\cdot \left\| {{{\vec{v}}}_{1}} \right\|\cdot \sin \theta \),除以两个向量的长度,得到\(\cos \theta \)或者\(\sin \theta \),再使用反三角函数,计算出角度\(\theta \)。显然,大量的除法、求根、反三角形函数的操作,不仅会降低算法的效率,而且也会降低浮点数的精度,所以这种算法效率和稳定性相对较差。但是,Weiler[6](1994)提出了增量角度法,对角度累计法做了进一步的改进,统计相邻两个顶点与待检测点之间的增量角度,避免了对三角函数的计算,大大提高了算法的效率和稳定性,而且该算法也可以应用在带洞的多边形和复杂多边形。接下来将详细的介绍该算法,这里只考虑非凸的、简单的、不带洞的多边形。

以待检测点为原点,水平的\(X\)轴和垂直的\(Y\)轴,可以把平面分成四个象限,分象限的方法可以采用如下伪码的方法:

int QUADRANT(Q, V)
/* 计算多边形上的顶点V相对于待检测点Q的象限 */
1.    if Q.x == V.x && Q.y == V.y, then                   //说明待检测点Q与顶点Q重合
2.           return 5;
3.    else if V.x > Q.x, then
4.           if V.y > Q.y, then
5.                  return 0;
6.           else
7.                  return 3;
8.           end if;
9.    else
10.          if V.y > Q.y, then
11.                 return 1;
12.          else
13.                 return 2;
14           end if;
15.   end if;

上面所述的分象限法是基于待检测点的位置,会分成如图6.33所示的四个象限,由于上面采用的是大于操作符,所以实际的象限边界如图中的实线所示,把穿过待检测点的虚线向右上位移很小的值。

2015-3-23 18-14-46

图6.33 四个象限

对于待检测点在边上的情况,在后面讨论,这里假设待检测点不在多边形的任意一条边上。对于边上的两个顶点\({{V}_{0}}\)和\({{V}_{1}}\),若相对于待检测点\(Q\),顶点\({{V}_{0}}\)在象限0,顶点\({{V}_{1}}\)在象限1,则增量角就是指顶点\({{V}_{0}}\)到顶点\({{V}_{1}}\)的象限的变化值,即1-0=1。遍历所有的边,累加增量角,若最终的增量角的和等于4或者-4,则说明待检测点\(Q\)在多边形内,否则,不在多边形内,这就是增量算法的基本思想。

2015-3-23 18-15-01

图6.34 存在边在对角象限上

注意,采用简单的减法计算增量角可能会产生错误,需要对增量角进行调整,主要分为两种情况:

  1. 两邻的两个顶点分别在象限0和象限3中。若顶点由象限3跳到象限0,采用简单的减法操作得到的增量角是-3,但是实际上增量角是1,同样对于由象限0跳到象限3时计算出的增量角也需要调整做类似的调整。
  2. 相邻的两个顶点在对角的两个象限上。如图1.34所示的两种情况,图(a)和图(b)中两个三角形所在的象限完全一样,但是图(a)的待检测点在三角形外,图(b)的待检测点在三角形内。因此,这里需要增加一个对增量角的调整,若对于一条边的两个顶点所在的象限采用减法操作,计算出的增量角是\(\Delta =2\)或者\(\Delta =-2\),就需要判断待检测点为原点、方向是(1,0)的射线与该边相交点的X分量,即图1.34中相交点\(T\)的X分量。若\({{T}_{x}}\succ {{Q}_{x}}\),则对增量角\(\Delta \)值取反,否则,增量角\(\Delta \)不变。\({{T}_{x}}\)和\({{Q}_{x}}\)的大小判断方法,可以参考交叉法中的描述。

对采用减法操作得到的增量角进行调整的示意图,如图6.35所示,对上述的两种情况的增量角做出了调整。

2015-3-23 18-15-16

图6.35 对增量角的调整

接下来讨论顶点在边上的情况,算法需要遍历多边形的每个顶点,求出顶点所在的象限,若顶点与待检测点重合的话,说明顶点在边上,算法结束。除此之外,当\(\left\| \Delta \right\|=2\)时,需要计算射线与边的交点\(T\),此时增加一个条件判断,若\({{T}_{x}}={{Q}_{x}}\),则顶点在边上;当\(\left\| \Delta \right\|=1\)时,需要增加待检测点与边关系的判定,若待检测点在边上,则边一定是水平方向或者垂直方向。

对增量角调整的算法伪码如下所示:

int ADJUST_DELTA(Q, v0, v1, delta)
/* Q表示待检测点,v0,v1是简单多边形上连续的两个顶点,delta表示顶点v0到顶点v1的增量角*/
1.    if delta == 3, then
2.           delta = -1;
3.    else if delta == -3, then
4.           delta = 1;
5.    else if (delta == 2 || delta == -2), then
6.           dy = v1.y – v0.y;
7.           d = (Q.y – v0.y) * (v1.x – v0.x) + (v0.x – Q.x) * dy) * dy;
8.           if d > 0, then
9.                  delta = –delta;
10.          else if d == 0, then
11.                 delta = 5;
12.          end if;
13.   else if (delta == 1 || delta == -1), then
14.          if 点Q在线段v0,v1上, then
15.                 delta = 5;
16.          end if;
17.   end if;
18.   return delta;

点与多边形关系判定的增量角算法的伪码如下所示:

判定点Q与简单多边形\({\rm P} = \{ {V_0},{V_1},…,{V_{n – 1}}\} \)的关系,,其中\({V_i} = ({x_i},{y_i})\),设简单多边形是逆时针顺序
return ({OUTSIDE, INSIDE, ON_EDGE})
1.    last = n – 1, current = 0;
2.    sum = 0;                                                   //增量角的和
3.    quadrant0 = QUADRANT (Q, \({{V}_{last}}\));
4.    if quadrant0 == 5, then
5.           return ON_EDGE;
6.    end if;
7.    while( current < n )
8.           quadrant1 = QUADRANT (Q, \({{V}_{current}}\));
9.           if quadrant1 == 5, then
10.                 return ON_EDGE;                      //点Q与点\({{V}_{current}}\)重合
11.          end if;
12.          delta = quadrant1quadrant0;
13.          delta = ADJUST_DELTA(Q,\({{V}_{last}}\),\({{V}_{current}}\),delta);
14.          if delta == 5, then
15.                 return ON_EDGE;                      //点Q在边\(\overline{{{V}_{last}}{{V}_{current}}}\)上
16.          end if;
17.          sum += delta;
18.          last = current;
19.          current += 1;
20.          quadrant0 = quadrant;
21.   end while;
22.   return (sum == 4 || sum == -4 ) ? INSIDE:OUTSIDE;

6.5.3. 二分法

问题描述:判定点\(Q\)与凸多边形\({\rm P} = \{ {V_0},{V_1},…,{V_{n – 1}}\} \)的关系。

对于点与凸多边形关系的判定,Schneider[5]等人(2002)描述了一种复杂度是\(O(\log n)\)的算法,采用了二分查找的思想。

设凸多边形是逆时针排列的,多边形的边可以表示为\({{\vec{e}}_{i}}={{V}_{i+1}}-{{V}_{i}}\),该边上指向外面的法向量可以表示为\({{\vec{n}}_{i}}=({{\vec{e}}_{i,y}},-{{\vec{e}}_{i,x}})\)。若对于待检测点\(Q\),有\({{\vec{n}}_{i}}\cdot (Q-{{V}_{i}})\succ 0\),则它在多边形外;若对于所有的边,都有\({{\vec{n}}_{i}}\cdot (Q-{{V}_{i}})\le 0\),且至少存在一条边使得等号成立,则它在多边形的边上;否则,对于所有的边,都有\({{\vec{n}}_{i}}\cdot (Q-{{V}_{i}})\prec 0\)成立,它在多边形内。判定点\(Q\)是否在多边形内,需要检测所有的边,因此算法的复杂度是\(O(n)\)。

2015-3-23 18-15-28

图6.36 沿着顶\({{V}_{i}}\),\({{V}_{j}}\)的连线分割凸多边形

如图6.36所示的逆时针方向的凸多边形,可以发现连接凸多边形上任意两个顶点\({{V}_{i}}\]和\({{V}_{j}}\),其中\(i < j\),求出边\(\vec{e}=\overline{{{V}_{i}}{{V}_{j}}}\)指向外面的法向量\(\vec{n}=({{\vec{e}}_{y}},-{{\vec{e}}_{x}})\),若\({{\vec{n}}_{i}}\cdot (Q-{{V}_{i}})\succ 0\),说明待检测点Q可能在\({{V}_{i}},{{V}_{i+1}},…,{{V}_{j}}\)形成的子凸多边形内,即灰色区域;否则,可能在凸多边形的另一半子凸多边形上。由凸多边形的这种特性,可以用二分搜索的方法进行点与多边形关系的判定,对于凸多边形\({\rm P} = \{ {V_0},{V_1},…,{V_{n – 1}}\} \),先检测选取顶点\({{V}_{0}}\)和\({{V}_{n/2}}\)作一条边,求出它的法向量,使用上面描述的方法,确定点可能位于哪个子凸多边形中。再在子凸多边形做进一步分割,直到算法的结束。这种二分搜索的方法,可以使用算法的复杂度降低到\(O(\log n)\)。

采用二分搜索的方法判定点与凸多边形关系的伪码如下所示:

      判定点Q与凸多边形\({\rm P} = \{ {V_0},{V_1},…,{V_{n – 1}}\} \)的关系
return ({OUTSIDE, INSIDE, ON_EDGE})
1.    if 凸多边形P是逆时针方向排列, then
2.           isCCW = TRUE;
3.    else
4.           isCCW = FALSE;
5.    end if;
6.    left = right = 0;
7.    while( TRUE )
8.           if (rightleft + n) % n == 1, then
9.                  edge = Vright –Vleft;
10.                 if isCCW, then
11.                        normal = (edge.y , –edge.x);
12.                 else
13.                        normal = (-edge.y, edge.x);
14.                 end if;
15.                 d = normal∙(QVleft);
16.                 if d > 0, then
17.                        return OUTSIDE;
18.                 else if d < 0, then
19.                        return INSIDE;
20.                 else if 点Q在边Vleft,Vright上, then
21.                        return INSIDE;
22.                 else
23.                        return OUTSIDE;
24.                 end if;
25.          end if;
26.          if left < right, then
27.                 mid = (right + left) / 2;
28.          else
29.                 mid = ((right + left + n) / 2) % 2
30.          end if;
31.          edge = VmidVleft;
32.          if isCCW, then
33.                 normal = (edge.y , –edge.x);
34.          else
35.                 normal = (-edge.y, edge.x);
36.          end if;
37.          d = normal∙(Q Vleft);
38.          if d ≥ 0, then
39.                 right = middle;
40.          else
41.                 left = middle;
42.          end if;
43.   end while;

算法实现时,除以2的操作,可以使用移位操作来替换。判断多边形的顶点是逆时针还是顺时针,可以采用求点积和的方法,但是凸多边形的判断会更加简单些。对于连续的三个顶点\({{V}_{i}},{{V}_{i+1}},{{V}_{i+2}}\),令\({{\vec{e}}_{0}}={{V}_{i+1}}-{{V}_{i}},{{\vec{e}}_{1}}={{V}_{i+2}}-{{V}_{i}},\vec{n}=({{\vec{e}}_{0,y}},-{{\vec{e}}_{0,x}})\),若\(\vec{n}\cdot {{\vec{e}}_{1}}\prec 0\),则凸多边形的顶点是逆时针方向排列;否则,是顺时针方向排列。

2015-3-23 18-15-40

图6.37 采用二分搜索法判定点与凸多边形的关系

简单的模拟下算法的过程,如图6.37所示,设凸多边有7个顶点,按照逆时针方向排列,现在对点\({{Q}_{0}}\)进行检测。第1次,left = 0,right = 0,因为left = right,则middle = 3,连接顶点\({{V}_{0}}\)和\({{V}_{3}}\),它的法向量是\({{\vec{n}}_{0}}\),点积\({{\vec{n}}_{0}}\cdot ({{Q}_{0}}-{{V}_{0}})\prec 0\),因此有left = middle = 3;第2次,left = 3,right=0,因为left > right,则middle = 5,连接顶点\({{V}_{3}}\)和\({{V}_{5}}\),它的法向量是\({{\vec{n}}_{1}}\),点积\({{\vec{n}}_{1}}\cdot ({{Q}_{0}}-{{V}_{3}})\succ 0\),则有right = middle = 5;第3次,left = 3,right = 5,因为left < right,则middle = 4,连接顶点\({{V}_{3}}\)和\({{V}_{4}}\),它的法向量是\({{\vec{n}}_{3}}\),点积\({{\vec{n}}_{3}}\cdot ({{Q}_{0}}-{{V}_{3}})\prec 0\),则有left = middle = 4;第4次,left = 4,right = 5,因为(right left + n) % n = 1,进入最后的判定阶段,连接顶点\({{V}_{4}}\)和\({{V}_{5}}\),它的法向量是\({{\vec{n}}_{4}}\),易知\({{\vec{n}}_{3}}\cdot ({{Q}_{0}}-{{V}_{4}})\prec 0\),所以点\({{Q}_{0}}\)在凸多边形内。对于\({{Q}_{1}}\)的算法过程也是相似,读者可以自行模拟其算法过程。

 

spacer

Leave a reply