Intersection of Lines and Planes

The intersection of geometric primitives is a fundamental construct in many computer graphics and modeling applications ([Foley et al, 1996], [Mortenson, 1990], [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. Other representations are discussed in Algorithm 2 about the Distance of a Point to a Line, Ray or Segment. 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 P₀ and P₁ can be written as

where the parameter s is a real number and u=P₁-P₀ is a line direction vector. Using this representation P(0)=P₀, P(1)=P₁, and when 0<s<1, P(s) is a point on the finite segment between P₀ and P₁ where s is the fraction of P(s)'s distance along the whole P₀P₁ line segment. That is, s = d(P₀,P(s)) / d(P₀,P₁). Further, if s<0 then P(s) is outside the segment on the P₀ side, and if s>1 then P(s) is outside the segment on the P₁ side.

Let two lines be given by equations: P(s) = P₀+ s (P₁-P₀) and Q(t) = Q₀+ t (Q₁-Q₀), either or both of which could be a finite segment or a ray. These lines are parallel when and only when their directions are collinear, namely when the two vectors u=P₁-P₀ and v=Q₁-Q₀ are linearly related as u=Cv for some real number C. For u=(u_k) and v=(v_k), this means that all ratios u_k/v_k have the same value, or that u₁/v₁=u_k/v_k for all k, which is equivalent to the conditions that all u₁v_k-u_kv₁=0. In 2D, with u=(u₁,u₂) and v=(v₁,v₂), this is the perp product [Hill, 1994] condition u^{^}· v = u₁v₂-u₂v₁ = 0 where u^{^}=(-u₂,u₁) is perpendicular to u ("^{^}" is called the perp operator). This condition says that two lines in the Euclidean plane are parallel when they are both perpendicular to the same direction u^{^}. When true, the two lines are either coincident or do not intersect at all.

Coincidence is easily checked by testing if a point on one line, say P₀, also lies on the other line Q(t). That is, there exists a t₀ such that: P₀=Q(t₀)=Q₀+t₀v. If w=(w_k)=P₀-Q₀, then this is equivalent to the condition that w=t₀v for some t₀ which is the same as w₁v_k-w_kv₁=0 for all k. In 2D, this is the perp product condition: w^{^}· v = w₁v₂-w₂v₁= 0. If these conditions hold, one can solve for t₀, and the infinite lines are coincident. 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 t₀ and t₁ such that P₀=Q(t₀) and P₁=Q(t₁). If the segment intervals [t₀,t₁] and [0,1] are disjoint, then there is no intersection. Otherwise, intersect the intervals (using max and min operations) to get [r₀,r₁] = [t₀,t₁] Ç [0,1]. Then the intersection segment is Q(r₀)Q(r₁) = P₀P₁ Ç Q₀Q₁. This works in any dimension.

When the two lines are not parallel, they may intersect in a unique point. In 2D Euclidean space, they always intersect. In higher dimensions they may miss each other and not intersect. But if they do intersect, then their linear projections onto a 2D plane will also intersect. So, one can simply restrict to two coordinates (for which u or v are nonzero), compute the 2D intersection point I at P(s_I) and Q(t_I) for those two coordinates, and then test if P(s_I)=Q(t_I) for all coordinates. To compute the 2D intersection point I=P(s_I)=Q(t_I), consider the two lines and the associated vectors in the diagram:

To determine s_I, we have the vector equality P(s)-Q₀= w + suwhere w=P₀-Q₀. At the intersection point, the vector P(s)-Q₀ is perpendicular to v^{^}, and this is equivalent to the perp product condition that v^{^}· (w+su) = 0. Solving this equation, we get:

Note that the denominator v^{^} · u = 0 only when the lines are parallel as previously discussed. Similarly, solving for Q(t_I), we get:

The denominators are the same up to sign, since u^{^}·v = -v^{^}·u, and should only be computed once if we want to know both s_I and t_I.

If one of the two lines is a finite segment (or a ray), say P₀P₁, then the intersect point is inside the segment only when 0<=s_I<=1 (or s_I>=0 for a ray). If both lines are segments, then both solution parameters, s_I and t_I, 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

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) = P₀+ s(P₁-P₀) = P₀+ su, and the plane p be given by the point V₀ on it and it's normal vector n. We first check if L is parallel to p by testing if n·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 specific point of L, say P₀, is contained in p, that is whether it satisfies the implicit line equation: n·(P-V₀) = 0.

If the line and plane are not parallel, then L and p intersect in a unique point P(s_I) which is computed using a method similar to the one for the intersection of two lines in 2D. Consider the diagram:

At the intersect point, the vector P(s)-V₀= w+su,where w=P₀-V₀, is perpendicular to n. This is equivalent to the dot product condition: n · (w+su) = 0. Solving we get:

If the line L is a finite segment from P₀ to P₁, then one just has to check that 0<=s_I<=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>=0.

Intersection of 2 Planes

In 3D, two planes p₁ and p₂ are either parallel or they intersect in a single straight line L. Let p_i (i=1,2) be given by a point V_i plus a normal vector n_i, and have an implicit equation: n_i· P + d_i= 0. The planes p₁ and p₂ are parallel whenever their normal vectors n₁ and n₂ are parallel, and this is equivalent to the condition that: n₁×n₂= 0. In software, one usually tests if |n₁×n₂| < d where division by d would cause overflow, and uses this as the robust condition for n₁ and n₂ being parallel. When not parallel, u = n₁×n₂³ d > 0 is a direction vector for the intersection line L since u is perpendicular to both n₁ and n₂, 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.

After computing n₁×n₂ (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 P₀=(x₀,y₀,z₀) that lies in both planes. We can do this by finding a common solution of the implicit equations for p₁ and p₂. But there are only two equations in the 3 unknowns since the point P₀ can lie anywhere on the 1-dimensional line L which is another degree of freedom. So we need another constraint to solve for a specific P₀. There are a number of ways this could be done:

(A) Direct Linear Equation. One could set one coordinate to zero, say z₀=0, and then solve for the other two. But this will only work when L intersects the plane z₀=0. This is will be true when the z-coordinate u_z of u =(u_x,u_y,u_z) is nonzero. So, one must first select a nonzero coordinate of u, and then set that coordinate of P₀ = 0. One should choose the coordinate with the largest absolute value, as this will give the most robust computations. Suppose that u_z¹0, then P₀=(x₀,y₀,0) lies on L. Solving the two equations: a₁x₀+b₁y₀+d₁=0 and a₂x₀+b₂y₀+d₂=0, we get:

The denominator here is equal to the non-zero 3rd coordinate of u. So, ignoring the test for a nonzero coordinate, and equating division with multiplication, the total operations for this solution = 5 adds + 12 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, this will be a point on the line of intersection for the two planes, and the parametric equation of L is: P(s) = I + s (n₁×n₂). To compute n₁×n₂ and the intersection point (given the line) has a 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 n₂ on p₁ defines a line that intersects p₂ in the sought for point P₀ on L. More specifically, project the two points 0=(0,0,0) and n₂=(nx₂,ny₂,nz₂) to p₁(0) and p₁(n₂) respectively. Then the projected line in p₁ is L₁: Q(t) = p₁(0) + t p₁(n₂), and intersection of it with p₂ can be computed. In the most efficient case, where both n₁ and n₂ are unit normal vectors and the constant p₁(0) is pre-stored, the total operations = 17 adds + 22 multiplies.

(C) 3 Plane Intersect Point. One can select a third plane p₃ with an implicit equation n₃· P = 0 where the normal vector n₃= n₁×n₂ and d₃=0 (meaning it passes through the origin). This always works since: (1) L is perpendicular to p₃ and thus intersects it, and (2) the vectors n₁, n₂, and n₃ are linearly independent. Thus the planes p₁, p₂ and p₃ intersect in a unique point P₀ which must be on L. Using the formula for the intersection of 3 planes (see the next section), where d₃=0 for p₃, we get:

This formula can be simplified when n₁ and n₂ are unit normal vectors, and n₁· n₂ = cos q where q is the angle between these vectors. Then, n₁· n₁ = n₂· n₂ = 1, and |n₁×n₂|² = sin² q = 1 - cos² q = 1 - (n₁· n₂)² for the denominator. Finally, we use the identity for right-association of the cross-product; namely, for any three (3D) vectors a, b, and c, it is always true that: a×(b×c) = (a· c)b - (a· b)c. Applying this to the numerator, and simplifying, we get: (d₂n₁-d₁n₂)×(n₁×n₂) = (d₂ cos q - d₁) n₁ + (d₁ cos q - d₂) n₂ . Thus,

This solution expresses the point P₀ on L as a linear combination of the vectors n₁ and n₂ that are a basis spanning the plane p₃ through the origin. The number of operations = 11 adds + 20 multiplies.

But the bottom line is that the most efficient method is the direct solution (A) that only uses 5 adds + 12 multiplies to compute the equation of the intersection line.

[Footnote: I would like to acknowledge a fruitful email exchange with Nick Polyayev in July 2001 discussing alternative methods for determining the intersection line for 2 planes. This resulted in a greatly improved presentation of this section.]

Intersection of 3 Planes

In 3D, three planes p₁, p₂ and p₃ can intersect (or not) in the following ways:


Geometric Relation	Intersection	Algebraic Condition
All 3 planes are parallel		n_j×n_k=0 for all j,k=1,3
They coincide	A plane	n₁·V₁=n₁·V₂=n₁·V₃
They are disjoint	None	n₁·V₁¹n₁·V₂¹n₁·V₃¹n₁·V₁

Just two planes are parallel, and the 3rd plane cuts each in a line [Note: the 2 parallel planes may coincide]	2 parallel lines [coincide => 1 line]	Only one n_j×n_k=0 for j¹k

No two planes are parallel, so pairwise they intersect in 3 lines		n_j×n_k¹0 for all j¹k
The intersection lines are parallel		n₁·(n₂×n₃)=0
They coincide	1 line	test a point from one line
They are disjoint	3 parallel lines
No intersection lines are parallel	A unique point	n₁·(n₂×n₃)¹0

One should test for the most frequent case of a unique intersect point first, namely that n₁·(n₂×n₃)¹0, since this excludes the other cases. When the intersection is a unique point, it is given by the formula:

which can verified by showing that this P₀ satisfies the parametric equations for p₁, p₂ and p₃.

However, there can be a problem with the robustness of this computation when the denominator n₁·(n₂×n₃) is very small. In this 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

// Copyright 2001, softSurfer (www.softsurfer.com)
// 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(): the intersection of 2 finite 2D 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_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);

    // 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(): intersect 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(): the 3D intersect 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

Francis Hill, "The Pleasures of 'Perp Dot' Products" in Graphics Gems IV (1994)
[Note: the first critical definition has a typo, and should be: a^{^}= (-a_y, a_x).]

Intersections of Lines, Segments and Planes (2D and 3D)

Line and Segment Intersections

Plane Intersections

Line-Plane Intersection

Intersection of 2 Planes

Intersection of 3 Planes

Implementations

References

Intersections of Lines, Segments
and Planes (2D and 3D)