Surface triangulation WIP

24 Jun 2014

It’s been a while since the last post, so I thought this would be a good time for an update, even though not much worth showing got done in the mean time. With the World Cup soccer underway, I had to re-prioritize a few activities, which meant work on 2K14: The Game fell a little below the ‘threshold for progress’ over the last weeks ;-)

In the previous post I mentioned temporarily suspending development on the core engine, the basic functions of which are mostly in place, to start work on graphics rendering. The first thing I’d like to see working is rendering the planet surface the way it appears in the original game, as a solid, interlaced mass. This requires creating a polygon representing the planet mass (everything below the surface), triangulating it, and writing a GPU shader that reproduces the interlacing effect.

The triangulation problem

When setting up a planet, the shape of the planet surface is represented by a left-to-right, unconnected line segment. To render the planet mass, we need to create a closed polygon representing it, which is easy enough. The planet surface is bounded on the left and right, and the bottom of the polygon can be chosen more or less arbitrarily, as long as the resulting polygon is ‘deep’ enough that the empty space below it will never become visible. The hard part is rendering this polygon, which is almost by construction non-convex (a convex surface would make for a very boring planet). The thing is, graphics cards and the API’s used to program them don’t like non-convex polygons. In fact, they don’t like anything if it’s not a triangle. Triangulating a convex polygon is trivial, but as it turns out, triangulating a non-convex polygon isn’t. In fact, it’s surprisingly hard if you want to do it efficiently without resorting to naive methods that produce many more polygons than strictly necessary, or create very narrow/long triangles (which, without going into details here, are bad for rendering).

The figure above illustrates what we need: a function that creates a non-convex polygon representing the planet surface, and a function that efficiently triangulates this polygon. It’s debatable whether efficiency really matters all that much for the simple surface geometry expected for typical planets, probably a simple naive algorithm (for example trapezoidal decomposition of the polygon) would work equally well for our application, but I decided to go for The Right Solution™ anyway. After some reading up I found enough papers and websites pointing to what seems to be the most commonly used algorithm for general-case non-convex polygon triangulation, and it seemed a fun algorithm to implement.

The triangulation algorithm

For years, efficient polygon triangulation was a very popular topic in algorithm research. A naive exponential algorithm is trivial (keep adding non-intersecting edges inside the untriangulated area until only triangles remain), which basically amounts to brute-forcing the solution. Practical as well as ideological considerations rule out this algorithm ;-). An improved quadratic algorithm is to search for ‘ears’ of the polygon. A polygon ear is defined by two connected outline edges that can be joined by a non-intersecting edge completely contained in the polygon interior, to form a triangle that can then be clipped from the input polygon). By construction, any polygon will always have at least 2 ears, and applying the ear cliping iteratively is guaranteed to eventually eventually leave a single triangle. The Wikipedia page on polygon triangulation has a little more information about this method.

In 1978, a paper Triangulating a Simple Polygon (BiBTeX: GJPT-78) by Garey, Johnson, Preparata, and Tarjan was published in Information Processing Letters, which describes an O(n log n) algorithm for polygon triangulation. This algorithm has been improved on afterwards, and algorithms are known that have O(n log log n) complexity. In 1991 a completely different algorithm was discovered by Chazelle, that (in theory) even has O(n) complexity, but apparently it is almost impossible to implement, and not suitable for most practical purposes because of the high constant factor involved.

Because the GJPT-78 algorithm is well documented and efficient enough for our purposes, that’s the algorithm I’ve started to implement. Unfortunately the original article is behind a pay-wall everywhere, so I can’t link to it. Several ‘algorithm survey’ papers and CS course material can be found that describe the algorithm though. I mainly used the following two: UC Santa Barbara course material (no author mentioned), and University of Maryland course material (by Dave M. Mount, relevant parts in ‘lecture 6’ of the linked PDF).

The algorithm can be categorized as a two-pass sweep line algorithm algorithm, where the first pass decomposes the polygon into monotone sub-components (parts of the polygon interior with a bottom and top edge chain that both strictly run in one direction), and the second pass triangulates these monotone sub-components. It would go far beyond the scope of this post to explain the algorithm here, so for anyone interested I would suggest reading the documents linked above.

DCEL’s

An interesting sub-problem of the polygon triangulation algorithm worth mentioning is how to represent and modify the polygon during the algorithm. Typically, polygons are defined explicitly using a set of sorted vertices forming an edge loop, but this representation is not very suitable for the triangulation algorithm. For one, it does not allow splitting a single polygon into sub-components (faces) without breaking them up into disjoint vertex/edge sets and losing the connection information between these faces. Additionally, the explicit representation does not allow fast lookup of e.g. all incident edges of vertices and edges, the adjacent faces of an edge, or the edge loop surrounding a face. The triangulation algorithm depends heavily on such queries, and its runtime complexity entirely depends on the premise that these queries can be done in constant time (finding adjacent edges/faces) or (sub-) linear time (finding the edge loop around a face).

The data structure typically used for geometric problems like this is what is most commonly referred to as a Doubly Connected Edge List (DCEL) or half-edge representation. A DCEL is a data structure that stores separate vertex, edge and face sets, and links vertices with edges and edges with faces using pointers. Edges in a DCEL are split into two half-edges running in opposite directions, one along the interior of a face, and the other (its twin edge) along the interior. By linking each half edge with both its twin, its starting vertex, and with the next edge in edge loop of the same face of the polygon, topological queries can be implemented very efficiently.

Just like the triangulation algorithm itself, a full treatment of DCEL’s is beyond the scope of this post, and countless explanations can be found online that are much better than what I would be able to come up with. Have look here, here or here if you want to know more about DCEL’s.

I’ve implemented a simple DCEL representation of my own, which supports basic operations such as creating a DCEL from a simple polygon, adding edges, and retrieving incident edges, face edge loops, etc. The implementation of the DCEL data structure is short and concise, yet it took me quite a bit of time to get right. The problem with a DCEL is that while it is great for querying, it’s a bit of a drag to keep consistent under modification. There’s a lot of invariants to maintain, and debugging the data structure is an absolute nightmare because everything is tied together using pointers. Often you need to follow multiple pointers to figure out the other side of some kind of relation between vertices, edges and faces, since almost nothing is stored explicitly in the DCEL. The number of invariants and relations that need to be maintained is probably the reason why the data structure itself is usually implemented using the absolute minimum of pointers required to represent all relevant information (ie: only store a single outgoing edge with each vertex, only store the next edge in a loop and not the previous one, etc). Eventually, after a good few hours of debugging and writing an elaborate unit test to exercise the DCEL, my implementation now appears to be clean and consistent.

Bonus protip: mixing Objective-C++ with STL containers

One stupid mistake I made while working on the sweep-line part of the triangulation algorithm ate at least 3 hours of my time, so in the hope of maybe saving someone else some time let me warn you about using STL containers in Objective-C++ code here: in case you want to use an STL container (std::set in my case) as part of the implementation of an Objective-C class (I would not recommend using STL containers in the interface, since that would mean any Objective-C file including the interface has to be changed to Objective-C++), do not create a property typed std::set (or std::vector or whatever) to access the container! Because I suffer from a mild case of OCD when programming, I use Objective-C properties instead of direct ivar references whenever possible, even for private properties. This left me with a bit of code looking somewhat like this:

// K14SweepLine private interface
@interface SweepLine()
{
  ...

  @property std::set<K14DCELHalfEdge> edges;

  ...
}

@end

...

@implementation K14SweepLine 

-(NSArray *) activeEdges
{
  NSMutableArray *edges = [NSMutableArray new];

  std::set<K14DCELHalfEdge> it = self.edges.begin();

  while (it != self.edges.end())
  {
    [edges addObject:(*it)];

    it++;
  }

  return edges;
}

@end

I could of course leave it as an exercise for the reader to figure out the problem here, but I’ll go ahead and spoil it. Every time the edges property is retrieved using self.edges, behind the scenes the std::set will be copied in its entirety. This is not just wasteful, but also wrong, and dangerous. For example inserting an element into self.edges will not actually update the underlying set, but the copy returned by the property accessor. Even worse, the iterator returned by self.edges.begin() will refer to the first element of a different copy of the edges set than the self.edges.end() iterator, meaning the while loop will never terminate.

Note that directly accessing the edge set using the (auto-generated) _edges ivar would have been perfectly valid, no copies would be made this way, and the ivar itself is private so external code would never be able to modify it anyway. Moral of the story: if you want use STL containers in your Objective-C++ classes, either don’t use Objective-C properties, or make them pointers and create the container explicitly inside your initializer. Note that if you chose the latter, you also need to explicitly delete them in your class destructor, because Objective-C++ will does not automatically call destructors of C++ ivars when an object is deallocated!

Current status

The polygon triangulation algorithm itself is still very much work-in-progress, right now I have a DCEL data structure, and a simple sweep-line class that will be used to implement the 2 passes of the triangulation algorithm. Hopefully the algorithm will be done by the time of the next update, but if I run into the same amount of stupid bugs and problems as with the DCEL, it could take a little longer, we’ll have to see ;-)

Development scoreboard

Implementing and debugging the DCEL took me about 5 hours, and work on the sweep line classes somewhere around 2 hours up to now. This makes for a total of 51 hours. I’ll update the SLOC counts later, as they don’t make a lot of sense considering the unfinished state of the current code.