7.6 凸多边形交集
(Intersection of Convex Polygons)

两个顶点数分别为 $n$ 和 $m$ 的任意多边形的交集复杂度可能是二次的，即 $\Omega(nm)$。例如图 7.12 中多边形的交集包含了 25 个正方形。但两个凸多边形的交集复杂度仅仅是线性的，即 $O(n + m)$。凸多边形的交集是许多算法的关键组成部分，例如判定两组点是否可以被一条直线分割，以及求解二元线性规划问题(Shamos 1978)。Shamos(1978) 提出了第一个线性算法，此后，人们开发了多种不同的算法，均实现了 $O(n + m)$ 的时间复杂度。本节介绍的算法是我与三名本科生共同开发的，这是他们的一项家庭作业所提出的解决方案的综合(O'Rourke, Chien, Olson & Naddor 1982)。我认为这是目前可用的最简单的算法，但这并非是一个客观的观点。

该算法的基本思想很简单，但将其转化为代码却有些棘手。假设两个多边形 $P$ 和 $Q$ 的边界方向与往常一样是逆时针方向，并设 $A$ 和 $B$ 分别为它们各自的有向边。该算法让 $A$ 和 $B$ 互相"追逐(chasing)"，并调整它们的移动速度，使它们在 $\partial P$ 和 $\partial Q$ 的每个交叉点处相遇。基本结构如 Code7.2 所示。图 7.13 展示了该算法运行过程。图中，边 $A$ 和 $B$ 以向量的形式表示。关键点显然在于 $A$ 和 $B$ 的移动规则，我们接下来将对此进行探讨。

设 $a$ 为 $A$ 的表头的索引，$b$ 为 $B$ 的头部的索引。如果 $B$ 指向包含 $A$ 的直线，但并未与其相交，如图 7.14 中所有实线向量所示，那么我们需要推进 $B$，使其接近与 $A$ 的可能交点。这就是推进规则的核心所在。图中的情形可以概括如下: 设 $H(A)$ 为 $A$ 左侧的开半平面。我将使用符号 $A \times B > 0$ 来表示叉积向量的 $z$ 轴坐标 $z > 0$，这意味着从 $A$ 到 $B$ 的最短转向是逆时针方向:

$$ \begin{aligned} & \text{if } A \times B > 0 \text{ and } b \notin H(A), \text{ or} \\ & \text{if } A \times B < 0 \text{ and } b \in H(A), \text{then advance } B. \end{aligned} $$

我们暂时忽略 $P[a]$ 与 $B$ 或 $P[b]$ 与 $A$ 的共线情况。因为 $B \times A = -A \times B$，交换 $A$ 和 $B$ 的角色也适用类似的规则:

$$ \begin{aligned} & \text{if } A \times B < 0 \text{ and } a \notin H(B), \text{ or} \\ & \text{if } A \times B > 0 \text{ and } a \in H(B), \text{then advance } A. \end{aligned} $$

如果两个向量都指向对方，则可以推进其中任何一个。当 $A$ 和 $B$ 都不指向对方时，我们推进位于对方半平面之外的那个。如果两个都在对方半平面之外则推进任意一个。需要一些思考才能意识到，如果 $a \in H(B)$ 且 $b \in H(A)$，那么其中一个向量必然指向另一个。因此上述规则涵盖了所有情况。这些情况可以整理成下表:

$A \times B$	$a \in H(B)$	$b \in H(A)$	Advance Rule
$> 0$	T	T	A
$> 0$	T	F	A or B
$> 0$	F	T	A
$> 0$	F	F	B
$< 0$	T	T	B
$< 0$	T	F	B
$< 0$	F	T	A or B
$< 0$	F	F	A

这些规则通过缩减为下表，它利用了那些任意条目的自由度。这将是我们在下文中采用的推进规则。

$A \times B$	Halfplane Condition	Advance Rule
$> 0$	$b \in H(A)$	A
$> 0$	$b \notin H(A)$	B
$< 0$	$a \in H(B)$	B
$< 0$	$a \notin H(B)$	A

7.6.1. 实现 (Implemetation)

该实现的核心是一个 do-while 循环，它实现了上表中的推进规则。虽然该表的转换很简单直接，但我还没有找到一种简洁的方法，来处理所有外围问题，这些问题有可能是算法的核心功能不堪重负。我们将先讨论一般情况，然后再讨论特殊情况。

多边形顶点存储在两个数组 $P$ 和 $Q$ 中，分别以 $a$ 和 $b$ 索引，顶点数分别为 $n$ 和 $m$。决策所依据的三个主要几何变量均使用 AreaSign计算: $A \times B$ 由 AreaSign(O, A, B) 提供，其中 $O$ 是原点。 $a \in H(B)$ 对应 AreaSign(Q[b1], Q[b], P[a])，其中 $b_1 = (b - 1) \mod n$。$b \in H(A)$ 是类似的表达式。 Code7.17 展示了局部变量、初始化和整体结构。核心的 do-while 循环，如 Code7.18 所示，仅适用于一般情况。特殊情况将在后面的 Code7.21 中讨论。

展开 Code 7.17. Convexlntersect, part (a): top-level structure.>>>

        /* P has n vertices, Q has m vertices. */
        void ConvexIntersect( tPolygoni P, tPolygoni Q, int n, int m )
        {
            int     a, b;           /* indices on P and Q (resp.) */
            int     a1, b1;         /* a-1, b-1 (resp.) */
            tPointi A, B;           /* directed edges on P and Q (resp.) */
            int     cross;          /* sign of z-component of A x B */
            int     bHA, aHB;       /* b in H(A); a in H(b). */
            tPointi Origin = {0,0}; /* (0,0) */
            tPointd p;              /* double point of intersection */
            tPointd q;              /* second point of intersection (for 'e') */
            tInFlag inflag;         /* {Pin, Qin, Unknown}: which inside */
            int     aa, ba;         /* # advs. on a & b indices (after 1st inter.) */
            bool    FirstPoint;     /* Is this first point? (used to initialize).*/
            tPointd p0;             /* The first point. */
            int     code;           /* SegSegInt return code. */
        
            /* Initialize variables. */
            a = 0; b = 0; aa = 0; ba = 0;
            inflag = Unknown; FirstPoint = TRUE;
        
            do {
        
                /* BODY of do-while: see part (b) below. */
        
            /* Quit when both adv. indices have cycled, or one has cycled twice. */
            } while ( ((aa < n) || (ba < m)) && (aa < 2*n) && (ba < 2*m) );
        
            if ( !FirstPoint )      /* If at least one point output, close up. */
                LineTo( p0 );
        
            /* Deal with remaining special cases: not implemented. */
            if ( inflag == Unknown)
                printf("The boundaries of P and Q do not cross.\n");
        }

在 Code7.17 中，变量 inflag 是一个枚举类型，用于跟踪当前哪个多边形在内部。它可能有三个取值 {Pin, Qin, Unknown}。在第一次相交之前，其值为 Unknown。在检测到 $A$ 与 $B$ 的交叉点后，inflag 会记住在刚刚越过交叉点之后，哪个多边形处于局部内部。该标志用于在推进边向量时输出相应的多边形顶点。如果 inflag 在计数器循环过程中始终保持 Unknown，则表明 $\partial P$ 和 $\partial Q$ 没有正常相交。它们要么不相交，要么仅在一点相交，要么一个包含另一个。

展开 Code 7.18. Convexlntersect, part (b): do-while loop. MoveTo_d prints a Postscript "moveto" command for double coordinates.>>>

        do {
            /* Computations of key variables. */
            a1 = (a + n - 1) % n;
            b1 = (b + m - 1) % m;

            SubVec( P[a], P[a1], A );
            SubVec( Q[b], Q[b1], B );

            cross = AreaSign( Origin, A, B );
            aHB   = AreaSign( Q[b1], Q[b], P[a] );
            bHA   = AreaSign( P[a1], P[a], Q[b] );

            /* If A & B intersect, update inflag. */
            code = SegSegInt( P[a1], P[a], Q[b1], Q[b], p );
            if ( code == 'i' || code == 'v' ) {
                if ( inflag == Unknown && FirstPoint ) {
                    aa = ba = 0;
                    FirstPoint = FALSE;
                    p0[X] = p[X]; p0[Y] = p[Y];
                    MoveTo_d( p0 );
                }
                inflag = InOut( p, inflag, aHB, bHA );
            }

            /*---Advance rules-----*/
            /* SPECIAL CASES: see part (c) below */
            /* Generic cases. */
            else if ( cross >= 0 ) {
                if ( bHA > 0)
                    a = Advance( a, &aa, n, inflag == Pin, P[a] );
                else
                    b = Advance( b, &ba, m, inflag == Qin, Q[b] );
            }
            else /* if (cross < 0) */ {
                if ( aHB > 0)
                    b = Advance( b, &ba, m, inflag == Qin, Q[b] );
                else
                    a = Advance( a, &aa, n, inflag == Pin, P[a] );
            }

        } while ( ((aa < n) || (ba < m)) && (aa < 2*n) && (ba < 2*m) );

循环的终止条件在概念上很简单，但在实践中很棘手。可以等待第一个输出点的返回(练习 7.6.2[1])，但此处所示的版本基于边向量索引来终止。当 $a$ 和 $b$ 都绕各自的多边形循环一圈后，结束。在某些退化相交的情况下，其中一个索引不会循环。因此，当任一索引在第一次相交后循环两次时，while 语句也会终止，此时计数器 aa 和 ba (后缀 a 代表 "advances") 将被重置。

除了初始化细节外，循环(Code7.18)内的基本操作包括: 使用 SegSegInt 计算边向量 $A$ 和 $B$ 的交点，如果确实相交，则调用函数 InOut 打印交点 $p$ 并切换 inflag 标志。最后调用 Advance 根据规则推进 $a$ 或 $b$，可能还会打印一个顶点。当 SegSegInt 返回 '1' 或 'v' 时，即认为产生了交点。返回 'e' 表示共线重叠，除了稍后讨论的特殊情况外，不会切换标志，也不会产生任何输出。函数 InOut(Code7.19)打印交点，然后根据哪个边向量的头部位于另一个向量的半平面内来决定如何设置 inflag。如果情况无法确定，则不切换标志。

展开 Code 7.19. InOut. LineTo_d prints a Postscript "lineto" command..>>>

        tInFlag InOut( tPointd tp, tInFlag inflag, int aHB, int bHA )
        {
            LineTo_d( p );
        
            /* Update inflag. */
            if      ( aHB > 0 )
                return Pin;
            else if ( bHA > 0 )
                return Qin;
            else    /* Keep status quo. */
                return inflag;
        }

函数 Advance(Code7.20) 会递增计数器。如果刚经过的顶点在内部inside，则将其打印出来。注意，当 inflag` 被设置为例如 Pin 时，我们无法确定 $a$ 是否真的在内部。只知道除了一些特殊情况，在刚过交点之后的位置，$A$ 是在内部的。但是当 Advance 递增 $a$ 时，就可以确定 $a$ 确实位于交集之中。

展开 Code 7.20. Advance. LineTo_i prints a "lineto" command for ints.>>>

        int Advance( int a, int *aa, int n, bool inside, tPointi v )
        {
            if ( inside )
                LineTo_i( v );
        
            (*aa)++;
            return (a+1) % n;
        }

在讨论特殊情况之前，我们先展示一个相对通用的示例上的输出，如图 7.15 所示。代码为这些多边形生成以下 Postscript 输出：

                     5.00    8.00   moveto
                     5.00    8.00   lineto
                    13.00    0.00   lineto
                    19       2      lineto
                    24      10      lineto
                    24.00   21.33   lineto
                    17.80   29.60   lineto
                    14.67   31.17   lineto
                     1.62   23.01   lineto
                     0.72   19.12   lineto
                     2.50   12.00   lineto
                     5.00    8.00   lineto

图 7.13 中的"动画"展示了此示例中 $A$ 和 $B$ 的推进过程。该示例并非完全通用，因为 $P$ 的边 $(1, 2)$ 与 $Q$ 的边 $(4, 0)$ 共线重叠。

输出中整数和实数的交替出现，表示该点是顶点(在 Advance 中输出打印)还是计算出的交点(在 InOut 中输出打印)。列表开头和结尾的点 $(5,8)$ 重复出现是 Postscript 初始化和闭合处理的产物，很容易就可以将其去除。下面，我们将看到更棘手的点重复问题。

特殊情况 (Special Cases)

最后，我们来看一些特殊情况，这些情况种类繁多，令人眼花缭乱。图 7.16 展示了一些示例。

所有特殊情况都是由三个特殊的几何关系造成的: $A \times B = 0$，$a$ 位于包含 $B$ 的直线上时，以及 $b$ 位于包含 $A$ 的直线上。这三种情况都通过 AreaSign 返回 0 来表示。一些组合可以通过"通用情况(generic cases)"的推进规则(Code7.18)得到适当处理。但在 Code7.21 中，有三种情况被单独拎出来处理:

展开 Code 7.21. Convexlntersect, part (c): special cases.>>>

        /*....Advance rules (continued from part (a)) */

        /* Special case: A & B overlap and oppositely oriented. */
        if ( ( code == 'e' ) && ( Dot( A, B ) < 0 ) )
            PrintSharedSeg( p, q ), exit(EXIT_SUCCESS);
        
        /* Special case: A & B parallel and separated. */
        if ( (cross == 0) && ( aHB < 0) && ( bHA < 0 ) )
            printf("%%P and Q are disjoint.\n"), exit(EXIT_SUCCESS);
        
        /* Special case: A & B collinear. */
        else if ( (cross == 0) && ( aHB == 0) && ( bHA == 0 ) ) {
            /* Advance but do not output point. */
            if ( inflag == Pin )
                b = Advance( b, &ba, m, inflag == Qin, Q[b] );
            else
                a = Advance( a, &aa, n, inflag == Pin, P[a] );
        }
        
        /* Generic cases (continued in part (b) above)..... */

如果 $A$ 和 $B$ 重叠(返回 'e')且方向相反($A \cdot B < 0$)，那么它们的重叠部分就是总交集。因为 $P,Q$ 必须位于包含 $A,B$ 的直线两侧。稍后将描述如何找到相交线段。
如果 $A$ 和 $B$ 平行($A \times B = 0$)且 $a$ 严格在 $B$ 的右侧，同时 $b$ 严格在 $A$ 的右侧，那么我们可以得出结论 $P \cap Q = \emptyset$。

inside

回顾 7.2 节，即使当返回码 'e' 表示重叠，我们的 SegSegInt 函数也只返回一个交点 $p$。现在我们需要实际的重叠部分，因此修改 SegSegInt 以返回第二个点 $q$，该点由 ParallelInt 计算得出，使得 $pq$ 为相交线段。对 ParallelInt 的修改如 Code7.22 所示，系统地识别出六种可能的重叠情况： $cd \subseteq ab$, $ab \subseteq cd$$, $c \in ab$(两种情况)，以及 $d \in ab$(两种情况)，并适当地设置 $p$ 或 $q$。

展开 Code 7.22. ParallelInt; See Code 7.3 for Assigndi.>>>

        char ParallelInt( tPointi a, tPointi b, tPointi c, tPointi d, tPointd p, tPointd q )
        {
            if ( !Collinear( a, b, c ) )
                return '0';
        
            if ( Between( a, b, c ) && Between( a, b, d ) ) {
                Assigndi( p, c );
                Assigndi( q, d );
                return 'e';
            }
        
            if ( Between( c, d, a ) && Between( c, d, b ) ) {
                Assigndi( p, a );
                Assigndi( q, b );
                return 'e';
            }
        
            if ( Between( a, b, c ) && Between( c, d, b ) ) {
                Assigndi( p, c );
                Assigndi( q, b );
                return 'e';
            }
        
            if ( Between( a, b, c ) && Between( c, d, a ) ) {
                Assigndi( p, c );
                Assigndi( q, a );
                return 'e';
            }
        
            if ( Between( a, b, d ) && Between( c, d, b ) ) {
                Assigndi( p, d );
                Assigndi( q, b );
                return 'e';
            }
        
            if ( Between( a, b, d ) && Between( c, d, a ) ) {
                Assigndi( p, d );
                Assigndi( q, a );
                return 'e';
            }
        
            return '0';
        }

最后，我们列举了所展示代码的诸多不足，并提出了改进建议:

当循环结束时 inflag 仍为 Unknown 时，可能是 $P \subset Q$ 或 $Q \subset P$ 或 $P \cap Q = \emptyset$，甚至可能是 $P \cap Q = v$，其中 $v$ 是一个顶点。后一种情况是因为 InOut 有时会维持现状。处理这些情况需要额外的工作。虽然这些情况本身并不难处理，但如果能自动区分它们会更好。
使用两个计数器来终止循环，显得很笨拙。检查第一个点是否重复出现也是同样的问题。
虽然从图 7.15 的例子中看不出来，但退化可能导致点被多次输出，有时既作为顶点又作为交点输出。例如，两个共享一个角的嵌套正方形的输出就包含了这个序列:
```
                0.00    50.00   lineto
                0       50      lineto
                0.00    50.00   lineto
```
这种重复点可能会干扰另一个期望所有多边形顶点都不相同的程序。虽然抑制重复输出并不难(练习 7.6.2[2])，但这直接引发了整数与浮点数的问题。例如，必须判定整数 50 是否等于浮点数 50.00，而后者是在 SegSegInt 中通过双精度计算得出的。不精确的计算可能导致相同的点被视为不同，使用"模糊fuzz"因子将不可避免地导致某些不同的点被误判为相等。
最大的缺点是需要特殊处理许多情况，通常涉及一些微妙的逻辑。如果通用代码能自然地特化(specialize naturally)会更令人满意。我将其作为一个开放性问题(练习 7.6.2[4])。

7.6.2. 练习

[Programming] 循环终止 (Peter Schorn) Loop termination。修改代码，使循环终止条件为，是否经过了第一个输出顶点，而不依赖于循环计数器 aa 和 ba。在两个三角形 $(3, 4), (6, 4), (4, 7)$ 和 $(4, 7), (2, 5), (6, 2)$ 上测试你的代码。
[Programming] 重复点 Duplicate points。修改代码以抑制重复点的输出。不要关心上面讨论过的浮点数与整数的问题，而直接使用 == 比较 int 和 double 类型的数据，这将强制转换为 double 进行比较。如果你的代码能够运行，试着找一个能破坏它的案例。这可能相当困难，甚至在某些机器上是不可能的。
有理数 Rationals。
1. 讨论(但不要实现)如何将交点表示为精确的有理数，即两个整数的比值。
2. 设计一个布尔函数来判断两个这样的有理数是否相等。例如，$2/6 = 127131/381393$。
3. [Programming] 实现有理数相等的函数。
[Open] 整洁代码 Clean code。设计一组推进规则，要能够比本文的例程更优雅地处理特殊情况。理想情况下，图 7.16 中的所有情况，以及 $P \subset Q$, $Q \subset P$, $P \cap Q = \emptyset $，都能被自然地处理。

\(A \times B\)	\(a \in H(B)\)	\(b \in H(A)\)	Advance Rule
\(> 0\)	T	T	A
\(> 0\)	T	F	A or B
\(> 0\)	F	T	A
\(> 0\)	F	F	B
\(< 0\)	T	T	B
\(< 0\)	T	F	B
\(< 0\)	F	T	A or B
\(< 0\)	F	F	A

\(A \times B\)	Halfplane Condition	Advance Rule
\(> 0\)	\(b \in H(A)\)	A
\(> 0\)	\(b \notin H(A)\)	B
\(< 0\)	\(a \in H(B)\)	B
\(< 0\)	\(a \notin H(B)\)	A

7.6 凸多边形交集(Intersection of Convex Polygons)

7.6.1. 实现 (Implemetation)

特殊情况 (Special Cases)

7.6.2. 练习

7.6 凸多边形交集
(Intersection of Convex Polygons)