7.2 线段之间的交点
(Segment-Segment Intersection)

在1.5节中，我们花了一些时间编写代码来检测两条线段的相交，用于三角分割(Intersect，代码1.9)，但我们并未真正计算过交点。三角分割算法不需要它，而且计算交点会迫使我们离开舒适的整数坐标世界。然而，在许多应用中，我们需要交点的浮点坐标。在7.6节和7.8节中，我们将需要用它来计算两个多边形的交点。尽管存在一些潜在的坑，幸运的是，计算交点并不难，而且必要的浮点运算不像在其他问题下那样麻烦。在本节中，我们将编写代码来完成这项任务。

虽然利用第1章中的布尔交(Boolean Intersect)可以稍微简化计算，但这里选择进行独立计算。设两条线段的端点分别为 $a,b$ 以及 $c,d$，并设 $L_{ab}$ 和 $L_{cd}$ 为包含这两条线段的直线。计算交点的一种常用方法是联立求解 $L_{ab}, L_{cd}$ 的斜截式方程，即包含两个未知数的两个方程，得到交点坐标 $x, y$。但是我们将使用参数来表示两条线段，因为这样的变量含义似乎更直观。我们将在7.3节中看到，参数化方法可以很好地推广到更复杂的相交计算中。

令向量 $A = b - a，C = d - c$表示线段方向。那么直线 $L_{ab}$ 上的任意一点都可以表示为 $p(s) = a + sA$。其几何意义是，从 $L_{ab}$ 上的点 $a$ 出发，将 $A$ 缩放 $s$ 倍，沿直线移动一段距离，如图7.1所示。变量 $s$ 被称为该方程的参数。当 $s = 0, s = 1$, 以及 $s = \frac{1}{2}$ 时，分别可得: $p(0) = a, p(1) = a + A = a + b - a = b$，以及 $p(\frac{1}{2}) = \frac{a + b}{2}$。这些例子表明，$s \in [0, 1]$时，$p(s)$ 代表了线段 $ab$ 上的所有点，其中 $s$ 的值表示端点之间距离的比例。特别地，$s$ 的极值对应于端点。

同样地，我们可以将第二条线段上的点表示为 $q(t) = c + tC$，其中 $t \in [0, 1]$。线段之间的交点，由使得 $p(s)$ 等于 $q(t)$ 的 $s$ 和 $t$ 值指定，$a + sA = c + tC$。这是包含两个未知数的 $x, y$ 方程，两者都以 $s$ 和 $t$ 为未知数。通常我们用下标 0 和 1 分别表示 $x$ 和 $y$ 坐标，有

$$ \begin{equation}\label{f7.1} s = [a_0(d_1 - c_1) + c_0(a_1 - d_1) + d_0(c_1 - a_1)]/D \end{equation} \tag{7.1} $$ $$ \begin{equation}\label{f7.2} t = [a_0(c_1 - b_1) + b_0(a_1 - c_1) + c_0(b_1 - a_1)]/D \end{equation} \tag{7.2} $$ $$ \begin{equation}\label{f7.3} D = a_0(d_1 - c_1) + b_0(c_1 - d_1) + d_0(b_1 - a_1) + c_0(a_1 - b_1) \end{equation} \tag{7.3} $$

在这些方程中，除零是有可能的。当且仅当两条直线平行时，分母 $D$ 才会为零，这一命题留作练习 7.3.2[1]。正如我们在1.5.4节中详细说明的那样，有些平行线段涉及相交，而有些则不相交。我们暂时将平行线段视为不相交。根据上述方程，可以粗略的写出代码 7.1 中所示的代码。我们将首先描述这段代码，然后挑毛病，最后改它。

展开 Code 7.1 Segment-segment intersection code: rough attempt. >>>

        #define X 0
        #define Y 1
        typedef enum {FALSE, TRUE} bool;
        #define DIM 2               /* Dimension of points */
        typedef int tPointi[DIM];   /* Type integer point */
        typedef double tPointd[DIM];/* Type double point */
        
        bool SegSegInt( tPointi a, tPointi b, tPointi c, tPointi d, tPointd p )
        {
            double s, t;           /* Parameters of the parametric eqns. */
            double num, denom;     /* Numerator and denominator of eqns. */
        
            denom = a[X] * ( d[Y] - c[Y] ) +
                    b[X] * ( c[Y] - d[Y] ) +
                    d[X] * ( b[Y] - a[Y] ) +
                    c[X] * ( a[Y] - b[Y] );
        
            /* If denom is zero, then segments are parallel. */
            if (denom == 0.0)
                return FALSE;
        
            num = a[X] * ( d[Y] - c[Y] ) +
                  c[X] * ( a[Y] - d[Y] ) +
                  d[X] * ( c[Y] - a[Y] );
            s = num / denom;
        
            num = -( a[X] * ( c[Y] - b[Y] ) +
                     b[X] * ( a[Y] - c[Y] ) +
                     c[X] * ( b[Y] - a[Y] ));
            t = num / denom;
        
            p[X] = a[X] + s * ( b[X] - a[X] );
            p[Y] = a[Y] + s * ( b[Y] - a[Y] );
        
            if ( (0.0 <= s) && (s <= 1.0) &&
                 (0.0 <= t) && (t <= 1.0) )
                return TRUE;
            else
                return FALSE;
        }

上述代码以四个端点的整数坐标作为输入，有两种输出。它返回一个布尔值，指示线段是否相交，同时在 $p$ 中返回交点的双精度坐标。注意 $p$ 的类型是 double tPointd。分子和分母的计算完全对应式 ($\ref{f7.1}, \ref{f7.2}, \ref{f7.3}$)，相交的判据是 $0 \le s \le 1$ 且 $0 \le t \le 1$。

这段代码至少存在三个缺点:

不能处理平行线段。大多数应用程序需要知道线段是否重叠。
许多应用程序需要区分正常相交(proper)和非正常相交(improper)，就像我们在第 1 章中对三角分割所做的那样。在输出中区分这些情况会很有用。
尽管使用了浮点数，但在转换为双精度数之前，乘法仍然是用整数算术执行的。我们用一个简单的例子，说明此代码怎样溢出的。设四个端点为 $$ a = (-r, -r) $$ $$ b = (+r, +r) $$ $$ c = (+r, -r) $$ $$ d = (-r, +r) $$ 这些线段形成一个"X"形，相交于 $p = (0, 0)$。可以算出式$\ref{f7.1}, \ref{f7.2}$ 的分子均为 $-4r^2$。若 $r = 10^5$，该值为 $-4 \times 10^{10}$，这超出了 32 位所能表示的数据范围。在这种情况下，我的机器返回 $p = (-267702.8, -267702.8)$ 作为交点！

我们来逐一解决这三个问题。首先，我们将函数的返回值从布尔型改为 char，让它返回一个表示相交类型的编码。那些需要区分交点是 proper 还是 improper 的程序可以根据这个编码来判定。虽然具体使用的代码应根据应用程序而定，但以下几种基本能满足大多数需求

：

'e' 表示线段共线重叠，且共享一个点；'e' 代表 'edge'（边）。
'v' 表示某一线段的端点位于另一条线段上，但不满足 'e' 的情况；'v' 代表 'vertex'（顶点）。
'1' 表示线段正常相交，即它们共享一个点，且不满足 'v' 或 'e' 的情况；'1' 代表 TRUE（真）。
'0' 表示线段不相交，即它们没有共享任何点；'0' 代表 FALSE（假）。

需要注意的是，在这种分类方案中，两条共线线段仅共享一个端点的情况被归类为 'e'，尽管在某些语境下将其归为 'v' 可能更为合适。这种返回字符有点奇奇怪怪，个人更喜欢枚举类型。

然后，我们通过强制转换保证将乘法运算采用浮点数，使用 double 类型来扩大代码的适用范围。这便得到了代码 7.2 所示的代码。在讨论平行线段的情况之前，让我先指出这段代码的几个特点。检查 'v' 的情况时，使用的是 num 而不是除法后的 s 和 t。这避免了除法运算中可能出现的浮点误差。检查是否正确交集的条件是，$0 < s < 1$ 且 $0 < t < 1$。不等式不成立则不产生交集。

展开 Code 7.2 SegSegInt. >>>

        char SegSegInt( tPointi a, tPointi b, tPointi c, tPointi d, tPointd p )
        {
            double s, t;             /* The two parameters of the parametric eqns. */
            double num, denom;       /* Numerator and denominator of equations. */
            char code = '?';         /* Return char characterizing intersection. */
        
            denom = a[X] * (double)( d[Y] - c[Y] ) +
                    b[X] * (double)( c[Y] - d[Y] ) +
                    d[X] * (double)( b[Y] - a[Y] ) +
                    c[X] * (double)( a[Y] - b[Y] );
        
            /* If denom is zero, then segments are parallel: handle separately. */
            if (denom == 0.0)
                return ParallelInt(a, b, c, d, p);
        
            num = a[X] * (double)( d[Y] - c[Y] ) +
                  c[X] * (double)( a[Y] - d[Y] ) +
                  d[X] * (double)( c[Y] - a[Y] );
            if ( (num == 0.0) || (num == denom) ) code = 'v';
            s = num / denom;
        
            num = -( a[X] * (double)( c[Y] - b[Y] ) +
                     b[X] * (double)( a[Y] - c[Y] ) +
                     c[X] * (double)( b[Y] - a[Y] ) );
            if ( (num == 0.0) || (num == denom) ) code = 'v';
            t = num / denom;
        
            if ( (0.0 < s) && (s < 1.0) &&
                 (0.0 < t) && (t < 1.0) )
                code = '1';
            else if ( (0.0 > s) || (s > 1.0) ||
                      (0.0 > t) || (t > 1.0) )
                code = '0';
        
            p[X] = a[X] + s * ( b[X] - a[X] );
            p[Y] = a[Y] + s * ( b[Y] - a[Y] );
        
            return code;
        }

由于必须使用双精度浮点数，计算范围大大扩展。我只能在坐标值超过十亿时才能使其失效，在之前的溢出示例中，$r = 1234567809 ≈ 10^9$。这时溢出是正常的，因为 $-4r^2$ 已经超过了 $10^{18}$，这需要 60 位精度，超过了大多数机器上双精度尾数的长度。

最后，再来看平行线段，它由单独的函数 ParallelInt 处理。第 1 章中，我们使用函数 Between(Code 1.8)处理过共线重叠的情况。正好满足我们这里的需求。当两个线段重叠时，一个线段的端点位于另一个线段的端点之间。这里有一个小的简化。在三角分割代码中，我们使用 Between 检查共线，但在这里我们可以只检查一次。如果 $c$ 与 $ab$ 不共线，则平行线段 $ab$ 和 $cd$ 不相交。代码 7.3 给出了很直观的代码。其中，端点作为交点 $p$ 返回。可能，某些应用程序可能更喜欢返回重叠的中点。在 7.6 节中，我们将需要整个重叠线段。

展开 Code 7.3 ParallelInt. >>>

        char ParallelInt( tPointi a, tPointi b, tPointi c, tPointi d, tPointd p )
        {
            if ( !Collinear( a, b, c ) )
                return '0';
        
            if ( Between( a, b, c ) ) {
                AssignDi( p, c ); return 'e';
            }
            if ( Between( a, b, d ) ) {
                AssignDi( p, d ); return 'e';
            }
            if ( Between( c, d, a ) ) {
                AssignDi( p, a ); return 'e';
            }
            if ( Between( c, d, b ) ) {
                AssignDi( p, b ); return 'e';
            }
            return '0';
        }
        
        void AssignDi( tPointd p, tPointi a )
        {
            int i;
            for ( i = 0; i < DIM; i++ )
                p[i] = a[i];
        }
        
        bool Between( tPointi a, tPointi b, tPointi c )
        {
            tPointi ba, ca;
        
            /* If ab not vertical, check betweenness on x; else on y. */
            if ( a[X] != b[X] )
                return ((a[X] <= c[X]) && (c[X] <= b[X])) ||
                       ((a[X] >= c[X]) && (c[X] >= b[X]));
            else
                return ((a[Y] <= c[Y]) && (c[Y] <= b[Y])) ||
                       ((a[Y] >= c[Y]) && (c[Y] >= b[Y]));
        }

显然，通过改变的 $s$ 和 $t$值域，对该交点代码进行少量修改即可找到射线-线段(ray-segment),射线-射线(ray-ray),射线-直线(ray-line),直线-射线(line-ray)的交点。例如，允许第一个线段的 $s$ 可以取非负数，就可得到射线-线段的交点。

7.2 线段之间的交点(Segment-Segment Intersection)

7.2 线段之间的交点
(Segment-Segment Intersection)