直线和线段相交
平面相交
直线-平面相交
两平面相交
三个平面相交
实现
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+ su,
当P(s)是有限线段P0P1上的一个点且01,P(s) 在线段P1的外侧。
让两条直线满足方程: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的值,和无限的直线式重合的。如果一条直线(而不是其它)是一条有限线段,它们重合相交。然而,如果两条直线都是有限线段,它们就可能重叠也可能不重叠。在这种情况下,解出t0和t1 的值使得满足条件 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),考虑如下图中的两条直线和相关向量:
为了确定sI,在w=P0-Q0的情况下,我们得出向量等式P(s)-Q0 = w + su。在相交点,向量 P(s)-Q0 垂直于v⊥,这等同于v⊥· (w+su) = 0。求解等式,我们得到:
注意,只有像之前讨论过的直线平行时,v⊥· u = 0。相似地,求解Q(tI),我们得出:
当u⊥· v = – v⊥· u时,分母和标志相同,只有在我们想同时知道sI 和 tI时,才计算分母。
如果两条直线中的一条是有限线段(或是一条射线),设为P0P1,然后当且仅当0<=sI<=1 (或对于射线有 sI>=0)时相交点在线段的内部。如果两条线都为有限线段,就要求解出参数sI 和 tI,它们都必须在区间[0,1]之间,线段才有可能相交。尽管这听起来足够简单,但是两条线段相交的代码需要一些技巧,因为有很多特殊的情况需要验证(参见我们的2D线段相交实现例子)。
平面相交
直线–平面相交
在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)的值。如图所示:
在相交点,当w=P0-V0,向量P(s)-V0 = w+su垂直于n。这等同于点乘:n · (w+su) = 0。求解得出:
如果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是单位方向向量。
在计算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,得出:
还有L的参数方程:
这里的分母和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次相乘运算。
(C)3平面相交于一点。当法向量n3 = n1×n2和d3=0 (意味着通过原点)时,可以选择用隐式方程式表达的第三个平面π3。这一直成立因为:(1) L垂直于π3这样就和π3相交,(2) 向量n1,n2和n3线性无关。也就是说平面π1,π2和π3相交于唯一点P0,P0必然也在L上。当π3上的d3=0时,用3平面相交线的方程式,我们可以得出:
还有L的参数方程:
.
这次求解共有11次相加和23次相乘。
当n1和n2为单位法向量且n1· n2 = cos q(q是两向量之间的夹角)时,这个方程式可以被简化。然后,n1· n1 = n2· n2 = 1,和|n1×n2|2 = sin2 q = 1 – cos2 q = 1 – (n1· n2)2 作为它的分母。最终,我们利用叉乘的右结合性,对于任意的3D向量a,b,和c,等式a×(b×c) = (a· c)b – (a· b)c始终成立。应用到分子上,经过简化可以得出:(d2n1–d1n2)×(n1×n2) = (d2 cos q –d1) n1 + (d1 cos q –d2) n2。于是,
这种解决方案描述了L上的点P0作为向量n1和n2的线性组合,这是通过原点的平面p3的基础生产算法。运算次数为11次相加+20次相乘。
最后要说的一点是,最有效率的解决方法是解决方法(A),这种方法计算出相交线的方程只用到了5次相加和12次相乘。
三平面相交
在3D中,三平面p1,p2 和p3有几下几种方式能够相交(或者不想交):
几何关系
交集
数学表达
三平面平行
nj×nk=0 (j,k=1,3)
三平面相交
一个平面
n1·V1=n1·V2=n1·V3
三平面不相交
无
n1·V1¹n1·V2¹n1·V3¹n1·V1
只有两个平面平行,第三个平面和其余的两个平面分别相交于一条直线 [注:两个平行的平面有可能相交]
2条平行的直线
[重合部分 => 1 条直线]
只有当nj×nk=0 (j¹k)成立
没有两个平面平行,它们两两相交于三条直线。
nj×nk¹0 (j¹k)
相交线互相平行
n1·(n2×n3)=0
相交线相交
1条直线
从一条直线上测试一点
相交线不想交
3条平行的直线
相交线都不平行
一个特定的点
n1·(n2×n3)¹0
如图所示:
应该首先测试相交点的最普遍的情况,也就是说n1·(n2×n3)¹0,它可以排除其它的情况。当相交点是一个特定点时,它由以下公式给出:
可以通过P0是否满足p1,p2 和p3的参数方程来验证。
然而,当分母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)