线段交点/geos实现的分析以及我的实现

两线段(直线/射线)求交点是几何计算中的一个最基本的算法。虽然它大部分情况都不是影响程序效率的主要因素,但是它的执行次数可能非常的高,因此提升它的执行效率仍然具有价值。本文将分为两部分,一部分以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]\)的是否相交与交点。

发表回复

您的电子邮箱地址不会被公开。 必填项已用*标注