两线段（直线/射线）求交点是几何计算中的一个最基本的算法。虽然它大部分情况都不是影响程序效率的主要因素，但是它的执行次数可能非常的高，因此提升它的执行效率仍然具有价值。本文将分为两部分，一部分以geos为例，介绍主流的直线交点的判断与计算算法，另一部分给出我的实现，我相信它是一种更加简单和高效的算法，并且可以直接应用于更高的维度。
geos 中执行线段求交的类是geos::algorithm::LineIntersector::computeIntersection(const geom::Coordinate& p1, const geom::Coordinate& p2, const geom::Coordinate& p3, const geom::Coordinate& p4) 。求交的过程如下：
1）判断 p1p2构成的矩形与p3p4 构成的矩形是否相交，如果不想交，两个线段不存在交点。这个称为“短路测试”。可以有效的过滤不存在交点的情况。
// first try a fast test to see if the envelopes of the lines intersect
if (!Envelope::intersects(p1,p2,q1,q2))
{
return NO_INTERSECTION;
}
2) 计算p1,p2在 有向直线\( \vec {p_3p_4}\) 的方向信息，如果p1,p2均在\( \vec {p_3p_4}\) 的同一侧，那么两线段不相交，如果分列直线两侧，那么继续判断p3,p4与\( \vec {p_1p_2}\)的方位信息，如果也是分列两侧，那么能够认定两个线段存在交点。
// for each endpoint, compute which side of the other segment it lies
// if both endpoints lie on the same side of the other segment,
// the segments do not intersect
int Pq1=CGAlgorithms::orientationIndex(p1,p2,p3);
int Pq2=CGAlgorithms::orientationIndex(p1,p2,p4);
if ((Pq1>0 && Pq2>0) || (Pq1<0 && Pq2<0))
{
return NO_INTERSECTION;
}
//p1,p2 与 p3p4的方位关系判断略。
3）接下来，可以做一个判断，即p1p2,p3p4是否是共线的。如果p1,p2 都在\( \vec {p_3p_4}\) 上并且 p3,p4都在 \( \vec {p_1p_2}\)上，那么可以判断两线段共线。那么线段有两个交点。
bool collinear=Pq1==0 && Pq2==0 && Qp1==0 && Qp2==0;
4）如果不共线，那么两直线的存在且仅存在一个交点。geos 通过 intersection(const Coordinate& p1,const Coordinate& p2, const Coordinate& q1, const Coordinate& q2, Coordinate &intPt) const 方法实现两直线求交。
因为在实际运算中，p1,p2,p3,p4的值通常有非常大的offset，比如高斯坐标，在计算中很有可能造成溢出和精度问题，因此在运算时会减掉一个与上述点相近的值，这个过程称为“normalize”。geos的实现是减掉上述点的平均值：
double intMidX = (intMinX + intMaxX) / 2.0; double intMidY = (intMinY + intMaxY) / 2.0;
这样就不太好了。虽然double的值域非常大，但是这仍然不是一个很好的编程实践。
下面给出我的实现。
考虑线段pq可以用以下方式表达：$$l = p + \vec {pq}t ;t\in[0,1];$$
线段\( \vec {p_1p_2} = p_1 + \vec vt\)，线段 \( \vec {p_3p_4} = p_3 + \vec wr\),两线段的交点可以组成方程 \(p_1 + \vec vt = p_3 + \vec wr\) 。即交点满足两个线段的表达式。
如果方程有解，并且 \(t \in [0,1] r\in[0,1]\)代表两线段有唯一交点。
方程的矩阵形式如下：$$\begin{bmatrix}\vec v & -\vec w\end{bmatrix}\begin{bmatrix}t \\ r \end{bmatrix}=\begin{bmatrix}\vec {p_1p_3}\end{bmatrix}$$
1）如果左侧矩阵奇异，即v,w共线（线性相关），那么方程\(Ax=b\)无解或者存在无穷多个解。表现在几何上为两线段所在直线无交点(2D下表现为平行)或者存在无数个交点（共线）。从线性代数角度上讲，
当\(rank(A)=rank(A|b)\)时，方程存在无穷多解。
当\(rank(A) < rank(A|b)\)时，方程无解。
2）如果左侧矩阵可逆，那么直接求解t,r，检查t,r是否在[0,1]区间内，并且任意将t,r带回，可以求线段交点。
下列代码给出线段\(\vec{p1p2},\vec{p3p4}\)交点的计算代码(不计算两线段共线的情况)，借助线性代数算法库Eigen 实现:
#include <Eigen/Dense>
bool Intersection::line_segment_intersection(const Coordinate &p1,
const Coordinate &p2,
const Coordinate &p3,
const Coordinate &p4,
Coordinate &p)
{
auto v = p2 - p1;
auto w = p4 - p3;
Eigen::Matrix2d A;
A<<v.x,-w.x,v.y,-w.y;
//A is not inversible
if(0 == A.determinant())
return false;
Eigen::Vector2d b;
b <<p3.x-p1.x,p3.y-p1.y;
Eigen::Vector2d x = A.colPivHouseholderQr().solve(b);
double t = x(0);
double r = x(1);
//if intersection is not on line segment
if(t<0||t>1||r<0||r>1)
return true;
p = p1 + v*t;
return true;
}
两种方法的对比：
1）geos方法的交点是否存在和求交点的逻辑是分开的。每次计算点与直线的方位关系需要一次矢量叉积计算，这一样的操作要进行四次(存在交点)。
我给出的方法只需要解一个线性方程，运用crammer法则，只需要解三个\(2\times2\)行列式。并且直接可通过解的情况判断 不想交/共线/有交点。
2）我给出的方法同样适合高维线段的交点计算。geos因为要判断点在直线的哪一侧，会受限于二维运算。
3）根据方向向量系数值域的不同，我给出的方法可以表示 线段\(t\in[0,1]\) 射线\(t\in[0,\infty]\),直线\(t\in[-\infty,+\infty]\)的是否相交与交点。