Determining whether a line segment intersects a 3 vertex facet Written by Paul Bourke February 1997 The following will find the intersection point (if it exists) between a line segment and a planar 3 vertex facet. The mathematics and solution can also be used to find the intersection between a plane and line, a simpler problem. The intersection between more complex polygons can be found by first triangulating them into multiple 3 vertex facets.

Source code will be provided at the end, it illustrates the solution more than being written for efficiency. The labeling and naming conventions for the line segment and the facet are shown in the following diagram

The procedure will be implemented given the line segment defined by its two end points and the facet bounded by its three vertices.

• Check the line and plane are not parallel
  
• Find the intersection of the line, on which the given line segment lies, with the plane containing the facet
  
• Check that the intersection point lies along the line segment
  
• Check that the intersection point lies within the facet
  

The intersection point P is found by substituting the equation for the line P = P₁ + mu (P₂ - P₁) into the equation for the plane Ax + By + Cz + D = 0.

Note that the values of A,B,C are the components of the normal to the plane which can be found by taking the cross product of any two normalised edge vectors, for example

Then D is found by substituting one vertex into the equation for the plane for example

This gives an expression for mu from which the point of intersection P can be found using the equation of the line.

If the denominator above is 0 then the line is parallel to the plane and they don't intersect. For the intersection point to lie on the line segment, mu must be between 0 and 1.

Lastly, we need to check whether or not the intersection point lies within the planar facet bounded by P_a, P_b, P_c

The method used here relies on the fact that the sum of the internal angles of a point on the interior of a triangle is 2pi, points outside the triangular facet will have lower angle sums.

If we form the unit vectors P_a1, P_a2, P_a3 as follows (P is the point being tested to see if it is in the interior)

P_a1 = (P_a - P) / |(P_a - P)|
P_a2 = (P_b - P) / |(P_b - P)|
P_a3 = (P_c - P) / |(P_c - P)|

the angles are

a1 = acos(P_a1 dot P_a2)
a2 = acos(P_a2 dot P_a3)
a3 = acos(P_a3 dot P_a1)

Source code

/*
   Determine whether or not the line segment p1,p2
   Intersects the 3 vertex facet bounded by pa,pb,pc
   Return true/false and the intersection point p

   The equation of the line is p = p1 + mu (p2 - p1)
   The equation of the plane is a x + b y + c z + d = 0
                                n.x x + n.y y + n.z z + d = 0
*/
int LineFacet(p1,p2,pa,pb,pc,p)
XYZ p1,p2,pa,pb,pc,*p;
{
   double d;
   double a1,a2,a3;
   double total,denom,mu;
   XYZ n,pa1,pa2,pa3;

   /* Calculate the parameters for the plane */
   n.x = (pb.y - pa.y)*(pc.z - pa.z) - (pb.z - pa.z)*(pc.y - pa.y);
   n.y = (pb.z - pa.z)*(pc.x - pa.x) - (pb.x - pa.x)*(pc.z - pa.z);
   n.z = (pb.x - pa.x)*(pc.y - pa.y) - (pb.y - pa.y)*(pc.x - pa.x);
   Normalise(&n);
   d = - n.x * pa.x - n.y * pa.y - n.z * pa.z;

   /* Calculate the position on the line that intersects the plane */
   denom = n.x * (p2.x - p1.x) + n.y * (p2.y - p1.y) + n.z * (p2.z - p1.z);
   if (ABS(denom) < EPS)         /* Line and plane don't intersect */
      return(FALSE);
   mu = - (d + n.x * p1.x + n.y * p1.y + n.z * p1.z) / denom;
   p->x = p1.x + mu * (p2.x - p1.x);
   p->y = p1.y + mu * (p2.y - p1.y);
   p->z = p1.z + mu * (p2.z - p1.z);
   if (mu < 0 || mu > 1)   /* Intersection not along line segment */
      return(FALSE);

   /* Determine whether or not the intersection point is bounded by pa,pb,pc */
   pa1.x = pa.x - p->x;
   pa1.y = pa.y - p->y;
   pa1.z = pa.z - p->z;
   Normalise(&pa1);
   pa2.x = pb.x - p->x;
   pa2.y = pb.y - p->y;
   pa2.z = pb.z - p->z;
   Normalise(&pa2);
   pa3.x = pc.x - p->x;
   pa3.y = pc.y - p->y;
   pa3.z = pc.z - p->z;
   Normalise(&pa3);
   a1 = pa1.x*pa2.x + pa1.y*pa2.y + pa1.z*pa2.z;
   a2 = pa2.x*pa3.x + pa2.y*pa3.y + pa2.z*pa3.z;
   a3 = pa3.x*pa1.x + pa3.y*pa1.y + pa3.z*pa1.z;
   total = (acos(a1) + acos(a2) + acos(a3)) * RTOD;
   if (ABS(total - 360) > EPS)
      return(FALSE);

   return(TRUE);
}