Here are some notes for code that I fear may confuse. (especially me)
The array sp directly implements BNF 2. There are no terminals. Routines must know that sp[0] is E and sp[255] is F. Just E and F are ‘pure’. There are probably no shared nodes except the first 28 which we initialize. Each of these is a node with 8 pure children. Note the recursive data structure; sp[0] and sp[255] each refer to themselves.
The recursive routine ball takes a 3D ball by radius and center coordinates and returns a node of depth d approximating that ball. The resulting octree denotes the intersection of the ball and [−1, 1]3.
The routine tst takes a ball spec, fills in the octree, and probes a full mesh and reports cubie centers that are reported on the wrong side of the boundary.