Algorithm 12
[Home] [Overview] [History] [Algorithms] [Books] [Web Sites] [Gift Shop]

Intersection of a Segment with a
Convex Polygon or Polyhedron

by Dan Sunday

One of the most common computer graphics computations is the clipping of a line segment with a convex (often rectangular) window.  Several efficient algorithms have been developed for doing this computation.  One that is frequently used in implementations, referred to as "parametric line-clipping" [Foley et al, 1996], was originally developed by [Cyrus & Beck, 1978] and later improved by [Liang & Barsky, 1984].  Their algorithm is usually described as a method for clipping a segment with a 2D rectangular window.  However, it easily generalizes to 2D convex polygons and 3D convex polyhedrons, and is an efficient algorithm in both cases.  This month, we describe these generalized algorithms.  We also note that [Haines, 1994] has used the same underlying ideas to compute the intersection of a ray with a 3D polyhedron.

Throughout, let S be a line segment between endpoints P₀ and P₁.  The extended line through P₀ and P₁ is given by the parametric equation: P(t) = P₀ + t (P₁-P₀) = P₀ + tv  with v = P₁-P₀ the line direction vector.  Then, the segment S contains those points P(t) with 0<=t<=1.  This representation is valid for both 2D and 3D lines.  [Note: a more detailed discussion about line equations is in Algorithm 2 about the Distance of a Point to a Line.]

Intersect a Segment and a Convex Polygon (2D)

Let a polygon W be given by n vertices V₀,V₁,...,V_n-1 going counterclockwise (ccw) around the polygon, and let V_n=V₀.  Also let e_i be the i-th edge (line segment) from vertex V_i to V_i+1 for i=0,n-1; and ev_i = (e_i1,e_i2) = (V_i+1-V_i) be the edge vector. Then, an outward-pointing normal vector for e_i is given by: n_i = -ev_i^{^} = (e_i2,-e_i1), where "^{^}" is the 2D perp-operator of [Hill, 1994] as described in Algorithm 5 about the Intersections of Lines.

Using the technique from Algorithm 5, we first compute the intersection of the (extended) line P(t) with the extended line for a single edge e_i.  Consider the following diagram:

As indicated, intersection occurs when (P(t)-V_i) · n_i = 0 since any vector in the direction of the edge is perpendicular to the edge normal vector.  Substituting for P(t) and solving for t, we get:

at the intersection point I_i = P(t_i).  The denominator here is nonzero unless the two lines are parallel, and we treat that as a special case.  In particular, when S and e_i are parallel, if P(t) is on the outside of e_i, then the segment cannot intersect the polygon at all, and we are done.  This condition can be checked by testing if the vector from V_i to P₀ points to the outside of the edge.  This is true when (P₀-V_i) · n_i > 0, in which case processing can stop since the segment is completely outside the polygon.

Next, for each edge e_i, classify the associated t_i as entering or leaving by the criteria:

as illustrated in the following diagram.

We now observe that the line P(t) must have crossed from the outside to the inside of all the edges where it is entering before it can be inside the whole polygon.  This happens when it reaches the maximum value of all the entering t_i's.  Conversely, once P(t) has crossed from the inside to the outside across any edge where it is leaving, then it can never re-enter the polygon [NOTE: this is only true for convex polygons!].  So, once t gets larger than the minimum of all the leaving t_i's, then it is outside the polygon for good.  Put:

Then, the subsegment of the line inside the polygon will start at P(t_E) and end at P(t_L) for increasing t.  But this will be non-empty only when 0 £ t_E £ t_L £ 1.  If  t_L < t_E, the segment essentially leaves the polygon before it fully enters it, and is thus is never completely inside all at once.

Pseudo-Code: Intersect Segment with Polygon

The following pseudo-code efficiently encapsulates the logic of this algorithm.

A C++ implementation intersect2D_SegPoly() is given below.

Intersect a Segment and a Convex Polyhedron (3D)

The 2D computations and algorithm extend to 3D with very few changes.  The main difference concerns specifying the polyhedron data structure.  We let a polyhedron W be given by a list of n (planar) polygon faces F_i for i=0,n-1.  We do not assume any relationship between adjacent faces F_i and F_i+1, such as having a common edge or vertex (like adjacent edges of a polygon).  It is just an unstructured list of all the distinct faces.  However, to make the algorithm work, we need some extra information about each face, namely:

1. Each face F_i lies completely in a plane.

2. V_i = a point in the plane of face F_i (e.g.: a vertex of F_i).

3. n_i = an outward-pointing normal vector to F_i ; that is, a normal vector pointing to the side of F_i  that is exterior to the polyhedron W.

This information can be precomputed from any decent data structure for a polyhedron.

Again the 3D line segment S = P₀P₁ is given by a parametric equation P(t).  For the intersection of the extended line segment with the plane of a specific face F_i, consider the following diagram.

This is exactly the same condition as in the 2D case; namely, the intersection point occurs when (P(t)-V_i) · n_i = 0 as shown.  Again, solving this gives:

at the intersection point I_i = P(t_i), and the segment is parallel to the face when the denominator is zero.  This condition is treated exactly the same.

Further, since n_i is an outward-pointing normal vector to the face plane, we have exactly the same criteria for classifying each t_i as entering or leaving; namely:

And again we compute:

Finally, the part of the segment inside the convex polyhedron is the subsegment from P(t_E) to P(t_L) when 0 £ t_E £ t_L £ 1.  If  t_L < t_E, the segment is completely outside the polygon.

[Haines, 1994] used this same method to compute the intersection of a ray and a 3D polyhedron.  He noted that t is entering the polyhedron when the intersection is with a front face plane, and is leaving when it intersects a back face plane.  However, except for the different semantics, the dot-product classification criteria are exactly the same.  Also, for a positive-extended ray starting at t=0, t_L is not bounded above by 1, and is simply the minimum of the t_i's that are leaving.

Pseudo-Code: Intersect Segment with Polyhedron

The 3D pseudo-code is almost the same as for 2D polygons, but with the obvious modifications.

Implementation

Here is a sample "C++" implementation of this algorithm. 

// 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;}
//        operators for:
//            == to test equality
//            != to test inequality
//            Point  = Point ± Vector
//            Vector = Point - Point
//            Vector = Vector ± Vector
//            Vector = Scalar * Vector    (scalar product)
//    Segment with defining endpoints {Point P0, P1;}
//===================================================================

#define TRUE    1
#define FALSE   0
#define SMALL_NUM  0.00000001 // anything that avoids division overflow

#define dot(u,v)   ((u).x * (v).x + (u).y * (v).y)   // 2D dot product
#define perp(u,v)  ((u).x * (v).y - (u).y * (v).x)   // 2D perp product
// Note: see the April 2001 Algorithm for the "perp product"
 

// intersect2D_SegPoly():
//    Input:  S = 2D segment to intersect with the convex polygon
//            n = number of 2D points in the polygon
//            V[] = array of n+1 vertex points with V[n]=V[0]
//      Note: The polygon MUST be convex and
//                have vertices oriented counterclockwise (ccw).
//            This code does not check for and verify these conditions.
//    Output: *IS = the intersection segment (when it exists)
//    Return: FALSE = no intersection
//            TRUE  = a valid intersection segment exists
int
intersect2D_SegPoly( Segment S, Point* V, int n, Segment* IS)
{
    if (S.P0 == S.P1) {        // the segment S is a single point
        // test for inclusion of S.P0 in the polygon
        *IS = S;               // same point if inside polygon
        return cn_PnPoly( S.P0, V, n );  // March 2001 Algorithm
    }

    float  tE = 0;             // the maximum entering segment parameter
    float  tL = 1;             // the minimum leaving segment parameter
    float  t, N, D;            // intersect parameter t = N / D
    Vector dS = S.P1 - S.P0;   // the segment direction vector
    Vector e;                  // edge vector
    // Vector ne;              // edge outward normal (not explicit in code)

    for (int i=0; i < n; i++)  // process polygon edge V[i]V[i+1]
    {
        e = V[i+1] - V[i];
        N = perp(e, S.P0-V[i]);// = -dot(ne, S.P0-V[i])
        D = -perp(e, dS);      // = dot(ne, dS)
        if (fabs(D) < SMALL_NUM) { // S is nearly parallel to this edge
            if (N < 0)             // P0 is outside this edge, so
                return FALSE;      // S is outside the polygon
            else                   // S cannot cross this edge, so
                continue;          // ignore this edge
        }

        t = N / D;
        if (D < 0) {           // segment S is entering across this edge
            if (t > tE) {      // new max tE
                tE = t;
                if (tE > tL)   // S enters after leaving polygon
                    return FALSE;
            }
        }
        else {                 // segment S is leaving across this edge
            if (t < tL) {      // new min tL
                tL = t;
                if (tL < tE)   // S leaves before entering polygon
                    return FALSE;
            }
        }
    }

    // tE <= tL implies that there is a valid intersection subsegment
    IS->P0 = S.P₀ + tE * dS;   // = P(tE) = point where S enters polygon
    IS->P1 = S.P₀ + tL * dS;   // = P(tL) = point where S leaves polygon
    return TRUE;
}

References

M. Cyrus & J. Beck, "Generalized Two- and Three-Dimensional Clipping", Computers and Graphics 3(1), 23-28 (1978)

James Foley, Andries van Dam, Steven Feiner & John Hughes, Computer Graphics (2nd Edition in C), Section 3.12.4 "A Parametric Line-Clipping Algorithm" (1996)

Eric Haines, "Fast Ray-Convex Polyhedron Intersection" in Graphics Gems II (1994)

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).]

Y-D. Liang & B. Barsky, "A New Concept and Method for Line Clipping", ACM TOG 3(1), 1-22 (1984)

Copyright © 2001-2006 softSurfer.  All rights reserved.
Email comments and suggestions to feedback@softsurfer.com