Surface triangulation WIP #2

11 Aug 2014

Just a really quick update to break the silence, fill the void, show some sign of life, or whatever else you could call the lack of recent updates. Progress has stalled a little over the last few weeks, for various reasons, most of which aren’t all that interesting to mention here. That said, allow me to at least point a finger at the triangulation algorithm, which turned out to be a bit of a hurdle for progress.

At the time of writing the previous update, I had a DCEL class and a simple sweep line class in place, and the ‘only’ thing left to do was to use them to implement the two passes of the triangulation algorithm: the monotone subdivision pass, and the triangulation pass. Neither of these passes are particularly hard to implement, but as it turned out the implementation requires quite a bit of boring bookkeeping, for example to keep a list of active edges sorted bottom-to-top (which is harder than you probably think, more on that later), to classify vertices as start / end / split / merge vertices, to maintain a table of ‘helper’ vertices that link the vertices and their classifications to active edges of the sweep line, calculating interior angles between incident edges, figuring out what edges and vertices constitute the upper-, and lower-edge chain of the vertex etc. Every time I sat down to with the intent to finish the triangulation algorithm I ran into some simple sub-problem that needed to be solved before I was able to piece together the main loop of each of the two subdivision passes, resulting in what is commonly called ego-depletion in social psychology (incidentally, I’m reading a great book by Daniel Kahneman on thinking processes, titled Thinking, Fast and Slow, which is where I got this fancy term from ;-) Definitely check out the book if you are interested in how your brain directs the way you think and act).

Without spending too much time (and words) on this post, progress was slow and erratic, but right now I think I’m almost at the point were I can round up the triangulation algorithm. I’ve extended the sweep line classes to classify vertices and maintain a sorted list of active edges along the sweep operation, to calculate interior angles for the edges entering and leaving vertices, and I’ve added edge properties to indicate if they are on the upper and lower edge chain. Using all these bits and pieces I was able to implement and debug the monotone subdivision pass (which seems to work), and to implement the triangulation pass (which I still have to test and debug). I spent about 20 hours on all this, which would probably could have been only 10 if I had been able to concentrate for a little longer than 2 hours each time I turned on the computer to work on the game, but I guess we’ll never know. I hope to finish and integrate the triangulation algorithm before the next post.

A note about sorting edges

As I’ve mentioned above, the sweep line class used in both passes of the triangulation algorithm needs to maintain a list of active edges, sorted top to bottom. The ordering is needed to quickly (O(log N) time) insert and delete edges during the sweep operation, or to find the edge directly above or below a certain y coordinate at the current sweep line x coordinate. Sorting edges top to bottom may not sound like a particularly hard problem to solve, but it took me surprisingly long to sort out (pardon the pun ;-).

Because the Objective-C foundation framework does not really provide a suitable ordered set class that allows you to specify a custom sort function (at least not one that doesn’t require re-sorting after each insertion), I used a std::set with a custom comparator class. The problem I ran into was that during the sweep operation, sometimes the std::set failed to insert edges, inserted them twice, or inserted them at the wrong location. No matter how sane my comparator class looked, it did not work reliably. The issue turned out to be my failure to recognize that you can’t actually sort edges top to bottom, because unless the edges are parallel, they will intersect at some point, inverting their relative orientation. All of my initial attempts at sorting the edges used a comparator function that either compared the y-coordinates of the midpoints of edges, projected the starting point of one edge onto the infinite extension of the other edge, etc. None of these approaches yield a complete partial ordering of the edges, because such an ordering is only possible if you restrict it to the range of x-coordinates where the edges overlap. When this simple realization finally dawned on me, fixing the comparator class was straightforward. Here’s the code for it:

// Custom comparator used to order edges from top to bottom
//
// This comparator is extremely specific to the sweep line algorithm, and will only guarantee a
// complete partial ordering on edges under the assumption that edges are 1) always running left-to-right,
// 2) are always added in left-to-right order, and 3) are always removed from the set that uses this
// comparator as soon as the sweep line sweeps their rightmost vertex.
struct EdgeCompare
{
  // Comparison function.
  bool operator() (const SweepLineEdge *lhs, const SweepLineEdge *rhs) const
  {
    bool swapped = false;
   
    // Swap edges if the right-hand side edge extends the left-hand side edge on the left, to
    // reduce the number of comparison cases we need to handle.
    if (rhs.v0.x < lhs.v0.x)
    {
      std::swap(lhs, rhs);
      swapped = true;
    }
    
    // Handle the special case where (possibly after swapping), the left-hand side edge is singular. By
    // construction, the singular edge will be to the left of the starting vertex of the right-hand side
    // edge otherwise the edges would have been swapped. Calculate (oblique) projection of the singular
    // edge onto the right-hand side edge, and compare on y-coordinate.
    if (CGPointEqualToPoint(lhs.v0, lhs.v1))
    {
      float t = (lhs.v0.x - rhs.v0.x) / (rhs.v1.x - rhs.v0.x);
      
      GLKVector2 p = GLKVector2Lerp(GLKVector2Make(rhs.v0.x, rhs.v0.y), GLKVector2Make(rhs.v1.x, rhs.v1.y), t);
      
      return (swapped ? p.y < lhs.v0.y : lhs.v0.y < p.y);
    }
    
    // Handle special case where the two edges are connected left-to-right, which means they would otherwise be
    // considered equal. In theory it isn't necessary to handle this case, under the assumption the left-hand edge
    // would already have been removed from the sweep line before the right-hand edge is added (which means this case
    // will never actually occur). We explicitly check for it anyway, to prevent accidental bugs in case the sweep
    // operation implementation is changed.
    if (CGPointEqualToPoint(lhs.v1, rhs.v0))
    {
      return (swapped ? rhs.v1.y < lhs.v1.y : lhs.v1.y < rhs.v1.y);
    }
    
    // Handle the special case where the 2 edges share the same start vertex. We just compare the y-coordinates
    // of the rightmost points of both edges to decide which of the two is the lower edge.
    if (CGPointEqualToPoint(lhs.v0, rhs.v0))
    {
      return (swapped ? rhs.v1.y < lhs.v1.y : lhs.v1.y < rhs.v1.y);
    }
    
    // Left-hand side edge is not singular, and by construction extends the right-hand side on the left.
    // Calculate (oblique) projection of right-hand side start vertex onto left-hand side edge and compare
    // y-coordinates of this start vertex and its projection
    float t = (rhs.v0.x - lhs.v0.x) / (lhs.v1.x - lhs.v0.x);
    
    GLKVector2 p = GLKVector2Lerp(GLKVector2Make(lhs.v0.x, lhs.v0.y), GLKVector2Make(lhs.v1.x, lhs.v1.y), t);
    
    return (swapped ? rhs.v0.y < p.y : p.y < rhs.v0.y);
  }
};

Note that this comparator only works under the assumption that any pair of edges in the set have non-zero overlap in x, in other words: edges completely to the left of the sweep line position have to be removed from the set, before any new edges are added.

Development scoreboard

I spent about 20 hours to get the triangulation algorithm to its current state, which is ‘mostly done, but mostly untested’. I don’t expect to spend more than an hour or 3 to finish it. Total development time is now at about 71 hours. The SLOC count is 1249, which is an increase of about 400, most of which went into the supporting code for the triangulation algorithm (the DCEL and sweep line classes, and all the bookkeeping structures and enumeration types). Sometimes I start wondering whether I shouldn’t just have taken the easy way out and implement some kind of suboptimal hack such as trapezoidal decomposition instead of a ‘real’ triangualation algorithm, it would have been done long ago, with much less lines of code… I’ll probably be a lot happier with my decision when the triangulation algorithm is finished ;-)