clipLineWithPlanes.cpp

// Clip a line segment defined by endpoints at p1[3], p2[3] by a
// collection of n planes defined by the equation ax+by+cz+d=0.
// The parametric coordinates t1,t2 of new endpoints are returned,
// as well as the indices i1,i2 of the entry and exit planes for the line.
// Result values of -1 for i1 or i2 signify that the line was not clipped
// at the p1 end or the p2 end respectively.  The function return value
// will be 0 if the whole line segment lies outside, or 1 if the line
// segment lies fully or partially inside the planes.

int clipLineWithPlanes(const double p1[3], const double p2[3],
                       const planes[][4], int n,
                       double& t1, double& t2,
                       int& i1, int& i2)
{
  // Initialize the entry and exit planes
  i1 = -1;
  i2 = -1;

  // Initialize the entry and exit parametric coordinates
  t1 = 0.0;
  t2 = 1.0;

  for (int i = 0; i < n; i++)
    {
    // Copy planes[i] into an auto variable to help compiler optimizations
    double plane[4];
    plane[0] = planes[i][0];
    plane[1] = planes[i][1];
    plane[2] = planes[i][2];
    plane[3] = planes[i][3];

    // Evaluate the equation of the plane for the two points
    double d1 = p1[0]*plane[0] + p1[1]*plane[1] + p1[2]*plane[2] + plane[3];
    double d2 = p2[0]*plane[0] + p2[1]*plane[1] + p2[2]*plane[2] + plane[3];

    // If both distances are negative, both points are outside
    if (d1 < 0 && d2 < 0)
      {
      return 0;
      }

    // If one of the distances is negative, the line crosses the plane
    else if (d1 < 0 || d2 < 0)
      {
      // We will compute fraction "t" of the crossing between p1 & p2
      double t = 0.0;

      // This "if" avoids the cpu-intensive division operation when possible
      if (d1 != 0)
        {
        // Since d1 and d2 have different signs, we know they aren't equal
        t = d1/(d1 - d2);
        }

      // If point p1 was clipped, adjust t1
      if (d1 < 0)
        {
        if (t >= t1)
          {
          t1 = t;
          i1 = i;
          }
        }
      // else point p2 was clipped, so adjust t2
      else
        {
        if (t <= t2)
          {
          t2 = t;
          i2 = i;
          }
        }

      // If this happens, there's no line left
      if (t1 > t2)
        {
        return 0;
        }
      }
    }

  return 1;
}