7.4 多边形中的点
(Point in Polygon)

每当我们用鼠标在屏幕上点击某个形状时，就要解决这么一个问题：给定一个固定的多边形 $P$ 和一个查询点 $q$，判断 $q \in P$ 是否成立? 尽管一些特定的硬件可能可以避开几何运算，但我们在这里还是从计算几何的角度来考虑这个问题。

如果 $P$ 是凸多边形，最直接的方法就是，对多边形的每条边执行一次 LeftOn 函数(Code1.6)。事实上，我们在3.7 节中的二维增量凸包算法使用的正是这种技术。这可以改进到 $O(\log n)$，我们将其留作练习 7.4.3[1]。

当 $P$ 为非凸多边形时问题就更有趣了。目前已经流行了两种截然不同的方法: 计算射线交叉次数(counting ray crossings)和计算绕卷数(computing winding numbers)。两者的时间复杂度都是 $O(n)$，但其中一种比另一种快得多。这些算法将是接下来两个小节的主题。

7.4.1. 绕卷数 (Winding Number)

我们首先介绍一种在数学上很优雅的方法，但遗憾的是，事实证明它在实践中表现很差。它基于多边形"绕卷数(winding number)"的概念。想象你站在点 $q$ 上。当一个点 $p$ 沿逆时针方向完整遍历 $\partial P$ 时，你一边注视着 $p$，一边随之转动身体，始终保持面向 $p$。如果 $q \in P$，你会转满一整圈，即 $2\pi$ 弧度。如果 $q \notin P$，你的总转角将恰好为零。通常我们约定，逆时针转动为正，顺时针转动为负。如果 $P$ 是凸多边形，这一点很容易看出。即使 $P$ 是任意多边形，直觉上我也觉得它是对的。毕竟，你最终回到了初始的朝向，所以总转角必须是整数圈的旋转。如图 7.6 所示。我们将不再停下来证明这一命题。点 $q$ 关于 $P$ 的绕卷数(winding number)被定义为 $\partial P$ 绕 $q$ 旋转的圈数。即总有符号转角(记为 $\omega$)除以 $2\pi$。我们将计算的细节留给练习 7.4.3[8]。

尽管绕卷数算法很吸引人，但它依赖于浮点运算，特别是三角函数运算。这使得它在一般硬件上比我们接下来要讨论的射线交叉算法慢得多。 Haines 1994 做过实验对比，它慢了二十多倍！这也顺便说明了在大 $O$ 符号表示法中盲目忽略常数项的风险。

7.4.2. 射线交叉法 (Ray Crossing)

从点 $q$ 沿任意方向画一条射线 $R$，假如就是 $+x$ 方向，并计算 $R$ 与 $\partial P$ 的交点数量。交点数量分别为奇数或偶数，表明点 $q$ 在 $P$ 内或外。例如图 7.7 所示，点 $q_2$ 的射线就有两个交点，想象沿 $R$ 从无穷远移动到 $q_2$。第一个交点穿入 $P$ 的内部，第二个交点移到外部。因此 $q_2 \notin P$。类似的推理表明，图中的 $q_1$，必然在 $P$ 内部，其射线有五个交点。

尽管这个想法很简单，但实现起来却很脆，因为必须处理的边界问题太多了。如图 7.7 中的点 $q_3$，射线可能击中顶点或与某条边共线。还有一种可能是 $q$ 直接位于 $\partial P$ 上，在这种情况下，我们希望得出结论 $q \in P$，因为 $P$ 是闭集。注意，即使是在"没有三个多边形顶点共线"的假设下，也无法排除所有这些"退化"情况。

固定 $R$ 水平向右。如下条件可以消除大部分困难: 一条边 $e$ 要算作 $R$ 的一次穿越(crossing)，$e$ 的一个端点必须严格在 $R$ 上方，而另一个端点在 $R$ 上或其下方。非正式地说，$e$ 被认为包含其较低的端点但排除其较高的端点。将这一约束应用于图中的 $q_3$ 上，边 $(1, 2)$ 和 $(2, 3)$ 不算穿越，因为两条边都没有严格在上方的端点； $(6, 7)$ 和 $(7, 8)$ 算作穿越，因为 $v_7$ 在射线上或下方。$(3, 4)$ 和 $(4, 5)$ 不穿越，而 $(5, 6)$、$(10, 11)$ 和 $(11, 12)$ 都穿越。总共五次穿越意味着 $q_3 \in P$。注意，与射线共线的边不被算作穿越，因为它没有任何点严格在其上方。

在揭示这一约束仍未解决的问题之前，我们先用一个简单函数 InPoly0 来实现这一思想，如 Code7.12 所示。首先平移整个多边形，使 $q$ 变为原点，且 $R$ 与 $x$ 轴正半轴重合。这一步是不必要的，但会使代码更清晰。在遍历所有边 $e = (i - 1, i)$ 的循环中，它根据上述定义检查 $e$ 是否跨越 $x$ 轴。如果跨越，则通过解线段 $e$ 与直线 $y = 0$ 的线性方程求得交点的 $x$ 坐标：

展开 Code 7.12 InPoly0. Typedefs in Code 7.4 and 1.2.>>>

        char InPoly0( tPointi q, tPolygoni P, int n )
        {
            int     i, il;          /* point index; il = i-1 mod n */
            int     d;              /* dimension index */
            double  x;              /* x intersection of e with ray */
            int     Rcross = 0;     /* number of right edge/ray crossings */
        
            /* Shift so that q is the origin. Note this destroys the polygon.
               This is done for pedagogical clarity. */
            for( i = 0; i < n; i++ ) {
                for( d = 0; d < DIM; d++ )
                    P[i][d] = P[i][d] - q[d];
            }
        
            /* For each edge e = (i-1, i), see if crosses ray. */
            for( i = 0; i < n; i++ ) {
                il = ( i + n - 1 ) % n;
        
                /* If e "straddles" the x axis... */
                if( ( ( P[i][Y] > 0 ) && ( P[il][Y] <= 0 ) ) ||
                    ( ( P[il][Y] > 0 ) && ( P[i][Y] <= 0 ) ) ) {
        
                    /* ... compute intersection with x axis. */
                    x = (P[i][X] * P[il][Y] - P[il][X] * P[i][Y])
                        / (double)(P[il][Y] - P[i][Y]);
        
                    /* e crosses ray if strictly positive intersection. */
                    if (x > 0) Rcross++;
                }
            } /* end for */
        
            /* q inside iff an odd number of crossings. */
            if( (Rcross % 2) == 1 ) return 'i';
            else return 'o';
        }

$$ \begin{equation}\label{f7.9} y - y_{i-1} = (x - x_{i-1})(y_i - y_{i-1})/(x_i - x_{i-1}) \tag{7.9} \end{equation} $$

令 $y = 0$，这里 $(x_{i-1}, y_{i-1})$ 和 $(x_i, y_i)$ 是 $e$ 的端点。注意代码中 $x$ 是 double 类型。虽然可以消除这种对浮点计算的依赖(练习 7.4.3[7])，但我们将保留它以使注意力集中在算法上。每当交点严格在原点右侧时，就计数一次穿越。函数返回字符 'i' 或 'o' 分别表示"在内(in)"或"在外(out)"。

这段代码存在一个除了浮点计算问题外的缺陷，虽然它对于任何严格位于 $P$ 内部的点都能返回一致的答案，但它不能正确地处理位于 $\partial P$上的点。如果在图 7.7 中将 $q_3$ 水平移动到 $v_4$，InPoly0 将返回 'i'，但如果将 $q_3$ 放置在 $v_5$，它返回 'o'。该代码对于 $\partial P$ 上点的输出很复杂，如图 7.8 所示。让我们分析一下为什么 $v_{27}$ 被认为是在内部。边 $(26, 27)$ 和 $(27, 0)$ 都不算作穿越，这是因为语句 if (x > 0) 中使用了严格不等式。但是，$v_{27}$ 与 $q_3$ 相同射线有五个穿越点，因此 $v_{27} \in P$。请注意，$v_{22}$ 这个表面上相似的顶点却被判定为外部。现在描述这段代码如何处理边就容易一些:左侧边和水平底边算作在内，右侧边和水平顶边算作在外。

稍微有点强迫症的人都会对这种混杂的处理方式不满。但对于某些应用，特别是 GIS(Geographic Information Systems 地理信息系统)，这种或类似的行为是首选的，因为它具有这样一个性质：在一个区域被划分成许多多边形时，每个点将恰好"在"一个多边形内。这并不显而易见，但我们将视其为事实(练习 7.4.3[4])。其他应用则要求对边界点进行一致的处理，或者区分严格内部、严格外部和边界上这三种情况。

虽然可以通过修改 InPoly0 增加特定的测试，来检查查询点是否位于边界的每条边上，但利用上述关于边处理的特征描述，可以用相对简单的代码，完全确定 $q$ 与 $P$ 的关系。其思路如下，如果我们关于原点将 $P$ 反射形成一个新的多边形 $P'$，那么所有的左边和底边将分别变成顶边和右边，反之亦然。因此，边上的点可以通过 InPoly0 对 $P$ 和 $P'$ 中的同一个 $q$ 给出不同的结果。实际上我们不需要执行反射操作，使用向左的射线，并采用偏向下方的跨越性测试，可以达到相同的效果。

此方法处理了位于边相对内部的点，但不适用于顶点。例如，由于 x > 0 测试，$v_{16}$ 的两条射线都有零个交点。通过显式地检查 $q$ 是否为顶点可以处理这种情况。函数 InPoly1 代码 7.13，它返回四个字符之一 ${i, o, e, v}$，代表这些不同的情况：

展开 Code 7.13 InPoly1. (Some portions shared with InPoly0 are deleted above.)>>>

        char InPoly1( tPointi q, tPolygoni P, int n )
        {
            /* Declarations of i, il, d, x same as in InPoly0; DELETED here. */
            int     Rcross = 0;       /* number of right edge/ray crossings */
            int     Lcross = 0;       /* number of left edge/ray crossings */
            bool    Rstrad, Lstrad;   /* flags indicating the edge straddles the x axis. */
        
            /* Shift so that q is the origin, same as in InPoly0; DELETED here. */
        
            /* For each edge e = (i-1, i), see if crosses rays. */
            for( i = 0; i < n; i++ ) {
                /* First check if q = (0, 0) is a vertex. */
                if ( P[i][X] == 0 && P[i][Y] == 0 ) return 'v';
                il = ( i + n - 1 ) % n;
        
                /* Check if e straddles x axis, with bias above/below. */
                Rstrad = ( P[i][Y] > 0 ) != ( P[il][Y] > 0 );
                Lstrad = ( P[i][Y] < 0 ) != ( P[il][Y] < 0 );
        
                if( Rstrad || Lstrad ) {
                    /* Compute intersection of e with x axis. */
                    x = (P[i][X] * P[il][Y] - P[il][X] * P[i][Y])
                        / (double)(P[il][Y] - P[i][Y]);
        
                    if (Rstrad && x > 0) Rcross++;
                    if (Lstrad && x < 0) Lcross++;
                }
            } /* end straddle computation */
            } /* end for */
        
            /* q on an edge if L/R cross counts are not the same parity. */
            if( ( Rcross % 2 ) != ( Lcross % 2 ) ) return 'e';
        
            /* q inside iff an odd number of crossings. */
            if( ( Rcross % 2 ) == 1 ) return 'i';
            else return 'o';
        }

'i': $q$ 严格在内部。
'o': $q$ 严格在外部。
'e': $q$ 在某条边上，但不是端点。
'v': $q$ 是一个顶点。

关于这段代码需要说明几点。InPoly0 中的跨越测试

        ( ( P[i] [Y] > 0 ) && ( P[i1] [Y] <= 0 ) ) || ( ( P[i1] [Y] > 0 ) && ( P[i] [Y] <= 0 ) )

在 InPoly1 中已被替换为将该表达式赋值给布尔变量 Rstrad:

        ( P[i] [Y] > 0 ) != ( P[i1] [Y] > 0 )

虽然不那么明显，但这两个表达式在逻辑上是等价的。当且仅当 $e$ 的一个端点严格在 $x$ 轴上方而另一个不在，即另一个在轴上或下方时，Rstrad 为真(TRUE)。这种更简洁的形式使得更容易看出 Lstrad 的正确定义，只需反转不等式以偏向下方。现在，只要这两个跨越变量中有一个为真，就需要计算 $x$，这仅排除了经过 $q = (0, 0)$ 的边，并且顺便防止了除以 0。最后，何时返回 'e' 的判定，简化为查看两个射线穿越计数，是否具有不同的奇偶性。

7.4.3. 练习

凸多边形内的点 Point in convex polygon。设计一个算法，在 $O(\log n)$ 时间内判断一个点是否在一个凸多边形内。
最坏情况下的射线交叉 Worst ray crossings。对于任意查询点，射线交叉算法能否在 $O(\log n)$ 时间内运行?
[Programming] 除法 vs.乘法 Division vs.multiplication。在某些机器上，浮点除法可能比乘法慢二十倍。修改 InPoly1 代码以避免公式(\ref{f7.9}) 中的除法，并计时，看看在你的机器上是否有显著差异。
多边形的镶嵌 Tessellation by polygons。论证在将平面区域划分为多边形的分区时，InPoly0 最多将一个点 $q$ 划分到一个多边形内。
[Programming] 加速 Speed-up。通过避免在跨越线段位于射线负侧时计算 $x$，来加速 InPoly0(Code7.12)。
[Programming, easy] 避免平移 Avoid translation。开发一个新版本的 InPoly1(Code7.13)，避免对 $P$ 和 $q$ 进行不必要的平移。
[Programming] 整数射线交叉 Integer ray crossing。
1. 修改 InPoly0以规避唯一的浮点计算。使用 AreaSign(Code4.23)。
2. 使用 AreaSign 的结果来判断 $q$ 是否位于边上，从而在不向两个方向发射射线的情况下实现 InPoly1 的功能。
[Programming] 绕卷数 Winding number。实现绕卷数算法。所需的基本例程是计算从点 $q$ 看多边形边 $e_i$ 所张成的角度。角度 $\theta_i$ 可以通过叉积求得，$v_i \times v_{i+1} = |v_i||v_{i+1}| \sin \theta_i$。回想一下，Code 1.5中的 Area2(q, P[i], P[i+1]) 就是这个叉积的模长。然后必须除以向量的长度，计算反正弦值，并将所有 $i$ 对应的角度求和。
编写代码来计算这个角度和。要特别注意库函数 asin 返回的角度范围。所有逆时针转向是正角度，顺时针转向是负角度。决定当 $|v_i| = 0$ 或 $|v_{i+1}| = 0$ 时应该做什么。

7.4 多边形中的点(Point in Polygon)

7.4.1. 绕卷数 (Winding Number)

7.4.2. 射线交叉法 (Ray Crossing)

7.4.3. 练习

7.4 多边形中的点
(Point in Polygon)