The Complexity Elephant

Thoughts about software engineering, and game development.

Topics

1. Splitting Polygons, part 1 : knowing where to cut

1.1. The Problem

Splitting a polygon along a straight line is a sub-problem of polygon clipping. Many known algorithms already do polygon clipping, starting from Sutherland-Hodgeman, to Vatti’s algorithm, Weiler-Atherton, or Greiner-Horman.

These algorithms take two polygons as input, and are more or less supposed to compute their intersection (which might not always be a polygon).

For instance, the Sutherland-Hodgeman algorithm may produce a degenerate polygon if the intersection consists of multiple connected components. Similarly, the Greiner-Hormann algorithm is not robust to degeneracies.

In this context, we focus on a more specific problem: partitioning a polygon along a straight line. This can be conceptualized as clipping the polygon against two complementary half-planes defined by the line.

Throughout this discussion, we will refer to the polygon to be partitioned and the partitioning line. A key requirement here is the robust handling of degenerate cases, such as when a vertex or an edge of the polygon lies exactly on the line. Also, we’ll try to stay strictly below quadratic complexity.

Some algorithms sidestep these degenerate cases by suggesting minor perturbations to the input vertices, so the input isn’t degenerated anymore. This approach merely avoids the issue rather than addressing it. As we will demonstrate, properly handling degenerate cases constitutes the majority of the challenge in this problem.

1.2. A Naive Approach

Let’s start with a simplistic approach.

The partitioning line is represented as a normal vector and a distance to the origin. This splits the space into two half-spaces, the "positive" side and the "negative" side. Testing for which side we’re on is trivial, as it’s a simple scalar product.

aDist and bDist represent the algebraic (signed) distance to the partitioning line.

Our simplistic algorithm consists in the following : for each edge, check for an intersection with the partitioning line, and split it if needed.

#include <vector>

struct Vec2f
{
  float x, y;
};

float dotProduct(Vec2f a, Vec2f b) { return a.x * b.x + a.y * b.y; }
Vec2f operator + (Vec2f a, Vec2f b) { return { a.x + b.x, a.y + b.y }; }
Vec2f operator - (Vec2f a, Vec2f b) { return { a.x - b.x, a.y - b.y }; }
Vec2f operator * (float f, Vec2f a) { return { a.x * f, a.y * f }; }

struct Line
{
  Vec2f normal;
  float dist;
};

// only works for simplistic cases!
void splitPolygon(const std::vector<Vec2f>& polygon, Line line, std::vector<Vec2f>& outPositive, std::vector<Vec2f>& outNegative)
{
  std::vector<Vec2f> result;

  const int N = polygon.size();

  for(int i = 0; i < N; ++i)
  {
    Vec2f a = polygon[(i + 0) % N];
    Vec2f b = polygon[(i + 1) % N];

    float aDist = dotProduct(a, line.normal) - line.dist;
    float bDist = dotProduct(b, line.normal) - line.dist;

    if(aDist < 0 && bDist < 0)
      outNegative.push_back(a);
    else if(aDist >= 0 && bDist >= 0)
      outPositive.push_back(a);
    else
    {
      if(aDist < 0)
        outNegative.push_back(a);
      else
        outPositive.push_back(a);

      const float ratio = aDist / (aDist - bDist);
      const Vec2f intersection = a + ratio * (b - a);

      outPositive.push_back(intersection);
      outNegative.push_back(intersection);
    }
  }
}

void dumpToSVG(const std::vector<Vec2f>& input, const std::vector<Vec2f>& pos, const std::vector<Vec2f>& neg, Line line, const char* filename);

int main()
{
  const Line line = { { 0.71, 0.71 }, 5 };
  std::vector<Vec2f> polygon{ { 0, 0 }, { 10, 0 }, { 10, 10 }, { 0, 10 } };
  std::vector<Vec2f> pos, neg;
  splitPolygon(polygon, line, pos, neg);

  dumpToSVG(polygon, pos, neg, line, "output.svg");
  return 0;
}

Let’s try to split a square along a slanted line. Here’s the result.

Not too bad, this algorithm can indeed split some polygons. So, what’s the catch?

As hinted above, a function which can only return one polygon per side … is doomed by design.

Let’s say we’re trying to split a C-shaped polygon along a centered vertical line. The left part of the line consists of one connected component, but the right part will need two. Here’s what the correct output looks like:

However, if we try to apply our naive function to this C-shaped polygon, we will get two polygons. The left one will be correct, but the right one will be degenerated:

Indeed, there is a zero-width "bridge" along the partitioning line, linking what should have been two connected components. Depending on your application, this might suffice ( this is the kind of degenerated polygon that Sutherland-Hodgeman might produce).

However, there is a more serious limitation: when computing the intersection point, ratio is computed using the normal distance between two vertices … as a divisor.

If one of the edges happens to be almost parallel to the partitioning line, we’re in trouble. The divisor will then be extremely close to zero, and the result will contain a non-negligible amount of numerical noise. ( In practice, this leads to intersection points at bogus positions ). This is a problem, because splitting concave polygons along one of their edges is a rather classic operation (think: BSP construction).

Moreover, if one the input vertices happens to be exactly on the partitioning line, we’re also in trouble. As each vertex gets associated to one side, a side change will eventually be detected, but the newly inserted vertex will overlap the existing input vertex, creating a zero-length edge. Worse, this new vertex will be randomly assigned to the positive polygon, or to the negative polygon, once again, depending on numerical noise. Which means scaling or translating our input might change the vertex repartition in the result. This is a problem!

We’re not in the eighties anymore, and can’t hope that these situations won’t occur, or can be dismissed. We need to make special cases for vertices on the line, and also, for edges on the line.

1.3. A Better Approach

1.3.1. Classifying vertices

First, we’re going to define a minimum absolute distance to the partitioning line. Let’s this distance epsilon. All vertices below this minimum distance will be considered to be on the line. Which means now, a vertex is either on the negative side, on the positive side, or on the line. This shields us from intempestive side changes triggered by numerical noise: now, at worst, when subject to noise, a vertex will simply alternate between "on the line" and "on one definite side".

For practical and efficiency reasons, we will classify each vertex once, store the result.

  struct PointWithSide
  {
    Vec2 pos;
    int side = 0;
  };

  // classify points
  for(Vec2 pos : inputPolygon)
  {
    const float epsilon = 0.001;
    const float x = dotProduct(pos, line.normal) - line.dist;

    int side;
    if(x > +epsilon)
      side = +1;
    else if(x < -epsilon)
      side = -1;
    else
      side = 0;

    polygon.push_back({pos, side});
  }

We’re in a slightly better situation, now, each vertex now has a tri-valuated (but fixed) side.

1.3.2. Inserting the intersection points

We also want to insert the intersection points once and for all. Doing this will come with two benefits:

the list of vertices will then be completely fixed, which mean we can identify vertices by index, in a stable way.
an input vertex exactly on the line will not be a special case anymore. If one input vertex happens to be exactly on the line, it will be handled in exactly the same way that we handle the newly inserted intersection vertices.

Vec2f computeIntersection(Vec2f a, Vec2f b, Line line)
{
  const float q0 = dotProduct(a, line.normal);
  const float q1 = dotProduct(b, line.normal);
  const float ratio = (line.dist - q0) / (q1 - q0);

  return a + ratio * (b - a);
}

  // insert the intersection points
  {
    std::vector<PointWithSide> newPolygon;
    const int N = polygon.size();
    for(int curr = 0; curr < N; ++curr)
    {
      // insert the current point
      newPolygon.push_back(polygon[curr]);

      // insert the intersection point, if needed
      const int next = (curr + 1) % N;

      // different sides?
      if(polygon[curr].side * polygon[next].side == -1)
        newPolygon.push_back({computeIntersection(p0, p1, line), 0});
    }

    // replace 'polygon' with the one with all the intersections
    polygon = std::move(newPolygon);
  }

For each vertex, we compare its side with the side of the next vertex. If both vertices are on different sides, the product of their sides will be -1, and an intersection vertex will be inserted. If one of both vertices is a zero-sided vertex, the product will be 0, and no intersection vertex needs to be inserted.

Here, the computeIntersection function uses normal distance between vertices as a divisor, as before. However we don’t have the divisor problem anymore, thanks to the introduction of epsilon: If both vertices have different non-zero sides, in the worst case, the normal distance between them is at minimum twice epsilon : so we don’t worry about dividing by a noisy zero.

We replace the original polygon with our "completed" polygon. At this point, things are looking pretty good.

1.3.3. Along the partitioning line

Now, we need to express the fact that the partitioning line is indeed, well, a line. In other words, that we can go from one intersection point to the next, by following the line. We do this by sorting these intersection points by their coordinate along the line. We define a tangent vector, by rotating CCW of a quarter of a turn the partitioning line’s normal vector. Then, we compute the dot product of each intersection point with this tangent vector, and thus we can sort the intersection points.

    // sort the intersection points along the splitting line
    struct IntersectionPoint
    {
      int index;
      bool isEntryPoint = false;
    };

    std::vector<IntersectionPoint> intersectionPoints;

    for(int i = 0; i < (int)polygon.size(); ++i)
      if(polygon[i].side == 0)
        intersectionPoints.push_back({i});

    {
      const Vec2 tangent = rotateLeft(input.normal);

      auto byAbscissa = [&](const IntersectionPoint& ip0, const IntersectionPoint& ip1)
      {
        auto p0 = polygon[ip0.index].pos;
        auto p1 = polygon[ip1.index].pos;
        return dotProduct(p0, tangent) < dotProduct(p1, tangent);
      };

      std::sort(intersectionPoints.begin(), intersectionPoints.end(), byAbscissa);
    }

If we follow the partitioning line, starting from minus infinity, going to plus infinity, we’re going to cross the polygon. At minus infinity, we’re outside the polygon. Now, things get messier. What happens when we reach the first intersection point? Does it mean we enter the polygon?

Well, not always. In most cases, yes:

There are, though, tricky situations:

Let’s step back, and ask ourselves what we’re trying to compute here. We want to compute the diagonals along which to cut the original polyon.

Note that these "cut" diagonals are always pairs of consecutive intersection points.

Why consecutive? Well, a cut diagonal composed of non-consecutive intersection points would imply that at least one intersection point would lie on the interior of the cutting segment. Yet, as cutting segments are meant to become facets of the polygons resulting from the split … one of those polygons would be self-intersecting, which would mean a degenerated output.

This means that our only remaining goal here is now simple : find the first point of each diagonal. We will call these "first points" entry points, as this is where a cut "enters" the polygon to cut.

1.3.4. Finding the entry points

This is where we need to consider, at each intersection point, how the polygon crosses the partitioning line. Three points are involved here: the intersection point itself, and the previous and next points along the polygon. Note that theses points might also be intersection points. We will refer to these points as "prev", "curr", and "next".

The side of each point can be -1, 0, or +1. Since "curr" is an intersection point, its side is always zero. This leaves us with nine possible configurations to analyze:

prev	curr	next
-1	0	-1
-1	0	0
-1	0	+1
0	0	-1
0	0	0
0	0	+1
+1	0	-1
+1	0	0
+1	0	+1

However, we’re still missing some information here.

Consider the "-1,0,+1" configuration, the simplest crossing configuration. Is the partitioning line entering the polygon here? In otherwords, is "curr" an entry point?

Now, we need to clarify an aspect that we have intentionally left ambiguous for the moment: the winding of the input polygon. To discriminate entry points, we need to formally determine, given an edge of the polygon, which side is considered the interior and which side is the exterior. Without loss of generality, we shall henceforth posit that our input polygon is wound counter-clockwise. This convention implies that if one were to traverse an edge from one vertex to the next, the interior of the polygon would lie to its left.

Now, let’s go back to our "-1,0,+1" configuration. If we were to traverse in order, prev,curr,next, the partitioning line would cross our path from our right to our left, which also is the interior of the polygon. So, in this configuration, "curr" is an entry point. (Similarly, the "+1,0,-1" configuration would make "curr" an exit point, but we’re not interested in exit points).

Do we have all the information we need to discriminate entry points? Not yet.

Consider the two below cases, both corresponding to the "1, 0, 1" configuration.

On the left, the partitioning line lies at the exterior of the polygon, both making contact at one point. We don’t want to create a cut here, so "curr" shouldn’t be an entry point.

On the right, the partitioning line is locally inside the polygon. We definitely want "curr" to be the beginning of a "cut" diagonal here. So "curr" should be an entry point.

The difference between the two situations is the turning direction. In the left situation, prev-curr-next is a left turn, in the right situation prev-curr-next is a right turn.

Now we have all the required information, we only need to list the 9 * 2 possible cases.

The cases where "curr" is an entry point are:

-1 0 -1 turn right
0 0 -1 turn right
-1 0 +1 turn left
-1 0 +1 turn right
0 0 +1 turn right
+1 0 +1 turn left

Here, we use a determinant to determine the direction of the turn, and a neat C++ trick to write a switch statement on this cases.

    {
      static auto encode = [](int prevSide, int, int nextSide) { return prevSide * 8 + nextSide; };

      for(auto& ip : intersectionPoints)
      {
        const auto curr = ip.index;
        const auto prev = (ip.index - 1 + N) % N;
        const auto next = (ip.index + 1 + N) % N;

        auto prevPos = polygon[prev].pos;
        auto currPos = polygon[curr].pos;
        auto nextPos = polygon[next].pos;

        auto prevSide = polygon[prev].side;
        auto currSide = polygon[curr].side;
        auto nextSide = polygon[next].side;

        auto prevEdge = currPos - prevPos;
        auto nextEdge = nextPos - currPos;

        assert(currSide == 0); // we're on an intersection point

        const int turn = det2d(prevEdge, nextEdge) >= 0 ? +1 : -1;

        switch(encode(prevSide, currSide, nextSide))
        {
        case encode(-1, 0, -1):
          if(turn == -1)
            ip.isEntryPoint = true;
          break;
        case encode(0, 0, -1):
          if(turn == -1)
            ip.isEntryPoint = true;
          break;
        case encode(-1, 0, +1):
          ip.isEntryPoint = true;
          break;
        case encode(0, 0, +1):
          if(turn == -1)
            ip.isEntryPoint = true;
          break;
        case encode(+1, 0, +1):
          if(turn == +1)
            ip.isEntryPoint = true;
          break;
        }
      }
    }

1.3.5. Computing the list of "cut" diagonals

We’re almost done ; we now need to pair each entry point with its next intersection point, to make a list of diagonals:

    struct Diagonal
    {
      int a, b;
    };

    std::vector<Diagonal> diagonals;

    for(int i = 0; i + 1 < (int)intersectionPoints.size(); ++i)
    {
      if(intersectionPoints[i].isEntryPoint)
      {
        Diagonal diag;
        diag.a = intersectionPoints[i + 0].index;
        diag.b = intersectionPoints[i + 1].index;
        diagonals.push_back(diag);
      }
    }

1.3.6. To be continued …

And this is it, we now know were cut. Most of the difficult part is behind us. However, we still need to do the actual cutting, and some subtetlies are involved. This will be the subject for the second part of this article. Thanks for reading this far, comments are welcome!

The Complexity Elephant

Topics

Software Development

Geometry

Procedural content generation

Software Security

Meta