If you mean that the root of your octree has bounding box

( 0.0000000000, 0.0000000000, 0.0000000000 ) - ( 1.0000100001, 1.0000100001, 1.0000100001 )

and the bounding boxes for its childs are

( 0.0000000000, 0.0000000000, 0.0000000000 ) - ( 0.5000050001, 0.5000050001, 0.5000050001 )

( 0.5000050000, 0.0000000000, 0.0000000000 ) - ( 1.0000100001, 0.5000050001, 0.5000050001 )

( 0.0000000000, 0.5000050000, 0.0000000000 ) - ( 0.5000050001, 1.0000100001, 0.5000050001 )

( 0.5000050000, 0.5000050000, 0.0000000000 ) - ( 1.0000100001, 1.0000100001, 0.5000050001 )

( 0.0000000000, 0.0000000000, 0.5000050000 ) - ( 0.5000050001, 0.5000050001, 1.0000100001 )

( 0.5000050000, 0.0000000000, 0.5000050000 ) - ( 1.0000100001, 0.5000050001, 1.0000100001 )

( 0.0000000000, 0.5000050000, 0.5000050000 ) - ( 0.5000050001, 1.0000100001, 1.0000100001 )

( 0.5000050000, 0.5000050000, 0.5000050000 ) - ( 1.0000100001, 1.0000100001, 1.0000100001 )

with the split point in the center of each bounding box, it looks right to me.

Hmm.. You do not enlarge the bounding box in the recursive call, do you? My argument about the necessity of enlarging the bounding box as the upper limit is exclusive and not inclusive applies only to the initial bounding box. In other words, if

x_{MIN} = *MIN*( x_{0}, x_{1}, ..., x_{n-1} )

y_{MIN} = *MIN*( y_{0}, y_{1}, ..., y_{n-1} )

z_{MIN} = *MIN*( z_{0}, z_{1}, ..., z_{n-1} )

x_{MAX} = *MAX*( x_{0}, x_{1}, ..., x_{n-1} )

y_{MAX} = *MAX*( y_{0}, y_{1}, ..., y_{n-1} )

z_{MAX} = *MAX*( z_{0}, z_{1}, ..., z_{n-1} )

with *x*_{i}, *y*_{i}, *z*_{i} being the exact coordinates for vertex *i*, then the *initial* bounding box is

Code:

boundingbox.min.x = x_{MIN};
boundingbox.min.y = y_{MIN};
boundingbox.min.z = z_{MIN};
boundingbox.max.x = nextafter(x_{MAX}, HUGE_VAL);
boundingbox.max.y = nextafter(y_{MAX}, HUGE_VAL);
boundingbox.max.z = nextafter(z_{MAX}, HUGE_VAL);

In the recursive calls, the bounding box coordinates are treated as exact, no adjustment with nextafter(). The split point is calculated as the center (mean) of the bounding box,

Code:

split.x = boundingbox.min.x + 0.5 * (boundingbox.max.x - boundingbox.min.x);
split.y = boundingbox.min.y + 0.5 * (boundingbox.max.y - boundingbox.min.y);
split.z = boundingbox.min.z + 0.5 * (boundingbox.max.z - boundingbox.min.z);

which is numerically more stable than just halving their sum.

Because our lower limit of bounding box is inclusive, but the upper limit exclusive, the coordinate range *min*,*split* does not include the split point, but *split*,*max* does.

That also means that when descending the tree (in the polytope point inclusion function), we use

Code:

child = 1 * (testpoint.x >= split.x)
+ 2 * (testpoint.y >= split.y)
+ 4 * (testpoint.z >= split.z);

to determine which child (0..7) to descend into.

Every unique point describable using floating-point coordinates is inside exactly one of the child nodes. The split point of the parent, will be inside the last child.

If we did not do the initial bounding box upper limit adjustment, we would exclude those points of the polytope, that had the exact maximum coordinates (maximum coordinates from any vertex) -- again, because the upper limit is *exclusive*.

While that slice is infinitesimally thin to mathematicians (0.999... being equal to 1), it is measurable for us programmers; and that is exactly what *nextafter()* was designed to help with. I've seen how ignoring this kind of seemingly irrelevant details cause nastily visible seams and discontinuities in raytracing for example. Me being me, I find it important to do these details *right*, if you do it at all.

You make your own choices. As I've demonstrated, I am often wrong anyway.

Originally Posted by

**std10093**
If I have understood well, the* split point is being computed from the current box* and the current box from the parent's data.

Exactly.

Originally Posted by

**std10093**
if I modify it to this:

Code:

if(cornersInside == 0 || cornersOutside == 8 || cornersOutside == 0)
then halfspace is relevant

Ah-ha!

What does your code do when there are no relevant halfspaces? I think you end up recursing forever in that case!

For the above tests, I have either expanded the polytope bounding box slightly, or perturbed the vertex coordinates slightly, so that my code does not find that three of the halfspaces are coplanar with the polytope bounding box -- for *both* the (0,0,0)-(1,0,0) cube, and the (0,0,0),(1,0,0),(0,1,0),(0,0,1) tetrahedron.

You see, in both cases the X=0, Y=0, and Z=0 faces are coplanar to the halfspace planes, and can therefore be ignored.

I am assuming that your code dismisses those halfspaces too, at least in the recursive case. You end up calling your recursive function for at least one of the sub-boxes, with no relevant halfspaces (because the logic has culled the halfspaces coplanar to the bounding box out).

If you change the test logic, you retain those otherwise culled halfspaces, and therefore avoid the bug.

Originally Posted by

**std10093**
The tree becomes a tree of one level only

That is absolutely correct, if the cube or tetrahedron coordinates are exact integers.

The algorithm notes that three of the cube or tetrahedron faces/halfspace planes are coplanar with the bounding box (lower limits), and therefore irrelevant. This yields only three (cube) or one (tetrahedron) relevant halfspaces *for the initial bounding box* (entire polytope), which is obviously storeable in one node; thus yielding a "tree" with exactly one node.

This is obviously desirable for real life code, as it speeds polytope inclusion tests up, and produces shallower/faster trees.

To avoid this effect, I either enlarge the bounding box, in which case all halfspaces are inside the initial bounding box, not coplanar with any part of it, or perturb the vertex coordinates by a fraction (typically by a few hundredths).

When I test my algorithm, I test all three cases, because all three are possible in real-life data. (Well, except enlarging the initial bounding box, as I calculate that directly from the vertex coordinates.)