I’m building a CSG library because eventually I’d like to have a native Macintosh editor drive my MakerBot. Building a CSG algorithm that uses BSP trees to handle the 3 dimensional math is pretty trivial; the algorithm I’m using came from the paper “Improved Binary Space Partition Merging” by Lysenko et.al., which itself uses the convex hull algorithm outlined in “Linear Programming and Convex Hulls Made Easy” by Seidel.

The problem is, once you have a BSP representation of your object, you need to turn it back into a collection of polygons. And that means a good 2D polygon manipulation engine. (I’m converting directly from a BSP tree to polygons; the way this works is to define a cube bounding the entire BSP tree (by picking one that is “big enough”), then recursing down the BSP tree structure, splitting that cube as you go. At the bottom, you throw away or return the split cube, and “glue it back together” as you come back up the tree. This implies that on return you have to glue arbitrary polygon shapes on both sides of a split plane together–and on the polygon surface, that means doing a 2 dimensional XOR operation on the polygons there.)

So I’ve been going through the book “Computational Geometry” trying to implement the intersection algorithm outlined in Chapter 2 there.

And I have a few observations.

(1) The exposition of the algorithm is geared more towards a general discussion of the algorithm rather than providing a complete description of the algorithm in a way that is suitable for implementation. That means implementing the algorithm requires a fundamental understanding of what’s going on; you can’t just look at the outline of the algorithm and convert it to code. And you need to discover using various data sets the gaps in your implementation and find where things go south.

(2) The algorithm fundamentally relies on the data structure in section 2.2, which uses half-edges to represent each ring. The half edges must be complete: that is, if you have just a naked polygon sitting in an otherwise empty space, you need both the inner ring representing the inside, and an outer ring of edges which represent the outer boundary. You can’t get away with unpaired half-edges and easily make this thing work. No shortcuts.

(That’s because when two half edges intersect (figure 2.5, section 2.3), you need to tie together the half edges that meet at the point–and sometimes an inner edge pairs with an outer edge and visa-versa. Imagine two polygons intersecting: the outer edges of the two polygons form the inner edges of the penetrating part of the two outer polygons of the intersection.)

You can be excused, by the way, in thinking the half planes could be naked and not have a matching half-edge for an outer defined polygon, given the discussion of monotone polygons and polygon triangulation in Chapter 3 works with naked half-edge polygon rings.

(3) The algorithm itself also relies on a sorted tree structure that you can search by edge or by point. In C++ or Objective C, that means some additional work. Fortunately the number of polygons are small enough that it’s not worth the overhead of using a balanced tree structure.

(4) The algorithm goes south when two edges are coincident. The trick appears to be throwing away identical half edges that are from S₁ and from S₂ in Algorithm MapOverlay in section 2.3, at step 1, creating the new doubly-connected edge list D. (In practice you can do this by dropping duplicate edges in step 2, when initializing the event queue Q: as each edge is added to each event in Q, determine if there is already an identical edge already in that event.)

Also, when computing the intersections, if an intersection test creates a new intersection edge that is identical to an existing edge (because the two intersecting edges were coincident but without the same start or end points), throw away the newly constructed edge intersection if they already exist. This means in practice if, while updating the event in Q in algorithm HandleEventPoint in section 2.1 as a result of finding an intersecting edge, when inserting the found segments, verify they aren’t already in the event in FindNewEvent.

(5) The “final step” referred to in the paragraph preceding MapOverlay, discovering which polygon in one overlay subdivision contains the vertices of the other, can be determined by using a sweep algorithm, sweeping on the vertices in both subdivisions S₁ and S2, storing the edges of each in their own tree list. At each point, you can then determine which edge of one subdivision is to the left of the point in the other; that leftmost edge determines which polygon contains our vertex.

I’m sure you could probably also do this during the sweep in step 2 of MapOverlay; in fact, we do capture the half edge to the left of the event point, though that half edge is of D, which is the combination of S₁ and S₂. So it may be a matter of continuing the leftward walk to find representatives of both sets, or simply maintaining two other tree structures at the same time so you can get the leftward edges for D, S₁ and S₂ at the same time.

At some point, after I successfully get an implementation of my XOR algorithm going correctly, I’m going to write an exposition of this algorithm that isn’t so damned obtuse and which can be more easily implemented as code. Because this article (along with the article describing triangulation of a polygon) feels like it was written by a professor who intuitively knew the algorithms, but was not working from a functional example of the algorithms he was describing when writing the paper.

That’s because the gaps and holes in the exposition are big enough to drive a truck through. (For example, in the description FindNewEvent in section 2.1, we have the values L(p), U(p) and C(p) defined; however, from the description of FindNewEvent it’s unclear what we do at step 2: “Insert the event point as an event into Q.” Given an event contains L(p), U(p) and C(p), what is the event we insert when we insert the event point? And let’s skip by the “However, we won’t explain this final step here. It’s a good exercise to test your understanding of the plane sweep approach to design the algorithm yourself” comment later in the chapter.)

So if you were trying to implement a polygon mesh intersection from this discussion and decided you were just dumb, well, no you’re not: I’ve wasted a couple of months of my spare time trying to sort all this out, and I can tell you: you’re not dumb. The description is incomplete.