直线射线线段的相交判断

本文主要介绍计算机图形学中线型对象的一些几何关系判断，包括直线、射线和线段。

1. 二维空间

1.1 直线与直线相交

在欧几里得二维空间，两条直线的关系有以下几种：
1. 相交（交点只有一个）
2. 平行（没有任何交点）
3. 共线（两条直线重叠在一起）

直线的表示一般有参数式和隐式的点法表达方式，二者之间很容易互相转换到另一种形式。一般来说在计算直线相交时，使用参数式的方式更加方便。参数式的表达方式如下：

$$L = P + \vec at$$
P点是直线上的某一点，向量 $\vec a$通过直线上两点的向量得到。事实上射线和线段也可以使用同样的参数表达式，只不过在t的取值范围上有所不同。

有了上面的基础，开始推导计算的公式，设两条直线的参数方程是：
$P_0 + s\vec d_0$ 和 $P_1 + t\vec d_1$， 如果它们存在交点，那么交点应该满足：
$P_0 + s\vec d_0 = P_1 + t\vec d_1$，这个等式通过x和y坐标展开，实际上是一个二元一次的方程组，可以解出来s和t这两个未知数。为了简化形式，定义操作 $Kross（\vec v_0, \vec v_1) = v_0.x*v_1.y - v_1.x*v_0.y$ 和 $\vec \Delta = P_1 - P_0$ 解方程可以得到：

$$Kross(\vec d_0, \vec d_1) * s = Kross(\vec \Delta, \vec d_1) \qquad Kross(\vec d_0, \vec d_1) * t = Kross(\vec \Delta, \vec d_0)$$
得到上述表达式后，可以展开讨论：

1. 如果 $Kross(\vec d_0, \vec d_1)$结果不为0， 那么可以知道两条直线相交且只有一个交点，可以得到$s = Kross(\vec \Delta, \vec d_1) / Kross(\vec d_0, \vec d_1)$ 以及 $t = Kross(\vec \Delta, \vec d_0) / Kross(\vec d_0, \vec d_1)$

2. 如果 $Kross(\vec d_0, \vec d_1)$结果为0， 那么这两条直线要么是平行，要么就是重合（同一条直线），由于$Kross(\vec d_0, \vec d_1) = 0$，如果 $Kross(\vec \Delta, \vec d_0) = 0$（或者 $Kross(\vec \Delta, \vec d_1) = 0$)，可以得到这两条直线的方程只差一个系数关系。也就是它们是同一条直线（或者说是重合在一起），否则如果$Kross(\vec \Delta, \vec d_0) ！= 0$（或者 $Kross(\vec \Delta, \vec d_1) ！= 0$），它们就是平行的关系（不重合）

理论上的计算上面已经讨论结束，但是在编码的过程中，还有一些可以探讨的地方。（本文仅实现直线与直线以及线段与选段相交的算法，其他算法的实现读者可以参考上述两者的实现完成）

1.1.1 相交判断与求交点

有一些应用需要判断两条直线是否相交，但是并不要求求出交点，这个时候我们就没有必要计算到求出交点那一步，只要当确实需要知道交点的时候才去求，减少计算量。考虑到两个三维向量的叉乘，
$v1（x_0, y_0, 0)， v2（x_1, y_1, 0)$ ，那么它们叉乘的结果是：$（0， 0， x_0y_1-x_1,y_0)$，根据叉乘求模长的公式，可以得到：

$|v1 \times v2| = |v1||v2|sin\theta = |(0,0,x_0y_1-x_1y_0)|$
$|Kross(v1,v2)| = |v1||v2||sin\theta|$

要判断$Kcross(v1,v2) = 0$，需要等式右边这三项如果有一项接近0，也就是说，当我们的向量的模长很小时，得到的Kcross也很有可能为0，即使这两个向量并不平行（或者重合）。可行的一种方案是在判断的时候将v1和v2单位化（这样只有当$sin\theta$很小时[需要注意在$\theta$很小时，$sin\theta 近似等于\theta$]，两条直线才可能判断出平行或者重合的关系，这正好是与实际相符的），另一种方案是使用相对的长度来判断，也就是

$$\frac {||Kross(\vec d_0, \vec d_1)||}{||\vec d1|| \ ||\vec d2||} = |sin\theta| < \epsilon$$
在计算中为了减少做开方的运算，可以对等式两边取平方：
$${||Kross(\vec d_0, \vec d_1)||}^2 < \epsilon^2 || {\vec d_0} ||^2 \ || {\vec d_1} ||^2$$
代码如下，分为相交判断以及求交点判断两个函数（如果仅仅是测试是否相交而不需要求交点位置，可以使用test函数）

//find intersect
//return value : 
//1(one intersection) 0 (parallel,no intersection) 2(same line)
//@param point: if return value is one, point is the intersection point.
int  findLine2LineIntersection2D(
    const Vec2d& p0, const Vec2d& p1, 
    const Vec2d& p2, const Vec2d& p3, Vec2d& point)
{
    const double epsilon = 1e-7;

    Vec2d d0 = p1 - p0;
    Vec2d d1 = p3 - p2;
    Vec2d diff = p2 - p0;

    double krossd0d1    = d0.x()*d1.y() - d1.x()*d0.y();
    double sqrkrossd0d1 = krossd0d1*krossd0d1;
    double sqrlen0 = d0.length2();
    double sqrlen1 = d1.length2();

    if (sqrkrossd0d1 > epsilon * sqrlen0 * sqrlen1)
    {
        double s = diff.x() * d1.y() - d1.x() * diff.y();
        point = p0 + s * d0;
    }

    double krossdiffd0 = diff.x()*d0.y() - d0.x()*diff.y();
    double sqrkrossdiffd0 = krossdiffd0*krossdiffd0;
    double sqrlendiff = diff.length2();
    if (sqrkrossdiffd0 > epsilon * sqrlen0 * sqrlendiff)
    {
        return 0;
    }

    return 2;
}

//test intersect
//return value : 
//1(one intersection) 0 (parallel,no intersection) 2(same line)
int  testLine2LineIntersection2D(
    const Vec2d& p0, const Vec2d& p1, 
    const Vec2d& p2, const Vec2d& p3)
{
    const double epsilon = 1e-7;

    //one intersection
    Vec2d d0 = p1 - p0;
    Vec2d d1 = p3 - p2;
    double krossd0d1 = d0.x()*d1.y() - d1.x()*d0.y();
    double sqrKrossd0d1 = krossd0d1 * krossd0d1;
    double sqrL0 = d0.length2();
    double sqrL1 = d1.length2();
    if (sqrKrossd0d1 > epsilon * sqrL0 * sqrL1)
    {
        return 1;
    }

    //parallel
    Vec2d diff = p2 - p0;
    double sqrdiff = diff.length2();
    double krossdiffd0 = diff.x()*d0.y() - d0.x()*diff.y();
    double sqrkrossDiffd0 = krossdiffd0 * krossdiffd0;
    if (sqrkrossDiffd0 > epsilon * sqrdiff * sqrL0)
    {
        return 0;
    }

    //same line
    return 2;
}

1.2 直线与射线相交

直线与射线相交的判断和1.1中直线与直线的判断是类似的。射线的方程同样可以用直线的方程来表达，只不过此时射线的t是有取值范围的，它的取值范围是$[0, +\infty]$，在求出t值之后，比较t值是否在$[0, +\infty]$之间，即可得到两者之间的关系。

1.3 直线与线段相交

直线和线段相交，同样可以将直线方程写成1.1中的直线方程，只不过线段的取值范围有所限制，如果已知线段的起点为$P_0$，终点为$P_1$，可以使用下面的表达形式：

$$线段方程： C + s* \vec d \\ 其中： C = \frac{(P_0 + P_1)}{2} \\ \vec d = \frac {(P_1 - P_0)}{|P1 - P0|} \\ s \in [\frac{|P_1 - P_0|}{2}, \frac{|P_1-P_0|}2]$$
在使用公式计算出1.1种的t和s之后，比较s是否落在这个区间，从而判断出是否相交。

1.4 射线与射线相交

射线与射线相交有下面的4种情形：
1. 两射线有唯一交点
2. 两射线平行（不重合），没有交点
3. 两射线重合但是方向相反，有一个公共的相交线段
4. 两射线重合但是方向相同
其中第3种情况，不存在于直线和直线相交的情况。
处理的方式还是在直线和直线相交的基础上，通过限制t和s的取值来讨论。

1.5 射线与线段相交

射线与线段相交，可以同样使用直线和直线相交的方式来判断，射线与线段相交的情况有以下几种：
1. 唯一交点
2. 没有交点
3. 重合，交点是公共的相交线段

1.6 线段与线段相交

线段和线段相交也是同样的方式，通过直线和直线相交判断，之后限制参数的取值范围，来计算出公共的交点。线段和线段相交在求重叠区域的代码的代码如下，对这段代码稍微解释一下：

    double s0 = d0*diff / sqrlen0;
    double s1 = s0 + d0*d1 / sqrlen0;

    double smin = std::min(s0, s1);
    double smax = std::max(s0, s1);

    double w[2];

    int imax = findOverlappedInterval(0.0, 1.0, smin, smax, w);

如果两条直线重合，那么对于第一条直线上的任意一点，它必然满足另一条直线的方程（毕竟二者是同一条直线），代码中的计算逻辑如下：
线段1的参数方程：$P_0 + \vec t d_0$
线段2的参数方程：$P_1 + \vec s d_1$
参数t和s的取值范围都是[0, 1]，也就是说两条当 t和s的取值小于0时，在线段的左侧（起点的延长线上），当t和s的取值大于1时，在线段的右侧（终点的延长线上）。现在我们知道过线段1和线段2的直线是重合的，那么在线段2所属直线上的任意一点必然在线段1所属直线上，我们现在计算两个点（线段2的起点和线段2的终点），看这两个点在线段1所属直线上的t的取值，也就是说我们计算：
$P_1 = P_0 + t_0 \vec d_0$，可以计算得到$t_0 = \frac {(P_1 - P_0) \cdot \vec d_0}{\vec d_0^2}$，同理可以计算得到 $t_1 = t_0 + \vec d_0 \cdot \vec d_1 / |d_0^2|$，也就是说线段二上面的两点在线段1所在直线上的区段是$[t_0，t_1]$，我们和线段一的区段[0,1]求一个公共的交集，就得到了两个线段的重叠区域了。

---

总结来说，二维线性对象（直线、射线、线段）的相交判断都是以直线和直线相交的基本理论为基础的，后续根据射线和线段是直线的某种特殊情形（限制了自变量的取值范围），来分析得出最后的交点。
完整的线段与线段求交的代码：

// return value:
// 0 (parallel)
// 1 (only one intersection)
// 2 (overlapped)
int testSegment2SegmentIntersection2D(
    const Vec2d& p0, const Vec2d& p1, 
    const Vec2d& p2, const Vec2d& p3)
{
    const double epsilon = 1e-7;
    Vec2d d0 = p1 - p0;
    Vec2d d1 = p3 - p2;
    Vec2d diff = p2 - p0;

    //only one intersection
    double krossd0d1 = d0.x()*d1.y() - d1.x()*d0.y();
    double sqrkrossd0d1 = krossd0d1*krossd0d1;
    double sqrlen0 = d0.length2();
    double sqrlen1 = d1.length2();
    if (sqrkrossd0d1 > epsilon*sqrlen0*sqrlen0)
    {
        double s = diff.x()*d1.y() - d1.x()*diff.y();
        if (s < 0 || s > 1)
        {
            return 0;
        }

        double t = diff.x()*d0.y() - d0.x()*diff.y();
        if (t < 0 || t > 1)
        {
            return 0;
        }

        return 1;
    }

    //parallel
    double sqrlendiff = diff.length2();
    double krossdiffd0 = diff.x()*d0.y() - d0.x()*diff.y();
    double sqrkrossdiffd0 = krossdiffd0*krossdiffd0;

    if (sqrkrossdiffd0 > epsilon*sqrlen0*sqrlendiff)
    {
        return 0;
    }

    //same line, possibly overlapped
    double s0 = d0*diff / sqrlen0;
    double s1 = s0 + d0*d1 / sqrlen0;

    double smin = std::min(s0, s1);
    double smax = std::max(s0, s1);

    double w[2];

    return findOverlappedInterval(0.0, 1.0, smin, smax, w);
}

// return value:
// 0 (parallel)
// 1 (only one intersection)
// 2 (overlapped)
int findSegment2SegmentIntersection2D(
    const Vec2d& p0, const Vec2d& p1, 
    const Vec2d& p2, const Vec2d& p3, Vec2d point[2])
{
    const double epsilon = 1e-7;
    Vec2d d0 = p1 - p0;
    Vec2d d1 = p3 - p2;
    Vec2d diff = p2 - p0;

    //only one intersection
    double krossd0d1 = d0.x()*d1.y() - d1.x()*d0.y();
    double sqrkrossd0d1 = krossd0d1*krossd0d1;
    double sqrlen0 = d0.length2();
    double sqrlen1 = d1.length2();
    if (sqrkrossd0d1 > epsilon*sqrlen0*sqrlen0)
    {
        double s = diff.x()*d1.y() - d1.x()*diff.y();
        if (s < 0 || s > 1)
        {
            return 0;
        }

        double t = diff.x()*d0.y() - d0.x()*diff.y();
        if (t < 0 || t > 1)
        {
            return 0;
        }

        point[0] = p0 + s*d0;

        return 1;
    }

    //parallel
    double sqrlendiff = diff.length2();
    double krossdiffd0 = diff.x()*d0.y() - d0.x()*diff.y();
    double sqrkrossdiffd0 = krossdiffd0*krossdiffd0;

    if (sqrkrossdiffd0 > epsilon*sqrlen0*sqrlendiff)
    {
        return 0;
    }

    //same line, possibly overlapped
    double s0 = d0*diff / sqrlen0;
    double s1 = s0 + d0*d1 / sqrlen0;

    double smin = std::min(s0, s1);
    double smax = std::max(s0, s1);

    double w[2];

    int imax = findOverlappedInterval(0.0, 1.0, smin, smax, w);
    for (int i = 0; i < imax; ++i)
    {
        point[i] = p0 + w[i] * d0;
    }
    return imax;
}

2. 三维空间

2.1 直线与直线相交

三维空间里直线与直线相交关系相比二维空间更复杂一些，三维空间的关系包括以下几种：

1. 共面
   1.1 相交
   1.2 平行(不重合)
   1.3 重合
2. 异面(Skew Lines)
   也就是说三维直线的相交关系判断里面多了异面这种情况

推导两三维空间直线相交关系的计算：
直线以参数形式描述：

$$L_1(t) = P_1 + \vec v_1 t \quad L_2(t) = P_2 + \vec v_2 s$$
如果两条直线相交，那么交点满足： $P_1 + \vec v_1 t = P_2 + \vec v_2 s$，也就是： $\vec v_1t = (P_2 - P_1) + \vec v_2 s$，等式两边叉乘 $\vec v_2$，有： $(\vec v_1 \times \vec v_2)t = (P_2 - P_1) \times \vec v_2$，两边再点乘 $(\vec v_1 \times \vec v_2)$， 得到：
$$(\vec v_1 \times \vec v_2) \cdot (\vec v_1 \times \vec v_2) = (（P_2-P_1) \times \vec v_2) \cdot (\vec v_1 \times \vec v_2)$$
根据叉乘和点乘运算的规则
$$(\vec P \times \vec Q) \cdot \vec R) = \begin{vmatrix} P_x & P_y & P_z \\ Q_x & Q_y & Q_z \\ R_x & R_y & R_z \\ \end{vmatrix}$$
正好是三者组成3x3行列式的值，也就是：
$$t = \frac { Det \{(P_2 - P_1), \vec v_2, \vec v_1 \times \vec v_2 \} }{ |\vec v_1 \times \vec v_2|^2 }$$
同理可得s的值
$$s = \frac { Det \{(P_2 - P_1), \vec v_1, \vec v_1 \times \vec v_2 \} }{ |\vec v_1 \times \vec v_2|^2 }$$

观察分母项目，当 $|\vec v_1 \times \vec v_2|^2$取值为0时，两条空间直线的关系是平行或者重合。
如果$（P_2 - P_1) \cdot (\vec v_2 \times \vec v_1) \ne 0$，那么两条直线就是异面的关系。否则两条直线就是相交的关系。

如果两条直线是异面的关系，那么通过上面s和t所求得的点，就是两条直线最短距离的两个点。详细证明参考参考文献6

参考文献：

1. 知乎：向量点乘的 x1x2 + y1y2 与 |a||b|cosθ 是如何联系起来的
2. Geometric Tool Engine
3. wiki:Line–line intersection
   4.MathWorld:Skew Lines
4. 《Graphics Gems I》 Intersection of two lines in three-space
5. 《Mathematics for3D Game Programming and Computer Graphics》5.1.2 Distance Between Two Lines
6. Intersections of Lines and Planes