Subscribe
to eAlerts
for Geom Site
Updates

CLICK HERE

by Dan Sunday

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 P₀ and P₁ can be represented as:

,

where the parameter s is a real number and is a line direction vector. Using this representation , , and when , P(s) is a point on the finite segment P₀P₁ where s is the fraction of P(s)'s distance along the 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: and , 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 and are linearly related as u = av for some real number a. For and , this means that all ratios have the value a, or that for all i. This is equivalent to the conditions that all . In 2D, with and , this is the perp product [Hill, 1994] condition that where , 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 . 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 P₀, also lies on the other line Q(t). That is, there exists a number t₀ such that: P₀ = Q(t₀) = Q₀ + t₀v. If w = (w_i) = P₀ – Q₀, then this is equivalent to the condition that w = t₀v for some t₀ which is the same as for all i. In 2D, this is another perp product condition: . If this condition holds, one has , 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 t₀ and t₁ such that P₀ = Q(t₀) and P₁ = Q(t₁). If the segment intervals [t₀,t₁] and [0,1] are disjoint, there is no intersection. Otherwise, intersect the intervals (using max and min operations) to get . Then the intersection segment is . 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(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, consider the two lines and the associated vectors in the diagram:

To determine s_I, we have the vector equality where . At the intersection, the vector is perpendicular to , and this is equivalent to the perp product condition that . Solving this equation, we get:

Note that the denominator 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 , and should only be computed once if we want to know both s_I and t_I. However, knowing either is enough to get the intersection point I = P(s_I) = Q(t_I).

Further, if one of the two lines is a finite segment (or a ray), say P₀P₁, then the intersect point is in the segment only when (or 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

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: , and the plane P  be given by a point V₀ on it and a normal vector . We first check if L is parallel to P  by testing if , 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 P₀, is contained in P , that is whether it satisfies the implicit line equation: .

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 is perpendicular to n, where . This is equivalent to the dot product condition: . Solving we get:

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

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 and a normal vector n_i, and have an implicit equation: n_i · P + d_i = 0, where P = (x, y, z). 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₁ x n₂ = 0. In software, one would test if where division by would cause overflow, and uses this as the robust condition for n₁ and n₂ to be parallel in practice. When not parallel, 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 (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 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. 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 is nonzero. So, one must first select a nonzero coordinate of u, and then set the corresponding coordinate of P₀ to 0. Further, one should choose the coordinate with the largest absolute value, as this will give the most robust computations. Suppose that , then lies on L. Solving the two equations: and , one gets:

and the parametric equation of L is:

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 (n₁ x n₂). To compute 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 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₂) – P₁(0)), 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. Another method selects a third plane P₃ with an implicit equation n₃ · P = 0 where n₃ = n₁ x 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:

and the parametric equation of L is:

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 P₁, P₂ and P₃ can intersect (or not) in the following ways:

As shown in the illustrations:

One should first test for the most frequent case of a unique intersect point, namely that , since this excludes all 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 all planes P₁, P₂ and P₃.

However, there can be a problem with the robustness of this computation when the denominator 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)

© Copyright 2012 Dan Sunday, 2001 softSurfer