Quadtrees

[std10093]

So, we use for the same exact reason. I just have a function inside the loop you have.

>If a halfspace intersects the box so that all points in the box are still inside the box, the halfspace is not interesting.
I suspect, you did not want to write the italics and you have made a typo, because I do not get it.

Obviously I thought that - as you said before - A halfspace is only relevant to the polytope inclusion test, if it intersects the current bounding box.

So, I am going to rename the function to ... hmm..I can't think of anything short and illustrative now..

[Nominal Animal]

I suspect, you did not want to write the italics and you have made a typo, because I do not get it.

D'oh. Make the latter "box" "halfspace": "If a halfspace intersects the box so that all points in the box are still inside the halfspace, the halfspace is not interesting."

Obviously I thought that - as you said before - A halfspace is only relevant to the polytope inclusion test, if it intersects the current bounding box.

Yes, but not all halfspaces that intersect the current bounding box are relevant to the polytope inclusion test.

If the halfspace plane is such that it coincides with a corner, edge, or face of the box, but otherwise the box is inside the halfspace, then that halfspace is irrelevant to the current box: There are no points that test "outside that halfspace" within the current box.

This is exactly the case (corners_outside == 0).

(It takes quite a bit of thought to grasp correctly. Like I said, I messed this detail up earlier. It is possible I still haven't thought it out right, so I recommend you don't take my word for it.)

So, I am going to rename the function to ... hmm..I can't think of anything short and illustrative now..

Some variation of "halfspace relevant" to box?

By the way, I personally firmly believe that stuff like that -- that the code uses names and comments that drive the reader and future programmers towards the correct concepts, and perhaps even their intuition towards the correct assumptions -- is an important part of that which makes great code great.

The more I work with code, the more weight I place on the long-term viability and maintainability of it. It is certainly the direction I wish to develop my own skills.

[std10093]

I think I have the code working..
However, I have as test case the unit cube. In your code, does it produce a swallow or normal-depth tree? Because in mine, it goes on and on! And I think that the reason is my logic and not some logical bug.
[EDIT]
Same for the tetrahedron.
[/EDIT]

---

You have an array of the halfspaces and you re-arrange them, thus you only use one array and the cost of re-arrangement.

I have my node as this right now:

struct Node {
    Point3d split;      // Split point
    BoundingBox* box;   // Pointer to bounding box. Free it after construction

    // -1 means that all halfspaces intersect the box (for example, in the root)
    std::vector<int> relevantHalfspaces;    // Indexes of halfspaces that intersect the current box. The indexes are to the halfspaces vector of the
                                                         // polytope
    Node* child[1<<Dimension];              // Pointers to children
};

Where the vector holds the indexes of the array(the same array as yours) which allows for getting away form the re-arrangement cost, but requires some space. And if you think that our tree might get big, then we will have some space multiplied by the number of nodes!
Now, I am thinking of it, I have to destroy the vectors, in all inner leaves.

What do you say about my approach? Are you happy from yours?

[Nominal Animal]

In your code, does it produce a swallow or normal-depth tree?

What are you using as the maximum number of halfspaces per leaf node? I only test with 3 or larger, as smaller limits generate unnaturally deep trees; the intersection of two or three halfspaces are quantized to least-significant-unit sized bounding boxes... I definitely consider 3 the smallest sane limit to test.

For a cube, I get one inner node with eight leaf nodes containing three halfspaces per leaf. (This is with perturbing the vertices very slightly, as otherwise my algorithm detects that three of the sides are coplanar with the box, and keeps only the three halfspaces, all in the root node.)

For the (0,0,0),(1,0,0),(0,1,0),(0,0,1) tetrahedron, I get two inner nodes; one with split point (0.5000000000000002,0.5000000000000002,0.500000000 0000002) with four leaf nodes completely outside the polytope, three leaf nodes with three halfspaces, and one inner node child. The inner node child has split point (0.2500000000000001,0.2500000000000001,0.250000000 0000001), with four leaf nodes with three halfspaces, three leaf nodes with two, and one with one halfspace.

The reason for the split point coordinates not being exactly half and a quarter is the nextafter() logic we discussed above: in my implementation, the upper limit is exclusive, i.e. nextafter(maxCoord, HUGE_VAL); due to rounding, the mean of the bounds is just above half.

struct Node {
    Point3d split;      // Split point
    BoundingBox* box;   // Pointer to bounding box. Free it after construction

Why would you store the bounding box in the node structure? It is only needed during construction! For the bounding box, I actually use six variables in the recursive function call, as it is easier to construct the bounding box using separate variables for the next recursive call.

    // -1 means that all halfspaces intersect the box (for example, in the root)
    std::vector<int> relevantHalfspaces;    // Indexes of halfspaces that intersect the current box. The indexes are to the halfspaces vector of the
                                                         // polytope
    Node* child[1<<Dimension];              // Pointers to children
};

Hmm. You need either the halfspaces, or the child nodes, never both.

To be exact, there are four different types of nodes you need to handle:

• Inner node -- has split point coordinates, and 2^dimensions child nodes
  
• Outside leaf node -- leaf node where all points in the bounding box are outside the polytope (I use NULL pointer for this)
  
• Inside leaf node -- leaf node where all points in the bounding box are inside the polytope (I use an easily-detected normal leaf node for this)
  
• Normal leaf node -- contains or refers to halfspaces which tell which points in the bounding box are inside the polytope, and which are outside

The root node, or the structure containing the root node to the tree, must also have the bounding box of the polytope itself, so that you don't need to test points you know for sure are outside the polytope.

What do you say about my approach?

Aside from the superfluous bounding box in the node structure, I worry that having indexes to the halfspaces (4 or 8 bytes per index) rather than the halfspace itself (16 or 32 bytes per halfspace) is going to be a slowdown.

The reason is cache behaviour. When you test for point inclusion in the polytope, you first need to descend through the tree. This is unavoidable, and a cost we just have to accept.

However, in my case, with the halfspaces contained within the leaf, that data will be cached, and immediately available to the CPU, to determine whether the point is inside the polytope or not.

In your case, the halfspaces are rarely stored near each other, so the CPU has to cache parts of the halfspace array too.

Remember: this implementation makes sense only if lots of halfspaces are needed to define the convex polytope, so we cannot assume all of the halfspace array will be cached from point inclusion test to another.

Obviously, only comparing the real-world results from both approaches would tell for sure.. but this is something I did consider when designing my node structures. So, I'm confident (that the indirect method is measurably slower for real-world use cases), but I don't know it for a fact.

---

As to my own algorithm implementation, I'm quite happy with it. I'm certain there is room for improvement, but it's efficient, clean, and easy-to-read as is. Before posting it, I would do a final pass of checking the comments and testing it with some random polytopes, but that's it.

Improving the split point logic from just halving would be interesting to do research on, but aside from a pure research topic, is not worth the effort.

(I suspect that first finding a case where the tree structure generated by the halving split point logic is clearly poor, would also provide enough information to devise a way to easily avoid these poor cases. I believe they occur only in cases where a small set of vertices defines many faces/halfspaces. For regular polytopes, the halving logic works really well.)

Methods of extending this to nonconvex polytopes are also very interesting, and perhaps more relevant to real-world use cases.

[std10093]

I see your points.

Except this one.

For the (0,0,0),(1,0,0),(0,1,0),(0,0,1) tetrahedron, I get two inner nodes; one with split point (0.5000000000000002,0.5000000000000002,0.500000000 0000002) with four leaf nodes completely outside the polytope, three leaf nodes with three halfspaces, and one inner node child. The inner node child has split point (0.2500000000000001,0.2500000000000001,0.250000000 0000001), with four leaf nodes with three halfspaces, three leaf nodes with two, and one with one halfspace.

For this I am getting this:
V = 4, epsilon = 1e-05, t = 3 <-- The goal is at most 3 halfspaces per leaf
Vertices
0,0,0
1,0,0
0,1,0
0,0,1
Halfspaces
-0 -0 -1 -0
0 -1 0 0
-1 -0 -0 -0
0.577350269189626 0.577350269189626 0.577350269189626 0.577350269189626

and then the octree, creates another level with eight children, as expected and the 6th child (i.e. the pre-last), creates more levels and goes one forever.
The 6th child has this data:
0.250000000050000,0.000000000000000,0.000000000000 000 //split point
Bounding box is: ( 0.000000000000000, 0.000000000000000, 0.000000000000000 ) - ( 0.500000000100000, 1.000000000100000, 1.000000000100000 )
and it intersects all halfspaces.
And then the 6th child of the upper 6th child has this:
0.125000000075000,0.000000000000000,0.000000000000 000
Bounding box is: ( 0.000000000000000, 0.000000000000000, 0.000000000000000 ) - ( 0.250000000150000, 1.000000000200000, 1.000000000200000 )
and it intersects all halfspaces too!

And the icosahedron and the unit cube had the same problem. In the 6th child, they would create children again and again, forever!

[Nominal Animal]

Vertices and halfspaces are as expected. I can even see from the split coordinates that even the initial bounding box should be correct (or at least the same as I have).

The 6th child has this data:
0.250000000050000,0.000000000000000,0.000000000000 000 //split point
Bounding box is: ( 0.000000000000000, 0.000000000000000, 0.000000000000000 ) - ( 0.500000000100000, 1.000000000100000, 1.000000000100000 )

Both that split point, and that bounding box, strike me as very odd.

If your original bounding box is [0,1) for all coordinates, then that bounding box is obviously wrong, right? The first-level children should have either [0,0.5) or [0.5,1) for each coordinate for the bounding box, with split point coordinates being either 0.25 or 0.75, never zero.

If you print out the first level child nodes -- their split point coordinates and bounding box coordinates, what do you get? Draw a 3D image (using e.g. the isometric perspective I mentioned above), it helps me a lot to visualize the situation correctly.

[std10093]

That's what I spotted too and that's why I asked you.

Here is where I am confused: When computing the bounding box of the 7th (for example) child, you have this formula: (split.x, split.y, split.z) - (max.x, max.y, max.z).
I guess, split point is from the 7th child and max from the parent, correct?
On the same child, we have to calculate the split point by this formula: split.x = min.x + 0.5*(max.x - min.x); (same for y and z).
The min and max come from the parent, correct?

I am testing, but something goes wrong..

The split point is calculated by this formula:

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

and the max box by these formulas:

0: (min.x, min.y, min.z) - (split.x, split.y, split.z)
1: (split.x, min.y, min.z) - (max.x, split.y, split.z)
2: (min.x, split.y, min.z) - (split.x, max.y, split.z)
3: (split.x, split.y, min.z) - (max.x, max.y, split.z)
4: (min.x, min.y, split.z) - (split.x, split.y, max.z)
5: (split.x, min.y, split.z) - (max.x, split.y, max.z)
6: (min.x, split.y, split.z) - (split.x, max.y, max.z)
7: (split.x, split.y, split.z) - (max.x, max.y, max.z)

and there seem to be no typos..

[Nominal Animal]

When computing the bounding box of the 7th (for example) child, you have this formula: (split.x, split.y, split.z) - (max.x, max.y, max.z).
I guess, split point is from the 7th child and max from the parent, correct?

Incorrect. Both are from parent.

Think about it for a minute. The parent is a box, with split point somewhere inside it; the common corner for all the child boxes. Where are the edges for the eight child boxes?

[std10093]

Incorrect. Both are from parent.

Think about it for a minute. The parent is a box, with split point somewhere inside it; the common corner for all the child boxes. Where are the edges for the eight child boxes?

Well, I have the code with "both are from parent", because of the thought you just said, but it will just not work. (giving the [0,05) and [0.5,1) you mentioned before.

[Nominal Animal]

Consider the following. I'll use the notation 1+ to mean the smallest number greater than 1.

Let's assume current bounding box is x=[0,1+), y=[0,1+), and z=[0,1+). We discover that there are more relevant halfspaces that can be referred to in a single node, so the current node needs to be split.

We split the current bounding box at the mean of bounding box coordinates (halving the bounding box along each axis), so the current split point will be (0.5+, 0.5+, 0.5+).

So, we have

min = ( 0, 0, 0 )
max = ( 1+, 1+, 1+ )
split = ( 0.5+, 0.5+, 0.5+ )

for the current box, right?

The bounding boxes for each child will be:

0: (  min.x,   min.y,   min.z) - (split.x, split.y, split.z)
1: (split.x,   min.y,   min.z) - (  max.x, split.y, split.z)
2: (  min.x, split.y,   min.z) - (split.x,   max.y, split.z)
3: (split.x, split.y,   min.z) - (  max.x,   max.y, split.z)
4: (  min.x,   min.y, split.z) - (split.x, split.y,   max.z)
5: (split.x,   min.y, split.z) - (  max.x, split.y,   max.z)
6: (  min.x, split.y, split.z) - (split.x,   max.y,   max.z)
7: (split.x, split.y, split.z) - (  max.x,   max.y,   max.z)

just like in your post #67; I only added some spaces to make it easier to read.

I recommend you modify your recursive function to just output the bounding box it receives, and not recurse any further. The above do produce

0: ( 0.0,  0.0,  0.0 ) - ( 0.5+, 0.5+, 0.5+ )
1: ( 0.5+, 0.0,  0.0 ) - ( 1.0+, 0.5+, 0.5+ )
2: ( 0.0,  0.5+, 0.0 ) - ( 0.5+, 1.0+, 0.5+ )
3: ( 0.5+, 0.5+, 0.0 ) - ( 1.0+, 1.0+, 0.5+ )
4: ( 0.0,  0.0,  0.5+) - ( 0.5+, 0.5+, 1.0+ )
5: ( 0.5+, 0.0,  0.5+) - ( 1.0+, 0.5+, 1.0+ )
6: ( 0.0,  0.5+, 0.5+) - ( 0.5+, 1.0+, 1.0+ )
7: ( 0.5+, 0.5+, 0.5+) - ( 1.0+, 1.0+, 1.0+ )

I definitely do not see your (0.0,0.0,0.0)-(0.5+,1.0+,1.0+) in the above list. (Which would be very surprising anyway, as it would fill the upper half of the box.)

The bounding box for the child 6 will be (0.0,0.5+,0.5+)-(0.5+,1.0+,1.0+), and if it is subdivided further (it shouldn't), it's split point will be (0.25+,0.75+,0.75+).

Don't get frustrated. I know this works.

[std10093]

I was modifying things yesterday and now I have the same boxes as you, but on the first child, the tree starts to get really deep. Is that normal? I think not! I have to investigate! I left the tree grow only one level and I saw that 1st, 2nd, 3rd and 4th child intersect all the halfspaces, thus need to be split.
Here is my octree, grown up until one level only.

Octree
0.250002500050000,0.250002500050000,0.250002500050000
Bounding box is: ( 0.000000000000000, 0.000000000000000, 0.000000000000000 ) - ( 0.500005000100000, 0.500005000100000, 0.500005000100000 )
/* That's the root */
0.500005000000000,0.500005000000000,0.500005000000000
Bounding box is: ( 0.000000000000000, 0.000000000000000, 0.000000000000000 ) - ( 1.000010000000000, 1.000010000000000, 1.000010000000000 )
0.750007500050000,0.250002500050000,0.250002500050000
Bounding box is: ( 0.500005000000000, 0.000000000000000, 0.000000000000000 ) - ( 1.000010000100000, 0.500005000100000, 0.500005000100000 )
0.250002500050000,0.750007500050000,0.250002500050000
Bounding box is: ( 0.000000000000000, 0.500005000000000, 0.000000000000000 ) - ( 0.500005000100000, 1.000010000100000, 0.500005000100000 )
0.750007500050000,0.750007500050000,0.250002500050000
Bounding box is: ( 0.500005000000000, 0.500005000000000, 0.000000000000000 ) - ( 1.000010000100000, 1.000010000100000, 0.500005000100000 )
0.250002500050000,0.250002500050000,0.750007500050000
Bounding box is: ( 0.000000000000000, 0.000000000000000, 0.500005000000000 ) - ( 0.500005000100000, 0.500005000100000, 1.000010000100000 )
0.750007500050000,0.250002500050000,0.750007500050000
Bounding box is: ( 0.500005000000000, 0.000000000000000, 0.500005000000000 ) - ( 1.000010000100000, 0.500005000100000, 1.000010000100000 )
0.250002500050000,0.750007500050000,0.750007500050000
Bounding box is: ( 0.000000000000000, 0.500005000000000, 0.500005000000000 ) - ( 0.500005000100000, 1.000010000100000, 1.000010000100000 )
0.750007500050000,0.750007500050000,0.750007500050000
Bounding box is: ( 0.500005000000000, 0.500005000000000, 0.500005000000000 ) - ( 1.000010000100000, 1.000010000100000, 1.000010000100000 )

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

EDIT: Wait, I found a bug.. No, nothing changed, have to look harder!

In the function that tests if a halfspace is relevant to box, if I modify it to this:

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

the tree becomes a tree of one level only, but this is wrong, as you have discussed...

[Nominal Animal]

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₀, x₁, ..., x_n-1 )
y_MIN = MIN( y₀, y₁, ..., y_n-1 )
z_MIN = MIN( z₀, z₁, ..., z_n-1 )
x_MAX = MAX( x₀, x₁, ..., x_n-1 )
y_MAX = MAX( y₀, y₁, ..., y_n-1 )
z_MAX = MAX( z₀, z₁, ..., z_n-1 )

with x_i, y_i, z_i being the exact coordinates for vertex i, then the initial bounding box is

boundingbox.min.x = xMIN;
boundingbox.min.y = yMIN;
boundingbox.min.z = zMIN;

boundingbox.max.x = nextafter(xMAX, HUGE_VAL);
boundingbox.max.y = nextafter(yMAX, HUGE_VAL);
boundingbox.max.z = nextafter(zMAX, 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,

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

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.

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.

if I modify it to this:

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.

---

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.)

[std10093]

>What does your code do when there are no relevant halfspaces? I think you end up recursing forever in that case!
I am afraid that is not the case... I am going to continue searching and if I found nothing I am going to compare it with your implementation and hope that everything will work ok.

[std10093]

The tree is becoming huge and I can not find the bug yet, so I just left the tree grow one level, thus 9 nodes in total, 1 root and 8 children and tested the inclusion test.
I have
Test point = 1.100000000000000,0.500000000000000,0.700000000000 000
and I am getting the answer yes. I descend to the 5th child, were the first three halfspaces are relevant

-0.000000000000000,-0.000000000000000,-1.000000000000000 <- that's the coordinates
 and -0.000000000000000 <- that's d

0.000000000000000,-1.000000000000000,0.000000000000000
 and 0.000000000000000

-1.000000000000000,-0.000000000000000,-0.000000000000000
 and -0.000000000000000

Shouldn't I receive no? Since the x of the inclusion test is out of the unit cube. I have the polytope of the tetrahedron.

[Nominal Animal]

Shouldn't I receive no?

Your test point seems to be outside the polytope bounding box, so it should not even descend into the tree.

Should I post my code yet?

I've been playing with implementing a better vertices-to-halfspaces algorithm than the naive O(n⁴) one currently used; especially one that calculates the convex hull of the given points (and handles degenerate cases well). I can distract myself with that, if you don't want to see my code (yet).