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.