计算几何算法5. 直线、线段和平面相交(2D和3D)

计算几何 everyinch 11648℃ 0评论

直线和线段相交
平面相交
直线-平面相交
两平面相交
三个平面相交
实现
intersect2D_2Segments()
inSegment()
intersect3D_SegmentPlane()
intersect3D_2Planes()
参考文献

几何图元的相交是许多计算机图形学和建模应用中的重要组成部分。在这里我们来看一下最简单的 2D 和 3D 线性图元:直线、线段和平面的算法。

直线和线段相交

为了计算2D和3D中直线和线段的相交点,使用参数方程表示法是最好的。其它表示法在《计算几何算法2. 关于线和点到线的距离》一节中已经讨论过。本节研究如何从其它表示法转换到参数表示法。

在任意维度中,由两点P0、P1定义的直线参数方程式,当参数s为实数且u=P1-P0为直线的方向向量时可以表示为

P(s) = P0 + s (P1-P0) = P0 + s u,

当P(s)是有限线段P0P1上的一个点且0<s<1时,用表达式P(0)=P0, P(1)=P1,s = d(P0,P(s)) / d(P0,P1). 更进一步, 如果s<0,P(s)在线段P0之外,如果if s>1,P(s) 在线段P1的外侧。

clip_image001

让两条直线满足方程:P(s) = P0 + s (P1-P0) 和 Q(t) = Q0 + t (Q1-Q0),它们中的一条或者全部是有限线段或者射线。当且仅当方向共线时,直线平行,也就是说当两个向量u=P1-P0 和v=Q1-Q0 成线性关系时,即u=Cv,C为实数。 对于u=(uk)和v=(vk),t意味着所有比值uk/vk 相等,或者对于所有的k,u1/v1=uk/vk 都成立,这等同于所有满足u1vk-ukv1=0的情况。在2D中,对于u=(u1,u2) 和v=(v1,v2),u· v = u1v2-u2v1 = 0 这里u=(-u2,u1)垂直于u(⊥称为正交运算符)。。这种情况说明在欧几里得平面上当两条都垂直于相同的方向向量u时,两条直线是平行的。当这种情况成立时,两条直线既不重合也不相交。

通过验证点P0是否同时在直线和另外一条直线Q(t)上就可以很容易的验证两条直线是否重合。也就是说,存在一个参数t0 使等式P0=Q(t0)=Q0+t0v成立。如果w=(wk)=P0-Q0,对于某个参数t0存在等式w=t0v,等同于对于任意的参数k 存在等式w1vk-wkv1=0。在2D中,这是正交运算情况:w· v = w1v2-w2v1 = 0。如果这些条件成立,可以解出参数t0的值,和无限的直线式重合的。如果一条直线(而不是其它)是一条有限线段,它们重合相交。然而,如果两条直线都是有限线段,它们就可能重叠也可能不重叠。在这种情况下,解出t0t1 的值使得满足条件 P0=Q(t0) 和P1=Q(t1)。如果线段区间[t0, t1] 和[0,1]不相交,直线就不相交。否则,对区间进行相交运算得出[r0, r1] = [t0, t1] ∩ [0,1]。相交线段就是Q(r0)Q(r1) = P0P1 ∩Q0Q1。这在任何维度空间里都成立。

当两条直线不平行时,它们可能相交于特殊的一点。在2D欧几里得空间,它们总会相交。在更高维度,它们可能彼此错过不相交。但是如果它们的确相交,它们在2D平面上的线性投影也会相交。所以,一条线段简单的约束到两个坐标点(u 或者v不为0),通过P(sI)和Q(tI)计算出2D相交点,然后验证是否P(sI)=Q(tI) 对所有的坐标点都成立。为了计算出2D相交点I=P(sI)=Q(tI),考虑如下图中的两条直线和相关向量:

clip_image002

为了确定sI,在w=P0-Q0的情况下,我们得出向量等式P(s)-Q0 = w + su。在相交点,向量 P(s)-Q0 垂直于v,这等同于v· (w+su) = 0。求解等式,我们得到:

clip_image003

注意,只有像之前讨论过的直线平行时,v· u = 0。相似地,求解Q(tI),我们得出:

clip_image004

当u· v = – v· u时,分母和标志相同,只有在我们想同时知道sItI时,才计算分母。

如果两条直线中的一条是有限线段(或是一条射线),设为P0P1,然后当且仅当0<=sI<=1 (或对于射线有 sI>=0)时相交点在线段的内部。如果两条线都为有限线段,就要求解出参数sItI,它们都必须在区间[0,1]之间,线段才有可能相交。尽管这听起来足够简单,但是两条线段相交的代码需要一些技巧,因为有很多特殊的情况需要验证(参见我们的2D线段相交实现例子)。

平面相交

平面相交的情况已经在《计算几何算法4. 关于平面以及点到平面的距离》一节中描述过了。

直线平面相交

在3D中,直线 L或者平行于平面π或者与平面π相交于一点。假设 L满足下列参数方程 P(s) = P0 + s(P1-P0) = P0 + su,给定点V0在平面p上,法向量为 n。我们首先检查L是否平行于平面π,如果n·u=0,那就意味着直线法向量 u垂直于平面法向量n。如果成立,直线 L和平面π平行,即或者不相交,或者完全在平面p内。通过验证L上是否存在点P在平面π内,可以判断L和平面π是否相交,也就是说是否满足隐式直线方程:n·(P-V0) = 0。

如果直线和平面不平行,即直线L和平面π相交于点P(sI),利用类似于2D中两条直线相交的方法,我们可以计算出点P(sI)的值。如图所示:

clip_image005

在相交点,当w=P0-V0,向量P(s)-V0 = w+su垂直于n。这等同于点乘:n · (w+su) = 0。求解得出:

clip_image006

如果L 是从P0到P1的有限线段,我们必须检查sI是否满足 0<= sI <=1,从而证明线段和平面是否相交。对于正方向的射线,当sI >=0时,和平面相交。

两平面相交

在3D中,两平面 π1和π2或者平行或者相交于直线L。假设 πi (i=1,2)是由点Vi 和法向量ni所确定, 满足方程:ni · P + di = 0。只要平面π1和π2的法向量平行,它们就平行,等同于n1×n2 = 0。在软件中,我们经常用是否|n1×n2| < δ当除以δ时是否引起溢出,这是判断 n1 和 n2 是否平行的条件。如果不平行的话,u = n1×n2 >= δ > 0 是相交直线L的方向向量,因为 u 同时垂直于n1 和 n2,也就如图中所示的那样平行于两个平面。 如果 |u| 是个小数,最好将它缩放至|u| = 1,以保证u是单位方向向量。

clip_image007

在计算n1×n2 (3次相加和6次相乘)后,为了完全确定相交直线,我们仍然需要在上面找到一个特殊点;也就是说,我们需要找到同时在两平面上的点P0=(x0,y0,z0)。我们可以通过求解关于p1和p2的隐式方程来得到P0。但是这只有两个方程,却需要求解3个未知数,既然点P0 能够在任意用于其它自由度的一维直线L上,所以我们需要另外一个条件来求解出点P0。有以下几种方法可以做到这一点:

A)直接使用线性方程。将其中一坐标点设为0,即为z0=0,然后解出另外两个点。但是这只在L和平面z0=0相交时适用。当u =(ux,uy,uz)中的z值不为0时为真。所以,我们必须选择u的非零坐标点,然后设P0 = 0。我们应该选择拥有最大绝对值的坐标点,因为这会是最健壮的。如果uz != 0,P0=(x0,y0,0)在L上。求解两方程:a1x0+b1y0+d1=0 和a2x0+b2y0+d2=0,得出:

clip_image008

还有L的参数方程:

clip_image009

这里的分母和u非零的第三个坐标点是相等的。所以,忽视对非零坐标点的验证,这次求解的所有运算等于有5次相加和12次相乘。

B)直线相交点。如果知道平面上一条特定的直线(例如,平面上的两个点),和这条直线和其他平面相交,那么相交点”I”同时位于两个平面上。如此,这个点就位于两平面相交线上,L 的参数方程可以表示为:P(s) = I + s (n1×n2)。为了计算n1×n2和相交点(给定直线的),我们要进行11次加法运算和19次乘法运算。

在平面内构建一条和其他平面相交的直线的一种方法是,将此平面的法向量投射到另外一个平面上。给定一条直线,令其始终垂直于两平面的相交线。因此,n2在平面π1上的投影确定了和π2相交的直线。更具体来说,将两点0=(0,0,0)和n2=(nx2,ny2,nz2)分别投影到π1 (0) 和π1 (n2) 上。投影到π1上的直线是L1: Q(t) = π1 (0) + tπ1(n2),它和π2的相交线能被计算出来。在最有效的情况下,也就是当n1和n2 都为单位法向量并且常数p1(0) 被预先存储起来的情况下,整个运算包括17次相加运算和22次相乘运算。

C3平面相交于一点。当法向量n3 = n1×n2和d3=0 (意味着通过原点)时,可以选择用隐式方程式表达的第三个平面π3。这一直成立因为:(1) L垂直于π3这样就和π3相交,(2) 向量n1,n2和n3线性无关。也就是说平面π1,π2π3相交于唯一点P0,P0必然也在L上。当π3上的d3=0时,用3平面相交线的方程式,我们可以得出:

clip_image010

还有L的参数方程:

clip_image011.

这次求解共有11次相加和23次相乘。

n1n2为单位法向量且n1 · n2 = cos q(q是两向量之间的夹角)时,这个方程式可以被简化。然后,n1 · n1 = n2 · n2 = 1,和|n1×n2|2 = sin2 q = 1 – cos2 q = 1 – (n1 · n2)2 作为它的分母。最终,我们利用叉乘的右结合性,对于任意的3D向量ab,和c,等式a×(b×c) = (a· c)b – (a· b)c始终成立。应用到分子上,经过简化可以得出:(d2n1d1n2)×(n1×n2) = (d2 cos q – d1) n1 + (d1 cos q – d2) n2。于是,

clip_image012

这种解决方案描述了L上的点P0作为向量n1和n2的线性组合,这是通过原点的平面p3的基础生产算法。运算次数为11次相加+20次相乘。

最后要说的一点是,最有效率的解决方法是解决方法(A),这种方法计算出相交线的方程只用到了5次相加和12次相乘。

三平面相交

在3D中,三平面p1,p2 和p3有几下几种方式能够相交(或者不想交):

几何关系 交集 数学表达
三平面平行 nj×nk=0 (j,k=1,3)
clip_image014[2]三平面相交 一个平面 n1·V1=n1·V2=n1·V3
clip_image014[3]三平面不相交 n1·V1¹n1·V2¹n1·V3¹n1·V1
只有两个平面平行,第三个平面和其余的两个平面分别相交于一条直线 [注:两个平行的平面有可能相交] 2条平行的直线
[重合部分 => 1 条直线]
只有当nj×nk=0 (j¹k成立
没有两个平面平行,它们两两相交于三条直线。 nj×nk¹0 (j¹k
clip_image014[4]相交线互相平行 n1·(n2×n3)=0
clip_image017相交线相交 1条直线 从一条直线上测试一点
clip_image017[1]相交线不想交 3条平行的直线
clip_image014[5]相交线都不平行 一个特定的点 n1·(n2×n3)¹0

如图所示:

clip_image018

应该首先测试相交点的最普遍的情况,也就是说n1·(n2×n3)¹0,它可以排除其它的情况。当相交点是一个特定点时,它由以下公式给出:

clip_image019

可以通过P0是否满足p1p2p3的参数方程来验证。

然而,当分母n1·(n2×n3)非常小时,在计算性能上存在问题。在这种情况下,最好的办法是取其中两个平面的相交线,然后计算这条直线和第三个平面的相交点。

程序实现

//    Point and Vector with
//        coordinates {float x, y, z;}
//        operators for:
//            == to test equality
//            != to test inequality
//            Point  = Point ± Vector
//            Vector = Point - Point
//            Vector = Scalar * Vector    (scalar product)
//            Vector = Vector * Vector    (3D cross product)
//    直线、射线和线段由点定义{Point P0, P1;}
//        (直线是无穷的,射线和线段开始于P0)
//        (射线从P1向外延伸,但是线段结束于P1)
//    平面由一个点和一个方向向量定义{Point V0; Vector n;}
//===============================================================
#define SMALL_NUM  0.00000001 // anything that avoids division overflow
// 点乘(3D)
#define dot(u,v)   ((u).x * (v).x + (u).y * (v).y + (u).z * (v).z)
#define perp(u,v)  ((u).x * (v).y - (u).y * (v).x)  //正交运算(2D)

// intersect2D_2Segments():2D有限线段相交
//   输入:两条有限线段S1和S2
//   输出: *I0 = 相交点(如果存在)
//            *I1 = 相交线段的终点 [I0,I1] (如果存在)
//   返回值: 0=disjoint (no intersect)
//          1=intersect in unique point I0
//          2=overlap in segment from I0 to I1
int
intersect2D_Segments( Segment S1, Segment S2, Point* I0, Point* I1 )
{
    Vector    u = S1.P1 - S1.P0;
    Vector    v = S2.P1 - S2.P0;
    Vector    w = S1.P0 - S2.P0;
    float     D = perp(u,v);
    // 验证是否平行(包含一个点的情况)
    if (fabs(D) < SMALL_NUM) { // S1 和 S2 平行 if (perp(u,w) != 0 || perp(v,w) != 0) { return 0; // 不共线 } // they are collinear or degenerate // check if they are degenerate points float du = dot(u,u); float dv = dot(v,v); if (du==0 && dv==0) { // 线段点集 if (S1.P0 != S2.P0) // 在特殊点上 return 0; *I0 = S1.P0; // 在相同点上 return 1; } if (du==0) { // S1 作为一个单独的点 if (inSegment(S1.P0, S2) == 0) // 不在 S2上 return 0; *I0 = S1.P0; return 1; } if (dv==0) { // S2 作为一个单独的点 if (inSegment(S2.P0, S1) == 0) // 不在 S1上 return 0; *I0 = S2.P0; return 1; } // 它们为共线线段 – 部分重叠 (或者不重叠) float t0, t1; Vector w2 = S1.P1 - S2.P0; if (v.x != 0) { t0 = w.x / v.x; t1 = w2.x / v.x; } else { t0 = w.y / v.y; t1 = w2.y / v.y; } if (t0 > t1) {                  // t0 必须比 t1小
                float t=t0; t0=t1; t1=t;    // 如果t0 不比 t1小就交换
        }
        if (t0 > 1 || t1 < 0) {
            return 0;     // 不重叠
        }
        t0 = t0<0? 0 : t0; // 控制最小值为 0 t1 = t1>1? 1 : t1;              // 控制最大值为 1
        if (t0 == t1) {                 // 交集是一个点
            *I0 = S2.P0 + t0 * v;
            return 1;
        }

        // 在有效子线段内部分重叠
        *I0 = S2.P0 + t0 * v;
        *I1 = S2.P0 + t1 * v;
        return 2;
    }

    // 线段倾斜,可能相交于一点
    // 取得S1的相交参数
    float     sI = perp(v,w) / D;
    if (sI < 0 || sI > 1)               // 和 S1没有相交
        return 0;
    // get the intersect parameter for S2
    float     tI = perp(u,w) / D;
    if (tI < 0 || tI > 1)               // 和 S2没有相交
        return 0;
    *I0 = S1.P0 + sI * u;               // 计算 S1 的相交点
    return 1;
}
//===============================================================
// inSegment(): 验证点是否在线段内
//    输入:  a point P, and a collinear segment S
//    返回值: 1 = P is inside S
//            0 = P is not inside S
int
inSegment( Point P, Segment S)
{
    if (S.P0.x != S.P1.x) {    // S is not vertical
        if (S.P0.x <= P.x && P.x <= S.P1.x) return 1; if (S.P0.x >= P.x && P.x >= S.P1.x)
            return 1;
    }
    else {    // S is vertical, so test y coordinate
        if (S.P0.y <= P.y && P.y <= S.P1.y) return 1; if (S.P0.y >= P.y && P.y >= S.P1.y)
            return 1;
    }
    return 0;
}
//===============================================================
// intersect3D_SegmentPlane(): 线段和平面相交
//    Input:  S = a segment, and Pn = a plane = {Point V0; Vector n;}
//    Output: *I0 = the intersect point (when it exists)
//    Return: 0 = disjoint (no intersection)
//            1 = intersection in the unique point *I0
//            2 = the segment lies in the plane
int
intersect3D_SegmentPlane( Segment S, Plane Pn, Point* I )
{
    Vector    u = S.P1 - S.P0;
    Vector    w = S.P0 - Pn.V0;

    float     D = dot(Pn.n, u);
    float     N = -dot(Pn.n, w);

    if (fabs(D) < SMALL_NUM) {          // 线段和平面平行
        if (N == 0)                     // 线段在平面内
            return 2;
        else
            return 0;                   // 没有交集
    }
    // they are not parallel
    // compute intersect param
    float sI = N / D;
    if (sI < 0 || sI > 1)
        return 0;                       // 没有交集

    *I = S.P0 + sI * u;                 // 计算线段的相交点
    return 1;
}
//===============================================================
// intersect3D_2Planes(): 两平面相交
//    Input:  two planes Pn1 and Pn2
//    Output: *L = the intersection line (when it exists)
//    Return: 0 = disjoint (no intersection)
//            1 = the two planes coincide
//            2 = intersection in the unique line *L
int
intersect3D_2Planes( Plane Pn1, Plane Pn2, Line* L )
{
    Vector   u = Pn1.n * Pn2.n;         // 叉乘
    float    ax = (u.x >= 0 ? u.x : -u.x);
    float    ay = (u.y >= 0 ? u.y : -u.y);
    float    az = (u.z >= 0 ? u.z : -u.z);

    // test if the two planes are parallel
    if ((ax+ay+az) < SMALL_NUM) { // Pn1 和 Pn2 近似平行 // test if disjoint or coincide Vector v = Pn2.V0 - Pn1.V0; if (dot(Pn1.n, v) == 0) // Pn2.V0 在 Pn1内 return 1; // Pn1 and Pn2 重合 else return 0; // Pn1 and Pn2 不相交 } // Pn1 和 Pn2 相交于一条直线 // 首先定义叉乘的最大坐标值 int maxc; // 最大坐标值 if (ax > ay) {
        if (ax > az)
             maxc = 1;
        else maxc = 3;
    }
    else {
        if (ay > az)
             maxc = 2;
        else maxc = 3;
    }

    // 接下来,取得相交线上的一点
    // 将最大值置为0,求解出另外两个值
    Point    iP;               // 相交点
    float    d1, d2;           // 两个平面方程中的常数
    d1 = -dot(Pn1.n, Pn1.V0);  // 注: 事先在平面内定义
    d2 = -dot(Pn2.n, Pn2.V0);  // 同上

    switch (maxc) {            // 选择最大坐标
    case 1:                    // 相交于 x=0
        iP.x = 0;
        iP.y = (d2*Pn1.n.z - d1*Pn2.n.z) / u.x;
        iP.z = (d1*Pn2.n.y - d2*Pn1.n.y) / u.x;
        break;
    case 2:                    // 相交于 y=0
        iP.x = (d1*Pn2.n.z - d2*Pn1.n.z) / u.y;
        iP.y = 0;
        iP.z = (d2*Pn1.n.x - d1*Pn2.n.x) / u.y;
        break;
    case 3:                    // 相交于 z=0
        iP.x = (d2*Pn1.n.y - d1*Pn2.n.y) / u.z;
        iP.y = (d1*Pn2.n.x - d2*Pn1.n.x) / u.z;
        iP.z = 0;
    }
    L->P0 = iP;
    L->P1 = iP + u;
    return 2;
}
//===============================================================

参考文献

Foley, van Dam, Feiner & Hughes, 计算机图形学( C语言第二版), 3.12节 “线裁减” (1996)
Francis Hill, “正交点乘运算的乐趣” in 图形几何IV (1994)
Michael Mortenson, 计算机图形学手册: 几何和数学 (1990)
Joseph O’Rourke, 基于C的计算几何(第二版), 第7章 “查询和相交” (1998)



转载请注明:陈童的博客 » 计算几何算法5. 直线、线段和平面相交(2D和3D)

喜欢 (1)
发表我的评论
取消评论

表情

Hi,您需要填写昵称和邮箱!

  • 昵称 (必填)
  • 邮箱 (必填)
  • 网址
(2)个小伙伴在吐槽

无觅相关文章插件,快速提升流量