Polygon plane partition

From CGAFaq

Overview

Partitioning a polygon with a plane is a matter of determining which side of the plane the polygon is on. This is referred to as a front/back test, and is performed by testing each point in the polygon against the plane. If all of the points lie to one side of the plane, then the entire polygon is on that side and does not need to be split. If some points lie on both sides of the plane, then the polygon is split into two or more pieces.

The basic algorithm is to loop across all the edges of the polygon and find those for which one vertex is on each side of the partition plane. The intersection points of these edges and the plane are computed, and those points are used as new vertices for the resulting pieces.

Implementation notes

Classifying a point with respect to a plane is done by passing the (x, y, z) values of the point into the plane equation, . The result of this operation is the distance from the plane to the point along the plane's normal vector. It will be positive if the point is on the side of the plane pointed to by the normal vector, negative otherwise. If the result is 0, the point is on the plane.

For those not familiar with the plane equation, The values A, B, and C are the coordinate values of the normal vector. D can be calculated by substituting a point known to be on the plane for x, y, and z.

Convex_polygons are generally easier to deal with in BSP tree construction than concave ones, because splitting them with a plane always results in exactly two convex pieces. Furthermore, the algorithm for splitting convex polygons is straightforward and robust. Splitting of concave polygons, especially self intersecting ones, is a significant problem in its own right.

Pseudo C++ code example Here is a very basic function to split a convex polygon with a plane:

 void Split_Polygon (polygon *poly, plane *part, polygon *&front,
                     polygon *&back)
{
  int count = poly->NumVertices ();
  int  out_c = 0, in_c = 0;
  point ptA, ptB, outpts[MAXPTS], inpts[MAXPTS];
  real sideA, sideB;
  ptA = poly->Vertex (count - 1);
  sideA = part->Classify_Point (ptA);
  for (short i = -1; ++i < count;)
  {
    ptB = poly->Vertex (i);
    sideB = part->Classify_Point (ptB);
    if (sideB > 0)
    {
      if (sideA < 0)
      {
        // compute the intersection point of the line
        // from point A to point B with the partition
        // plane. This is a simple ray-plane intersection.
        vector v = ptB - ptA;
        real   sect = - part->Classify_Point (ptA) / (part->Normal () | v);
        outpts[out_c++] = inpts[in_c++] = ptA + (v * sect); 
      }
      outpts[out_c++] = ptB;
    }
    else if (sideB < 0)
    {
      if (sideA > 0)
      {
        // compute the intersection point of the line
        // from point A to point B with the partition
        // plane. This is a simple ray-plane intersection.
        vector v = ptB - ptA;
        real   sect = - part->Classify_Point (ptA) / (part->Normal () | v);
        outpts[out_c++] = inpts[in_c++] = ptA + (v * sect);
      }
      inpts[in_c++] = ptB;
    }
    else
      outpts[out_c++] = inpts[in_c++] = ptB;
    ptA = ptB;
    sideA = sideB;
  }
  front = new polygon (outpts, out_c);
  back = new polygon (inpts, in_c);
}

A simple extension to this code that is good for BSP trees is to combine its functionality with the routine to classify a polygon with respect to a plane.

Note that this code is not robust, since numerical stability may cause errors in the classification of a point. The standard solution is to make the plane "thick" by use of an epsilon value.

Initial content for this page from the BSP tree FAQ, included with permission of its maintainer, Bretton Wade.