The intersection of geometric primitives is a fundamental construct in many computer graphics and modeling applications ([Foley et al, 1996], [O'Rourke, 1998]). Here we look at the algorithms for the simplest 2D and 3D linear primitives: lines, segments and planes.
Line and Segment Intersections
For computing intersections of lines and segments in 2D and 3D, it is best to use the parametric equation representation for lines. Other representations are discussed in Algorithm 2 about the Distance of a Point to a Line. It is shown there how to convert from other representations to the parametric one.
In any dimension, the parametric equation of a line defined by two points P0 and P1 can be represented as:
Ps=P0+s_P1-P0=P0+su,
where the parameter s is a real number and u=P1-P0 is a line direction vector. Using this representation P(0)=P0, P(1)=P1, and when s.ge-0.le-1, P(s) is a point on the finite segment P0P1 where s is the fraction of P(s)'s distance along the segment. That is, s = d(P0,P(s)) / d(P0,P1). Further, if s < 0 then P(s) is outside the segment on the P0 side, and if s > 1 then P(s) is outside the segment on the P1 side.
Pic_parametric-s
Let two lines be given by: Ps=P0+s_P1-P0=P0+su and Qt=Q0+t_Q1-Q0=Q0+tv, either or both of which could be infinite, a ray, or a finite segment.
Parallel Lines
These lines are parallel when and only when their directions are collinear, namely when the two vectors u=P1-P0 and v=Q1-Q0 are linearly related as u = av for some real number a. For u=(ui) and v=(vi), this means that all ratios ui_div_vi have the value a, or that u1_div_v1=ui_div_vi for all i. This is equivalent to the conditions that all u1vi-uiv1=0. In 2D, with u=(u1,u2) and v=(v1,v2), this is the perp product [Hill, 1994] condition that perp-u_dot_v=u1v2-u2v1=0 where perp_u=-u2,u1, the perp operator, is perpendicular to u. This condition says that two vectors in the Euclidean plane are parallel when they are both perpendicular to the same direction vector perp_u. When true, the two associated lines are either coincident or do not intersect at all.
Coincidence is easily checked by testing if a point on one line, say P0, also lies on the other line Q(t). That is, there exists a number t0 such that: P0 = Q(t0) = Q0 + t0v. If w = (wi) = P0 – Q0, then this is equivalent to the condition that w = t0v for some t0 which is the same as w1vi-wiv1=0 for all i. In 2D, this is another perp product condition: perp-w_dot_v=w1v2-w2v1=0. If this condition holds, one has t0=w1_div_v1, and the infinite lines are coincident. And if one line (but not the other) is a finite segment, then it is the coincident intersection. However, if both lines are finite segments, then they may (or may not) overlap. In this case, solve for t0 and t1 such that P0 = Q(t0) and P1 = Q(t1). If the segment intervals [t0,t1] and [0,1] are disjoint, there is no intersection. Otherwise, intersect the intervals (using max and min operations) to get r0,r1=t0,t1_inter_0,1. Then the intersection segment is Qr0Qr1=P0P1_inter_Q0Q1. This works in any dimension.
Non-Parallel Lines
When the two lines or segments are not parallel, they might intersect in a unique point. In 2D Euclidean space, infinite lines always intersect. In higher dimensions they usually miss each other and do not intersect. But if they intersect, then their linear projections onto a 2D plane will also intersect. So, one can simply restrict to two coordinates, for which u and v are not parallel, compute the 2D intersection point I at P(sI) and Q(tI) for those two coordinates, and then test if P(sI) = Q(tI) for all coordinates. To compute the 2D intersection point, consider the two lines and the associated vectors in the diagram:
Pic_line-line
To determine sI, we have the vector equality P(s)-Q0=w+su where w=P0-Q0. At the intersection, the vector P(sI)-Q0 is perpendicular to perp_v, and this is equivalent to the perp product condition that perp-v_dot_(w+su)=0. Solving this equation, we get:
Note that the denominator perp-v_dot_u=0 only when the lines are parallel as previously discussed. Similarly, solving for Q(tI), we get:
The denominators are the same up to sign, since perp-u_dot_v=-perp-v_dot_u=0, and should only be computed once if we want to know both sI and tI. However, knowing either is enough to get the intersection point I = P(sI) = Q(tI).
Further, if one of the two lines is a finite segment (or a ray), say P0P1, then the intersect point is in the segment only when s-I.ge-0.le-1 (or s-I.ge-0 for a ray). If both lines are segments, then both solution parameters, sI and tI, must be in the [0,1] interval for the segments to intersect. Although this sounds simple enough, the code for the intersection of two segments is a bit delicate since many special cases need to be checked (see our implementation intersect2D_2Segments()).
Plane Intersections
Planes are represented as described in Algorithm 4, see Planes.
Line-Plane Intersection
In 3D, a line L is either parallel to a plane P or intersects it in a single point. Let L be given by the parametric equation: P(s)=P0+s(P1-P0)=P0+su, and the plane P be given by a point V0 on it and a normal vector n=(a,b,c). We first check if L is parallel to P by testing if n_dot_u=0, which means that the line direction vector u is perpendicular to the plane normal n. If this is true, then L and P are parallel and either never intersect or else L lies totally in the plane P . Disjointness or coincidence can be determined by testing whether any one specific point of L, say P0, is contained in P , that is whether it satisfies the implicit line equation: n_dot_(P0-V0)=0.
If the line and plane are not parallel, then L and P intersect in a unique point P(sI) which is computed using a method similar to the one for the intersection of two lines in 2D. Consider the diagram:
Pic_line-plane
At the intersect point, the vector P(s)-V0=w+su is perpendicular to n, where w=P0-V0. This is equivalent to the dot product condition: n_dot_(w+su)=0. Solving we get:
If the line L is a finite segment from P0 to P1, then one just has to check that s-I.ge-0.le-1 to verify that there is an intersection between the segment and the plane. For a positive ray, there is an intersection with the plane when s-I.ge-0.
Intersection of 2 Planes
In 3D, two planes P1 and P2 are either parallel or they intersect in a single straight line L. Let P i (i = 1,2) be given by a point Vi and a normal vector ni, and have an implicit equation: ni · P + di = 0, where P = (x, y, z). The planes P1 and P2 are parallel whenever their normal vectors n1 and n2 are parallel, and this is equivalent to the condition that: n1 x n2 = 0. In software, one would test if n1_cross_n2.le.delta where division by delta.gt.0 would cause overflow, and uses this as the robust condition for n1 and n2 to be parallel in practice. When not parallel, u=n1_cross_n2 is a direction vector for the intersection line L since u is perpendicular to both n1 and n2, and thus is parallel to both planes as shown in the following diagram. If |u| is small, it is best to scale it so that |u| = 1, making u a unit direction vector.
Pic_2-planes
After computing n1_cross_n2 (3 adds + 6 multiplies), to fully determine the intersection line, we still need to find a specific point on it. That is, we need to find a point P0=(x0,y0,z0) that lies in both planes. We can do this by finding a common solution of the implicit equations for P1 and P2. But there are only two equations in the 3 unknowns since the point P0 can lie anywhere on the 1-dimensional line L. So we need another constraint to solve for a specific P0. There are a number of ways this could be done:
(A) Direct Linear Equation. One could set one coordinate to zero, say z0 = 0, and then solve for the other two. But this will only work when L intersects the plane z0 = 0. This is will be true when the z-coordinate uz of u=(ux,uy,uz) is nonzero. So, one must first select a nonzero coordinate of u, and then set the corresponding coordinate of P0 to 0. Further, one should choose the coordinate with the largest absolute value, as this will give the most robust computations. Suppose that uz.ne.0, then P0=(x0,y0,0) lies on L. Solving the two equations: a1x0+b1y0+d1=0 and a2x0+b2y0+d2=0, one gets:
Eqn_2planes-A-P0
and the parametric equation of L is:
Eqn_2planes-A
The denominator here is equal to the non-zero 3rd coordinate of u. So, ignoring the test for a large nonzero coordinate, and counting division as a multiplication, the total number of operations for this solution = 5 adds + 13 multiplies.
(B) Line Intersect Point. If one knows a specific line in one plane (for example, two points in the plane), and this line intersects the other plane, then its point of intersection, I, will lie in both planes. Thus, it is on the line of intersection for the two planes, and the parametric equation of L is: P(s) = I + s (n1 x n2). To compute n1_cross_n2 and the intersection point (given the line), the total number of operations = 11 adds + 19 multiplies.
One way of constructing a line in one plane that must intersect the other plane is to project one plane's normal vector onto the other plane. This gives a line that must always be orthogonal to the line of the planes' intersection. So, the projection of n2 on P1 defines a line that intersects P2 in the sought for point P0 on L. More specifically, project the two points 0 = (0,0,0) and n2 = (nx2, ny2, nz2) to P1(0) and P1(n2) respectively. Then the projected line in P1 is L1: Q(t) = P1(0) + t (P1(n2) – P1(0)), and intersection of it with P2 can be computed. In the most efficient case, where both n1 and n2 are unit normal vectors and the constant P1(0) is pre-stored, the total operations = 17 adds + 22 multiplies.
(C) 3 Plane Intersect Point. Another method selects a third plane P3 with an implicit equation n3 · P = 0 where n3 = n1 x n2 and d3 = 0 (meaning it passes through the origin). This always works since: (1) L is perpendicular to P3 and thus intersects it, and (2) the vectors n1, n2, and n3 are linearly independent. Thus the planes P1, P2 and P3 intersect in a unique point P0 which must be on L. Using the formula for the intersection of 3 planes (see the next section), where d3 = 0 for P3, we get:
Eqn_2planes-C-P0
and the parametric equation of L is:
Ps=2planes-line1
The number of operations for this solution = 11 adds + 23 multiplies.
The bottom line is that the most efficient method is the direct solution (A) that uses only 5 adds + 13 multiplies to compute the equation of the intersection line.
Intersection of 3 Planes
In 3D, three planes P1, P2 and P3 can intersect (or not) in the following ways:
Geometric Relation
Intersection
Algebraic Condition
All 3 planes are parallel
nj-x-nk=0 for all j.ne.k
They coincide
A plane
n1V1=n1V2=n1V3
They are disjoint
None
n1V1.ne.n1V2.ne.n1V3.ne.n1V1
Only two planes are parallel, and
the 3rd plane cuts each in a line
[Note: the 2 parallel planes may coincide]
2 parallel lines
[planes coincide
=> 1 line]
Only one nj-x-nk=0 for j.ne.k
No two planes are parallel, so pairwise they intersect in 3 lines
nj-x-nk.ne.0 for allj.ne.k
The intersect lines are parallel
n1-dot-n2xn3=0
They all coincide
1 line
Test a point of one line with another line
They are disjoint
3 parallel lines
Same test, but fails => do not coincide
No intersect lines are parallel
A unique point
n1-dot-n2xn3.ne.0
As shown in the illustrations:
Pic_3-planes
One should first test for the most frequent case of a unique intersect point, namely that n1-dot-n2xn3.ne.0, since this excludes all the other cases. When the intersection is a unique point, it is given by the formula:
Eqn_P0=3plane_pt
which can verified by showing that this P0 satisfies the parametric equations for all planes P1, P2 and P3.
However, there can be a problem with the robustness of this computation when the denominator n1-dot-n2xn3 is very small. In that case, it would be best to get a robust line of intersection for two of the planes, and then compute the point where this line intersects the third plane.
Implementations
Here are some sample "C++" implementations of these algorithms.
// Copyright 2001 softSurfer, 2012 Dan Sunday
// This code may be freely used and modified for any purpose
// providing that this copyright notice is included with it.
// SoftSurfer makes no warranty for this code, and cannot be held
// liable for any real or imagined damage resulting from its use.
// Users of this code must verify correctness for their application.
// Assume that classes are already given for the objects:
// 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)
// Line and Ray and Segment with defining points {Point P0, P1;}
// (a Line is infinite, Rays and Segments start at P0)
// (a Ray extends beyond P1, but a Segment ends at P1)
// Plane with a point and a normal {Point V0; Vector n;}
//===================================================================
#define SMALL_NUM 0.00000001 // anything that avoids division overflow
// dot product (3D) which allows vector operations in arguments
#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) // perp product (2D)
// intersect2D_2Segments(): find the 2D intersection of 2 finite segments
// Input: two finite segments S1 and S2
// Output: *I0 = intersect point (when it exists)
// *I1 = endpoint of intersect segment [I0,I1] (when it exists)
// Return: 0=disjoint (no intersect)
// 1=intersect in unique point I0
// 2=overlap in segment from I0 to I1
int
intersect2D_2Segments( 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);
// test if they are parallel (includes either being a point)
if (fabs(D) < SMALL_NUM) { // S1 and S2 are parallel
if (perp(u,w) != 0 || perp(v,w) != 0) {
return 0; // they are NOT collinear
}
// 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) { // both segments are points
if (S1.P0 != S2.P0) // they are distinct points
return 0;
*I0 = S1.P0; // they are the same point
return 1;
}
if (du==0) { // S1 is a single point
if (inSegment(S1.P0, S2) == 0) // but is not in S2
return 0;
*I0 = S1.P0;
return 1;
}
if (dv==0) { // S2 a single point
if (inSegment(S2.P0, S1) == 0) // but is not in S1
return 0;
*I0 = S2.P0;
return 1;
}
// they are collinear segments - get overlap (or not)
float t0, t1; // endpoints of S1 in eqn for S2
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) { // must have t0 smaller than t1
float t=t0; t0=t1; t1=t; // swap if not
}
if (t0 > 1 || t1 < 0) {
return 0; // NO overlap
}
t0 = t0<0? 0 : t0; // clip to min 0
t1 = t1>1? 1 : t1; // clip to max 1
if (t0 == t1) { // intersect is a point
*I0 = S2.P0 + t0 * v;
return 1;
}
// they overlap in a valid subsegment
*I0 = S2.P0 + t0 * v;
*I1 = S2.P0 + t1 * v;
return 2;
}
// the segments are skew and may intersect in a point
// get the intersect parameter for S1
float sI = perp(v,w) / D;
if (sI < 0 || sI > 1) // no intersect with S1
return 0;
// get the intersect parameter for S2
float tI = perp(u,w) / D;
if (tI < 0 || tI > 1) // no intersect with S2
return 0;
*I0 = S1.P0 + sI * u; // compute S1 intersect point
return 1;
}
//===================================================================
// inSegment(): determine if a point is inside a segment
// Input: a point P, and a collinear segment S
// Return: 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(): find the 3D intersection of a segment and a plane
// 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) { // segment is parallel to plane
if (N == 0) // segment lies in plane
return 2;
else
return 0; // no intersection
}
// they are not parallel
// compute intersect param
float sI = N / D;
if (sI < 0 || sI > 1)
return 0; // no intersection
*I = S.P0 + sI * u; // compute segment intersect point
return 1;
}
//===================================================================
// intersect3D_2Planes(): find the 3D intersection of two planes
// 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; // cross product
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 and Pn2 are near parallel
// test if disjoint or coincide
Vector v = Pn2.V0 - Pn1.V0;
if (dot(Pn1.n, v) == 0) // Pn2.V0 lies in Pn1
return 1; // Pn1 and Pn2 coincide
else
return 0; // Pn1 and Pn2 are disjoint
}
// Pn1 and Pn2 intersect in a line
// first determine max abs coordinate of cross product
int maxc; // max coordinate
if (ax > ay) {
if (ax > az)
maxc = 1;
else maxc = 3;
}
else {
if (ay > az)
maxc = 2;
else maxc = 3;
}
// next, to get a point on the intersect line
// zero the max coord, and solve for the other two
Point iP; // intersect point
float d1, d2; // the constants in the 2 plane equations
d1 = -dot(Pn1.n, Pn1.V0); // note: could be pre-stored with plane
d2 = -dot(Pn2.n, Pn2.V0); // ditto
switch (maxc) { // select max coordinate
case 1: // intersect with 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: // intersect with 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: // intersect with 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;
}
//===================================================================
References
James Foley, Andries van Dam, Steven Feiner & John Hughes, "Clipping Lines" in Computer Graphics (3rd Edition) (2013)
Joseph O'Rourke, "Search and Intersection" in Computational Geometry in C (2nd Edition) (1998)