A Healthy Fear of Threading

Continuing in the line of pithy quotes:

There are only two kinds of programmers: programmers with a healthy fear of threaded code and programmers who should fear code.

Now I'm not saying "never thread". I'm just saying "you better be getting something good for that threading, because it's driving up your development costs."

In particular, the effective execution order of threaded code can change with every run, and there is no guarantee that you have seen every combination of execution order by running your program a finite number of times.

Thus methods of checking your code quality by running your program (perhaps many times) won't detect bugs in threaded code. You may not find out until that user with one more core and a background program chewing up cycles hits an execution order that you haven't seen yet.

Instead for threaded code you have to prove logically that the execution order constraints applied (via locking, etc.) create a bounded set of execution combinations, and that each one is correct. This isn't quick or easy to do.

One way we cope with this development cost in X-Plane (where we need to use threads to fully utilize multiple cores) is to use threading design patterns with known execution limits. The most common one is a message queue, where ownership of data access flows with the message down a queue. This idiom not only guarantees serialized access to data without locks, but the implementation in C++ tends to make errors rare; if you have the message you have the pointer, and thus you have rights on the data. If you don't have the message, you have nothing to dereference.

How Does OpenGL Work?

If you are registered with Apple's developer sites, I strongly recommend the OpenGL and OpenGL ES video talks from WWDC 2010. The Apple engineers spell out in a fair amount of detail things that you had to infer previously, including:

• How state information is accumulated and then resynchronized at draw call time.
• How resources are synchronized and shared between host memory and the GPU.

The videos are in QuickTime format with subtitles, so you can play them back at 2x speed with captioning to get through the material faster.

Interval Sets With the STL

I spent some time today working on an interval set. The basic idea of an interval set is to record a set of disjoint ranges that partition a number space into a finite "included" area and an infinite "excluded" area. Or to put it more simply, [3, 6) is an interval, and [3, 6) [8, 10) is an interval set.

Googling around for this I found a few ideas based on using a sorted map, with the interval beginning as the key and the interval end as the value. My approach is different, and is closer to the original implementation of Macintosh regions: a vector of beginning and ending pairs, e.g. { 3, 6, 8, 10 }.

I'm not sure whether this approach is superior to a map-based approach; I think I'd have to code each one all the way to completion. The vector does have a few advantages:

• Compact storage, with minimal overhead.
• Reading the sorted array can usually be done in linear or log-N time.

The basic rules for the interval set are:

• Intervals are inclusive at the bottom and exclusive at the top. So the interval 3,6 includes the number 3 but excludes the number 6. (Thus given two intervals [3,6) and [6,9) the number 6 is included in exactly one of them.)
• Intervals must have non-zero length (so [3,3) is illegal).
• Intervals are finite - there is no notation to say "everything below 3" is included.

The heart of the algorithm is a "merge" operation. In a merge, the sequence of interval edges from two interval sets are traversed together (think of a merge sort) and each new smallest sub-interval is evaluated for inclusion by a boolean operation. This lets us perform a union, difference, intersection, or symmetric difference in linear time O(N+M) where N and M are the lengths of the vectors.

(The actual time is actually slightly worse because vector will need to periodically reallocate its memory during the creation of the new resulting vector. We could use a heuristic to pre-allocate some space at a loss of memory efficiency. If we used a set we could avoid memory costs, but we'd end up with O(NlogN) time to build the set anyway, and we'd pay node overhead, which is almost certainly worse than any extra on a vector of 32-bit floats or integers.)

When we search for an interval (using lower_bounds) we can tell whether we are "in" or "out" of the region by looking at whether the index of the returned region is even or odd - even regions are in the set and odd ones are outside of it.

The interval class is also heavily special cased for a number of optimizations:

• Separate operators on pair allow for the processing of a single interval (rather than a set). When we know the single interval, we can take a number of short-cuts, and we can perform "in-place editing" using log-N searches into the original interval set.
• Operations on sets can identify short cuts. For example, the intersection of two sets whose range is disjoint is always empty. (In other words, if the last value in A is less than the first value in B, intersecting A and B is an empty set.)

I haven't used the interval set class enough to profile it; real measurement will tell which of these optimizations is a win. One tricky aspect of the code is that vector is a leaky abstraction - it makes mid-vector insertion look cheap when really it is a linear operation (because all subsequent elements must be copied to their new locations).

As an example of why this might matter: consider symmetric difference (XOR) of an interval set and a single range. This operation can be computed simply by: deleting the range bounds from the set if they exist, otherwise inserting them. In other words, given the interval set [0,3) [6,9) [12,15) we can XOR this with the interval [6,8) by deleting 6 and inserting 8 - the new XOR is [0,3) [8,9) [12,15). This is a relatively fast operation: two log-N searches (for 6 and 8) and one delete and one insert.

Despite the simplicity of the algorithm, vector is going to require two mid-vector editing operations, so our average time complexity is O(N) - linear! (On average half the elements of the vector are after us, and we do two editing ops.)

For this reason, the special case of a disjoint XOR is special cased. If we XOR [-10, -8) into the above region, we can observe that -8 < 0, therefore the regions don't intersect, and -10, -8 simply needs to be pre-pended. This can be done with a single insert, and thus should run about twice as fast as a pair of individual inserts.

Importing Faces Into CGAL Arrangements

This problems occurs repeatedly in the X-Plane scenery tools code: we need to import a series of polygons into a single arrangement_2 structure, and we want to tag the faces contained by these polygons as they appear in the arrangement. This isn't entirely trivial for two reasons:

• The polygons may collide with each other and thus there may be a 1:many relationship between the original data and the final map.
• The polygons may be self-intersected or in other ways "hosed", so we need a particular strategy for handling this situation.

CGAL provides a number of built-in tools to deal with these situations. Here are 3 basic and useful building blocks:

1. If you are using the general polygon set code that is built on top of arrangements, you can simply perform unions and intersections of a large number of faces. For merging a large number of areas, this code is faster than anything else you might code, because it can do an N-way divide and conquer, where N is larger than 2.
2. For custom merging and handling of multiple polygons, you can use the overlay free function, which lets you specify how the combinations of each set of face from two maps are handled.
3. You can simply insert a set of curves into an empty arrangement and they will be "swept" together. This is a useful way to turn a messy polygon into something useful - it finds intersections, builds topology, runs quickly, and handles input no matter how degenerate.

For example, if the goal is to have any area contained by any piece of polygon as "inside" you can simply insert all polygon sides into an empty arrangement and then tag every bounded face.

Finding Polygon Internal Areas

When building arrangements out of polygons of dubious origin (or simply building an arrangement out of a large number of unrelated polygons) I use a bulk insert to "sweep" the curves into the arrangement. How do I then find the faces? Here are three techniques:

1. The contained area of a polygon can be found by simply checking whether the face is bounded or not. This is not useful though when importing multiple polygons at the same time. (When I need this technique, each polygon is individually imported into its own arrangement, then all arrangements are merged later, typically with general-polygon-set code.)

2. We can implement a "toggle" policy (e.g. each line toggles interior vs. exterior) by doing a search from the outside to the inside of the arrangement, toggling whether we are "inside" or "outside" each time we cross a halfedge. The halfedges can retain curve-based properties; typically I use a consolidated data curve so that halfedges retain every property attached to them.

   One danger: an antenna will produce incorrect results in this technique because it won't toggle the data property twice. This can be hard to work around because data from the insert is maintained per edge, not half-edge

   Unfortunately, topology of the final arrangement doesn't help us resolve this either. Imagine an antenna (in the original polygon that crashes into another polygon, thus becoming part of a real partition. Technically the original face is not split by the antenna, but the faces in the final arrangement are, so saying face()==twin()->face() doesn't tell us we have an antenna in the original.

   Two ways to work around this: don't insert antennas, or don't tag known antennas with any data. Both cases require knowing that we have an antenna ahead of time.

3. Sometimes we want to use a "winding rule" - that is, contain areas inside closed left turning contours. This is, for example, useful when calculating offset buffers and minkowski sums; the artifacts from strange shapes being offset too much turn out not to be left turning contours and get thrown out.

   To find the winding rule areas, we have to look at the direction of the curves inserted.

   The way my code does this is to look at the direction of the underlying curve vs. the direction of the half-edge and then mark the half-edge as being on the "inside" or "outside" of the winding with a dat a field on the half-edge itself. This is reconstructed after bulk insert, and then we can traverse the whole arrangement, counting windings.

   The limitation here is similar to above: if we have an antenna, the underlying curve can have only one direction, not two, and one half-edge will be incorrectly tagged. Fortunately antennas are not typically necessary to produce offset buffers.
   

It should be noted that if you insert curves incrementally (insert one curve into the arrangement) an observer of the arrangement returns all generated and overlapped half-edges, which gives you the contained bounds of the contour inserted. I use this technique when inserting a low-side-count face into a very complex arrangement, to avoid re-sweeping a huge amount of data. Overlays and bulk inserts do not produce "per curve" announcements via the observer mechanism.

CGAL: It's All About the Mantissa

In a past post I described CGAL as having no rounding errors. It does this by using number types of variable size (using dynamically allocated memory per number!) so that it never runs out of digits. (It also maintains the numerator and denominator of fractions separately to avoid problems with repeating decimals.)

The advantage of this is that geometric algorithms that rely on precise calculations never go haywire due to rounding errors. For example, when using fixed-precision math (e.g. IEEE floats) the intersection of two near-parallel lines will be calculated inaccurately - sometimes with the intersection showing up miles from the original lines. CGAL always has more precision, so it avoids this problem.

But there is one down-side: when you perform a series of intersections, the result is exact numbers whose mantissas (the number of actual digits) have grown very long. And CGAL won't blink about making them even longer as you do more calculations.

Instead CGAL will become insanely slow.

I hit this case the other day. The first piece of processing I do is to combine a whole pile of vector data from OSM into one integrated map. While OSM is not particularly high precision (from a bits standpoint) the resulting intersecting points are calculated "perfectly", sometimes with very large mantissas.

I then wrote a piece of code to take a city block from that OSM map and perform some calculations to find the sidewalk calculation. The problem: the four corners of the city block were already very long numbers since they were the result of a CGAL calculation. Thus a long calculation on a long calculation becomes very slow.

The original algorithm took about 36 minutes for a fully optimized build to find all sidewalks in San Diego. That is way too slow, and unusable for our project.

I the put a rounding stage in: fore each corner of the block, I would convert it to a regular 64-bit IEEE float and then back to CGAL, throwing out any "extra" precision that CGAL was saving. Note that the 64-bit float already gives me better than 1 millimeter precision, which is more than overkill for a road. The algorithm run on the "simplified" data ran in 67 seconds.

Now there is one danger: if, due to mismatched road locations in OSM or conflicting edits, some of the "blocks" were really tiny (less than 1 mm) CGAL would have correctly built that block using infinite precision, and my "rounding" would have incorrectly reshaped those blocks, perhaps turning them inside out or in some other way damaging them.

So a necessary step to productizing this 'resolution reduction' is to do a sanity check on each resulting block. Fortunately most of the time if the block contains too-small-to-use data, we don't need the data in the first place.

Constitutional Opposition

One part of this post by Daring Fireball on the iPhone SDK licensing agreement made me chuckle:

If you are constitutionally opposed to developing for a platform where you’re expected to follow the advice of the platform vendor, the iPhone OS is not the platform for you. It never was. It never will be.

It inspired me to come up with a new quotable:

If you are constitutionally opposed to developing for a platform where you’re expected to follow the advice of the platform vendor, you should not be a computer programmer.

See also basically every post by Raymond Chen: "just because, in Win98SE2, you could call SomeRandomWin32API with a combination of NULL, -1, and Bill Gate's IQ and get an undocumented behavior that violates all of Microsoft's guidelines for applications development doesn't mean it will continue to work in Windows 7."

Thank You Jeeves, That Will Be All

The other day I went in to discover why a new piece of scenery code had mysteriously stopped working. Eventually I came to this:

(p,path.size()/2,def,degree,inExtrudeFunc,
inObjectFunc,inChecker,ag_mode_draped_obj);

Ah! Now it all makes sense. The code should have read:

AG_extrude_string(p,path.size()/2,def,degree,inExtrudeFunc,
inObjectFunc,inChecker,ag_mode_draped_obj);

After having done a global search, clearly I had hit the space bar by accident, nuking my function call. The charming thing is that C++ doesn't question why I have a giant list of paranthetical "stuff", it just blissfully compiles it into an expression that does...well, pretty much nothing.

Some of my other favorite C++ isms:

case a: do_it(); break;
b: do_x(); break; // no case, not illegal - now "b'' is a label!
defaultl: do_more(); break; // typo in default? That's a label too!

Of course we are all familiar with the fun that emerges from swapping = and ==. And having a stray semi-colon never hurt anything.

Propsman had an apt characterization: C++ is like an overly polite butler. "A...hamburger on the rocks, Sir? Certainly, Sir, I'll bring you one directly..."