Quadtrees

[Nominal Animal]

Good catch, std10093. I'll investigate.

Edit 1: There was a typo when defining v3. The last component should be p3.z-p1.z, not p3.z-p3.z.

Edit 2: For the purposes of implementing this algorithm, we need to distinguish not only between "outside" and "inside" a half-space, but also "on the dividing plane". My code in #12 requires origin (0,0,0) to be inside the polyhedron, but you put one of the vertices there. That is not "inside", that's "on the dividing plane". That's why it fails.

In the mean time, you can use a much more robust half-space logic instead. This only requires the vertices to form a convex polyhedron, and that's it. (The origin does not need to lie within the polyhedron).

Loop over every triplet of vertex points p1, p2, p3:

    Calculate n' = (p2 - p1) × (p3 - p2)
    Calculate nsqr = n' · n'
    If nsqr < eps*eps, skip to next triplet (as p1, p2, p3 are [almost] collinear)

    Calculate n = n' / sqrt(nsqr)  ( == n' / || n' || )
    Calculate d = n · p1 ( == n · p2 == n · p3 )

    If half-space (n, d) is already known, skip to next triplet.

    Calculate the number of points inside, at, and outside the half-space.
    For each vertex p, p != p1, p != p2, p != p3:
        If (p · n - d > eps):
            p is outside the half-space
        Else, if (p · n - d < -eps):
            p is inside the half-space
        Else:
            p is on the space-dividing plane.

    If there was:
        at least one point inside the half-space, and
        no points outside the half-space, and
        at least three points on the space-dividing plane,
    then this half-space is one that defines the convex polyhedron;
    add it to the list of known half-spaces.

End loop

The corresponding halfspaces() function is

/* Given a list of 'points' vertices in 'point',
 * defining a convex polyhedron,
 * such that origin (0,0,0) is inside the polyhedron,
 * return the number of half-spaces needed to
 * define the polyhedron.
 *
 * Return 0 with errno set if an error occurs.
*/
size_t halfspaces(halfspace_t  **const listptr,
                  size_t        *const sizeptr,
                  const point_t *const point,
                  const size_t         points,
                  const double         epsilon)
{
    const double  esqr = epsilon * epsilon;
    halfspace_t  *list, curr;
    size_t        size, used;
    size_t        i, i1, i2, i3;

    /* Sanity check. */
    if (!listptr || !sizeptr || !point || points < 3) {
        errno = EINVAL;
        return (size_t)0;
    }

    /* If listptr points to NULL, then size is 0. */
    list = *listptr;
    if (list)
        size = *sizeptr;
    else
        size = 0;

    /* No halfspaces defined yet. */
    used = 0;

    for (i1 = 0; i1 < points; i1++) {
        const point_t  p1 = point[i1];

        for (i2 = 0; i2 < points; i2++) {
            if (i2 != i1) {
                const point_t  p2 = point[i2];

                for (i3 = 0; i3 < points; i3++) {
                    if (i3 != i1 && i3 != i2) {
                        const point_t  p3 = point[i3];

                        const point_t  v2 = { p2.x - p1.x, p2.y - p1.y, p2.z - p1.z };
                        const point_t  v3 = { p3.x - p1.x, p3.y - p1.y, p3.z - p1.z };

                        const point_t  n = { v2.y * v3.z - v2.z * v3.y,
                                             v2.z * v3.x - v2.x * v3.z,
                                             v2.x * v3.y - v2.y * v3.x };

                        const double   nn = n.x * n.x + n.y * n.y + n.z * n.z;

                        if (nn > esqr) {
                            const double nlen = sqrt(nn);
                            size_t       n_in = 0;  /* Vertices inside the half-space */
                            size_t       n_at = 0;  /* Vertices on the half-space plane */
                            size_t       n_out = 0; /* Vertices outside the half-space */

                            curr.x = n.x / nlen;
                            curr.y = n.y / nlen;
                            curr.z = n.z / nlen;
                            curr.d = p1.x * curr.x + p1.y * curr.y + p1.z * curr.z;

                            /* Already included in the list? */
                            for (i = 0; i < used; i++)
                                if (fabs(curr.x - list[i].x) <= epsilon &&
                                    fabs(curr.y - list[i].y) <= epsilon &&
                                    fabs(curr.z - list[i].z) <= epsilon &&
                                    fabs(curr.d - list[i].d) <= epsilon)
                                    break;
                            if (i < used)
                                continue;

                            /* Count points inside and outside the half-space. */
                            for (i = 0; i < points; i++) {
                                const double d = curr.x * point[i].x
                                               + curr.y * point[i].y
                                               + curr.z * point[i].z
                                               - curr.d;
                                if (d > epsilon)
                                    n_out++;
                                else
                                if (d < -epsilon)
                                    n_in++;
                                else
                                    n_at++;
                            }

                            /* There must be:
                             *  - at least one point inside,
                             *  - at least three points on the plane
                             *  - no points outside
                             * for this to be a valid half-space.
                            */
                            if (n_in < 1 || n_at < 3 || n_out > 0)
                                continue;

                            /* Make room for a new halfspace. */
                            if (used >= size) {
                                size = (used | 15) + 17;
                                list = realloc(list, size * sizeof *list);
                                if (!list) {
                                    errno = ENOMEM;
                                    return (size_t)0;
                                }
                                *listptr = list;
                                *sizeptr = size;
                            }

                            list[used++] = curr;
                        }
                    }
                }
            }
        }
    }

    errno = 0;
    return (size_t)used;
}

Edit 3: Here is the final version of halfspaces() function, using the same logic as outlined earlier in this message, but only using N³/6-N²/2+N/3 iterations (instead of N³).

/* Given a list of 'points' vertices in 'point',
 * defining a convex polyhedron,
 * such that origin (0,0,0) is inside the polyhedron,
 * return the number of half-spaces needed to
 * define the polyhedron.
 *
 * Return 0 with errno set if an error occurs.
*/
size_t halfspaces(halfspace_t  **const listptr,
                  size_t        *const sizeptr,
                  const point_t *const point,
                  const size_t         points,
                  const double         epsilon)
{
    const double  esqr = epsilon * epsilon;
    halfspace_t  *list, curr;
    size_t        size, used;
    size_t        i, i1, i2, i3;

    /* Sanity check. */
    if (!listptr || !sizeptr || !point || points < 3) {
        errno = EINVAL;
        return (size_t)0;
    }

    /* If listptr points to NULL, then size is 0. */
    list = *listptr;
    if (list)
        size = *sizeptr;
    else
        size = 0;

    /* No halfspaces defined yet. */
    used = 0;

    for (i1 = 0; i1 < points - 2; i1++) {
        const point_t  p1 = point[i1];

        for (i2 = i1 + 1; i2 < points - 1; i2++) {
            const point_t  p2 = point[i2];

            for (i3 = i2 + 1; i3 < points; i3++) {
                const point_t  p3 = point[i3];

                const point_t  v2 = { p2.x - p1.x, p2.y - p1.y, p2.z - p1.z };
                const point_t  v3 = { p3.x - p1.x, p3.y - p1.y, p3.z - p1.z };

                const point_t  n = { v2.y * v3.z - v2.z * v3.y,
                                     v2.z * v3.x - v2.x * v3.z,
                                     v2.x * v3.y - v2.y * v3.x };

                const double   nn = n.x * n.x + n.y * n.y + n.z * n.z;

                if (nn > esqr) {
                    const double nlen = sqrt(nn);
                    size_t       n_in = 0;  /* Vertices inside the half-space */
                    size_t       n_at = 0;  /* Vertices on the dividing plane */
                    size_t       n_out = 0; /* Vertices outside the half-space */

                    curr.x = n.x / nlen;
                    curr.y = n.y / nlen;
                    curr.z = n.z / nlen;
                    curr.d = p1.x * curr.x + p1.y * curr.y + p1.z * curr.z;

                    /* Already included in the list? Inverse included? */
                    for (i = 0; i < used; i++)
                        if (fabs(curr.x - list[i].x) <= epsilon &&
                            fabs(curr.y - list[i].y) <= epsilon &&
                            fabs(curr.z - list[i].z) <= epsilon &&
                            fabs(curr.d - list[i].d) <= epsilon)
                            break;
                        else
                        if (fabs(curr.x + list[i].x) <= epsilon &&
                            fabs(curr.y + list[i].y) <= epsilon &&
                            fabs(curr.z + list[i].z) <= epsilon &&
                            fabs(curr.d + list[i].d) <= epsilon)
                            break;
                    if (i < used)
                        continue;

                    /* Count points inside and outside the half-space. */
                    for (i = 0; i < points; i++) {
                        const double d = curr.x * point[i].x
                                       + curr.y * point[i].y
                                       + curr.z * point[i].z
                                       - curr.d;
                        if (d > epsilon)
                            n_out++;
                        else
                        if (d < -epsilon)
                            n_in++;
                        else
                            n_at++;
                    }

                    /* We need at least three points on the half-space plane,
                     * so that they form a true face. */
                    if (n_at < 3)
                        continue;

                    /* If there are points both inside and outside
                     * the half-space, it does not define the polyhedron. */
                    if (n_in > 0 && n_out > 0)
                        continue;

                    /* If we had points only outside the half-space,
                     * we need to reverse the half-space. Then,
                     * all those points were actually inside. */
                    if (n_out > 0 && n_in == 0) {
                        curr.x = -curr.x;
                        curr.y = -curr.y;
                        curr.z = -curr.z;
                        curr.d = -curr.d;
                        n_in = n_out;
                        n_out = 0;
                    }

                    /* If there are points outside the half-space,
                     * or no points inside the half-space,
                     * then this cannot be one that defines the polyhedron. */
                    if (n_in < 1 || n_out > 0)
                        continue;

                    /* Make room for a new halfspace. */
                    if (used >= size) {
                        size = (used | 15) + 17;
                        list = realloc(list, size * sizeof *list);
                        if (!list) {
                            errno = ENOMEM;
                            return (size_t)0;
                        }
                        *listptr = list;
                        *sizeptr = size;
                    }

                    list[used++] = curr;
                }
            }
        }
    }

    errno = 0;
    return (size_t)used;
}

The difference between this and the code in post #12 in this thread, is that this code does not rely on any extra rules for origin.

If you fix the line defining v3 in the code in post #12, that code will work; just remember that that code requires you to specify vertices that form a polyhedron with origin inside (and not on the face, edge, or vertex) of the polyhedron.

[std10093]

I agree with all you said. I am playing around with the code and check this

-1 1 1
1 1 -1
-1 -1 -1
1 -1 1

4 vertices: x y z
           -1.000000     1.000000     1.000000
            1.000000     1.000000    -1.000000
           -1.000000    -1.000000    -1.000000
            1.000000    -1.000000     1.000000
4 half-spaces: x y z d
        -0.577350  0.577350 -0.577350     0.577350
         0.577350  0.577350  0.577350     0.577350
        -0.577350 -0.577350  0.577350     0.577350
         0.577350 -0.577350 -0.577350     0.577350

RUN SUCCESSFUL (total time: 29s)

I want to feel the input and output of this (or any example)! I also now see, that I do not understand how we represent the halfspace..

[Nominal Animal]

I do not understand how we represent the halfspace.

Consider a cube (-1,-1,-1)-(+1,+1,+1):

-1 -1 -1
-1 -1 +1
-1 +1 -1
-1 +1 +1
+1 -1 -1
+1 -1 +1
+1 +1 -1
+1 +1 +1

The six half-spaces are

x = +1   y =  0   z =  0   d = 1
x =  0   y = +1   z =  0   d = 1
x =  0   y =  0   z = +1   d = 1
x =  0   y =  0   z = -1   d = 1
x =  0   y = -1   z =  0   d = 1
x = -1   y =  0   z =  0   d = 1

(x,y,z) of the half-space is the unit normal vector of the face.
d is the closest distance to origin, if the face plane is extended to infinity.
For the above cube, the normals are axis-aligned, and of course one unit away from origin.

(Technically, d is signed distance along the unit normal vector. That means that if the unit normal points away from the direction you're measuring, you measure "negative" distances.)

Now, consider your unit tetrahedron,

0 0 0      (vertex O, at origin)
1 0 0      (vertex A, on the X axis)
0 1 0      (vertex B, on the Y axis)
0 0 1      (vertex C, on the Z axis)

The half-spaces are now

x =  0   y =  0   z = -1   d = 0      (half-space defined by the OAB plane)
x =  0   y = -1   z =  0   d = 0      (half-space defined by the OAC plane)
x = -1   y =  0   z =  0   d = 0      (half-space defined by the OBC plane)
x = √⅓   y = √⅓   z = √⅓   d = √⅓     (half-space defined by the ABC plane)

The ABC plane is the interesting one. It's the one with the largest area (√¾≃0.866, the other three are ½=0.500). The plane is closest to origin at (⅓,⅓,⅓) where it is distance √[(⅓)²+(⅓)²+(⅓)²] = √⅓ from origin.

In the simplest terms, each half-space is just the face plane, but with the plane unit normal vector defined so that the normal points outwards.

If you start with a large enough ball of clay, you can cut it into any convex polyhedron. You need exactly one cut per face. The plane of that cut is obviously the face plane, but it also is exactly the half-space: outside is the clay you need to discard.

If you look at the Wikipedia article on (geometric) planes, it mentions that one way to define a plane is by knowing its normal vector n (i.e. a vector perpendicular to the plane), and knowing one point on the plane p₀. Then, a point p is on the plane if and only if

n·(p - p₀) = 0
⇔ n·p = n·p₀

The right side is exactly equivalent to our use of d, i.e. we use

d = n·p₀
⇔ n·p - d = 0

When the right side is nonzero, it can be either positive or negative. We use that to define our half-spaces. Using ε≥0 as a precision limit, I call

n·p - d > ε: Outside
n·p - d < ε: Inside
Otherwise ε ≥ n·p - d ≥ -ε: On the plane

Also note that it is not absolutely necessary to normalize the n vectors. I often need the normal vectors n for other purposes (visualization, describing distances to the user), so I "automatically" normalize them -- but if you only use them for the half-space inclusion tests, you can omit the square root and division (when constructing the half-spaces).

Does this help at all? If not, I recommend you brush up on your vector calculus and computational geometry. The intersection tests (line-plane, line-sphere, line-cylinder) and point distance tests (point-plane) are visual and not too dull to use as example cases.

[std10093]

I feel wonderful now with your explanation. However, I am going to go for a tour around the links you provided, but now I have to go, since at 6:30 I have a date with my girlfriend and now it is 6:33 and tomorrow she is leaving for vacation in a Greek - of course :P - island.

Thanks

[std10093]

Faces vs facets in a polytope. (polytope vs polyhedron is also a question, since some use it as the same thing. Is it?)
So, the second step of the algorithm says:
If the number of facets |P| ≤ t, then Q stores the
hyperplanes bounding P.

The problem is that I can not feel what the facets are (since I might get confused with the faces and the many links found in google. By the way, the most interesting was faces vs facets and the facets. )

I also looked for an algorithm to compute the number of facets -though that this would help me understand- but nothing. I am actually looking for an example.

[Nominal Animal]

Faces vs facets in a polytope. (polytope vs polyhedron is also a question, since some use it as the same thing.

In this thread, I have used:

• polyhedron == three-dimensional polytope == 3-polytope (as per Wikipedia polytope article)
• face == side of a polytope (as per Wikipedia articles on geometric edges and polytopes

In particular, "side" would be much better than "face" or "facet", since there is no ambiguity, and it also fits in well with the half-space terminology (and the analogy of cutting the polytope out of an D-dimensional ball of clay using flat cutting planes).

Note that one way to define a convex polytope is to define it as the intersection of half-spaces. Using this definition, face, facet, side, plane, and hyperplane at least all refer to the same thing, the (vector n and scalar d) that define the half-space.

The article relies heavily on that definition, and therefore I think it is safe to assume the authors meant roughly face == facet == side == the plane defined by the half-space (i.e. n·p - d = 0), plus or minus some details that do not affect the implementation of the algorithm.

It is important to remember that the plane equation, n·p - d = sum(n_i p_i, i=1..D) - d = 0, holds for any positive integer number of dimensions D. The data we use to define the half-space, also defines exactly one plane (or hyperplane or face or facet); no matter what the term, it refers to the same thing (to us programmers).

So, the second step of the algorithm says:
If the number of facets |P| ≤ t, then Q stores the
hyperplanes bounding P.

It means that when the convex polytope is defined as half-spaces, and the number of half-space planes intersecting the current box is at most t, the current box is a leaf, and should contain the definitions of those half-spaces. (If the current box does not intersect any of the planes defined by the half-spaces, then the entire box is either outside [if any point in the box is outside any of the half-spaces] or inside the convex polytope. This is a rare case, though.)

Since the half-space test is just D-1 multiplications and D+1 additions and two comparisons (assuming non-zero epsilon), in practical applications you can use a pretty large t. To me, something like 8 to 24 (for 3D) sounds sensible; I'd make t a compile-time constant, and have a couple of test cases that can be used to find out a good value for t on each new architecture.

The article is a bit funny how it switches between two dimensions and more than three dimensions. Don't let the terms fool or scare you.

I don't have the brainpower to correctly visualize more than three spatial dimensions, so what I do is I pick the three-dimensional case, and continuously compare it to the two-dimensional case to see how it would generalize to more than three dimensions. First, I try to get a general overview, without getting snagged on the specific meaning on each term. (That was hard to learn; I only learned to do that when I realized how "relaxed" or willy-nilly the use of these terms really is. Even if the terms have specific, detailed meanings, most people just don't use them that way. Often you just have to rely on the context instead, to determine what the author actually meant.)

After I get a general overview, I re-read the article again, to see if there are any "tricks" or important details hidden in the terms. Sometimes they may look like innocent details at first sight, but turn out to be very limiting for a practical implementation. Having a good nights sleep in between seems to help me a lot.

I didn't see any snags or tricks in this article, really. I only felt that the terminology jumped between 2D and >3D in a way unfamiliar to me. It almost felt like .. perhaps one of the reviewers had supplied conflicting comments, not realizing how the article applies to any positive integer number of dimensions, not just 2D, and some of the terms were changed to "make" that reviewer understand?

[std10093]

I want to be sure I understand the terms, that's why I asked.

The halfspace test you mention is referenced in post #4.

I see what you mean by the 2D and other dimensions, but since that does not affect the algorithm, we should not be concerned that much.

[Nominal Animal]

I see what you mean by the 2D and other dimensions, but since that does not affect the algorithm, we should not be concerned that much.

Exactly.

---

There is one implementation detail left we haven't discussed yet: how to construct the tree efficiently, when the defining half-spaces are known.

The initial bounding volume is known from the vertex coordinates; it's bounded by the minimum and maximum vertex coordinates. (In other words, the minimum and maximum coordinates of the bounding box are the minimum and maximum coordinates from the vertices.) All the half-space planes intersect this box.

If you use a recursive function, and start with one list of all half-spaces, then you only need to remove the half-spaces where the plane does not intersect the currect box when recursing into the sub-boxes. In other words, you do not need to scan the entire original list of half-spaces at every box, only those that did intersect the parent box. This very quickly narrows down the set of half-spaces, making the implementation faster.

It is also safe to assume it's best to start with a single-threaded algorithm for this. (You can always try to parallelize it later, or like I do, search for parallelizable alternate algorithms; having a known working one not only makes the development easier, but helps you catch bugs by comparing the results between the different implementations.)

Let's assume the parent box intersected with N_PARENT half-space planes listed in array H. We could always create a duplicate of H we pass to each child we recurse into, but I'm sure you can see how that would balloon the memory use. So, let's try to avoid duplicating the H array. Remember, we may have to recurse up to 2^D children in D dimensions.

Assume we trim H by moving the planes outside the current box to the end of the list. Then, the N half-space planes intersect the current box are listed first, followed by the other N_PARENT - N half-spaces.

In that case, after the recursive call to a sub-box, H still has the exact same N half-spaces in the array, only their order may have changed.

That turns out to be quite acceptable. If we return to our parent caller, the parent observes the same situation. If we do a recursive call to another sub-box, that call will first do the same "trimming" -- rearranging, really -- of the first N half-spaces in H, but not delete or remove any of them. In other words, having the order change during a recursive call is okay.

The recursive function therefore can use a single list of half-spaces; no termporary arrays or such are needed.

(You can make the above thread-safe, by giving each thread their own copy of the half-space list H. I suspect parallelization is not useful or practical for this -- because the data structure construction is not a bottleneck in any way --, but it definitely is possible. The constructed data structure itself can be accessed from multiple threads simultaneously, of course, since the data structure is not modified afterwards.)

At this point, I do have everything I need to write a full implementation of the algorithm, but as I don't need it myself for anything, I'd much rather you (either std10093, or another reader who wishes to implement this) wrote your own. If you run into troubles, I'm quite willing to write and post my version (after you, please ), so we can compare results, and/or see how to fix those troubles.

[std10093]

At this point, I do have everything I need to write a full implementation of the algorithm, but as I don't need it myself for anything, I'd much rather you (either std10093, or another reader who wishes to implement this) wrote your own. If you run into troubles, I'm quite willing to write and post my version (after you, please ), so we can compare results, and/or see how to fix those troubles.

You know that's my spirit too

I have a question, in the time that we are to compute the halfspaces that define the polytope. Won't they be, in number, equal with the number of faces of the polytope?
For example, the polytope defined by 4 vertices:
0 0 0
1 0 0
0 1 0
0 0 1
gives 4 halfspaces, thus we could use an array for the halfspaces or my idea is not good and we should use a list?

[Nominal Animal]

to compute the halfspaces that define the polytope. Won't they be, in number, equal with the number of faces of the polytope?

Yes.

One side of the convex polytope always generates one halfspace.

The plane of the side is the plane of the halfspace.

In fact, the side of the convex polytope ≡ the halfspace.

The only thing that the polytope side does not determine, is which half of the halfspace is inside and which one outside. (If you negate all components of the halfspace, it still specifies the exact same plane, only the "inside" and "outside" are swapped.)

Because the algorithm here needs halfspaces with normal vectors that point outside the polytope, the halfspace generator I've shown earlier in this thread checks the vertices. By definition, the vertices are either "inside" the halfspace, or "on" the plane of the halfspace (in which case that vertex is one of the vertices that defines the side of the convex polytope matching the halfspace).

we could use an array for the halfspaces

Yes, exactly.

The majority of my post #24 explains that you can use that single array when constructing the tree using a recursive function.

I'd use a recursive function, something like

struct node *construct_node(const point_t min, /* Bounding box minimum */
                            const point_t max, /* Bounding box maximum */
                            halfspace_t *const halfspace, /* array */
                            const size_t halfspaces /* size of array */);

where the function first reorders the halfspaces elements in halfspace so that the ones whose planes intersect the min-max bounding box are listed first. The number of those planes is supplied as halfspaces in the recursive calls, if this box is subdivided further.

Actually, when recursing into the sub-boxes, I'd use the vertex list as a basis on where to split the current box: if there is a vertex not at a corner in the current box, I'd split the box there. Thus, I'd need the vertex list, too. The same logic applies to the vertex list as applies to the half-space list -- I only need the vertices inside the current box (vertices at box corners are not needed) --, so really, my function declaration would most likely be

struct node *construct_node(const point_t min, /* Bounding box minimum */
                            const point_t max, /* Bounding box maximum */
                            halfspace_t *const halfspace, /* array of halfspaces */
                            const size_t halfspaces, /* number of halfspaces */
                            point_t *const point, /* array of vertices */
                            const size_t points /* number of vertices */);

You obviously do not need to do it this way. You can construct each node in the tree by checking the full array containing all halfspaces. But, because this one is significantly faster, and really no more difficult to implement, and I'm convinced it works right, I wanted to describe it.

(Darn it, I don't understand why my posts get this long.)

[std10093]

Since what you say is good, the longer your posts are the amount of good things for us to read get bigger. Maybe on the C and C++ thread, you need to be more laconic, but for sure, in these types of thread, where beginners do not hang out, it is great

You said
> it's bounded by the minimum and maximum vertex coordinates. (In other words, the minimum and maximum coordinates of the bounding box are the minimum and maximum coordinates from the vertices.) All the half-space planes intersect this box.
So, for these vertices, the minimum is the 0 0 0 and the maximum? All the other vertices are ok to be meant as maximum ones?
0 0 0
1 0 0
0 1 0
0 0 1

//I am now in the part that I am going to compute the halfspaces, that's why I do not focus a lot in the other steps

[Nominal Animal]

So, for these vertices,
0 0 0
1 0 0
0 1 0
0 0 1

the bounding box is (0,0,0) - (1,1,1).

The minimum x of the bounding box is the smallest x in any of the vertices,
the minimum y of the bounding box is the smallest y in any of the vertices,
and so on for all coordinate axes.

The maximum x of the bounding box is the largest x in any of the vertices,
the maximum y of the bounding box is the largest y in any of the vertices,
and so on for all the coordinate axes.

That's why it's called a bounding box: it describes the bounds (minimum and maximum) for each coordinate.

I am now in the part that I am going to compute the halfspaces, that's why I do not focus a lot in the other steps

Right; sorry for speeding ahead.

Note, you don't need the bounding box if you are generating the halfspace array based on the vertices.

[std10093]

I understood the bounding box concept, thank you (again).

I generate the halfspaces from the vertices, so as you said, I do not need these bounding boxes.

There is no problem in speeding, since everything you wrote will be read/understand by me, but just not now, but when I reach this level

I was thinking of parallel methods too :O // Good thing I had the same thought as you

Ok, now I think I am a bit confused.
The second step is: If the number of facets |P| ≤ t, then Q stores the hyperplanes bounding P.

I am constructing my code, based on this example (in order to test step by step)
0 0 0
1 0 0
0 1 0
0 0 1

which yields 4 halfspaces.
So, I generate the halfspaces, thus the root of the octree contains these four halfspaces. But you suggested t to be 8 or 24. Let's take 8.
Then, 4 < 8, thus in the 2nd step of the algorithm we stop already? Or the boxes do not represent what I think they are? Because of your previous posts, I think they don't.

I have the illusion (?), that the root holds the hyperplanes of the whole polytope and then we divide the polytope to 8 boxes, with every one box containing a part (the 1/8) of the polytope and so one. I think I am on the wrong track.

Moreover, a technical question:
How am I suppose to insert new data in an octree?
For example, in the binary tree, you simply do something like this

Treeptr addtree(Treeptr p, char *w)   /* Insert word w into tree p */
{ int cond;
  if (p == NULL) {                             /* If tree is empty */
    p = malloc(sizeof(Treenode));   /* Allocate space for new node */
                                    /* Allocate space to copy word */
    p->word = malloc((strlen(w)+1) * sizeof(char));
    strcpy(p->word, w);                /* Copy word w to tree node */
    p->left = NULL;           /* Left subtree of new node is empty */
    p->right = NULL;         /* Right subtree of new node is empty */
  }
  else if ((cond = strcmp(w, p->word)) < 0)
                      /* Does word w precede word of current node? */
    p->left = addtree(p->left, w);
                            /* If yes, insert it into left subtree */
  else if (cond > 0)       /* Does it follow word of current node? */
    p->right = addtree(p->right, w);
                           /* If yes, insert it into right subtree */
  /* If it is the same with word of current node, do not insert it */
  return p;                                         /* Return tree */
}

, thus you check if it has to go left or right, but in the octree?

[Nominal Animal]

which yields 4 halfspaces.
So, I generate the halfspaces, thus the root of the octree contains these four halfspaces.

Yup, and you end up with a tree with one node only.

You don't want to drop the limit t too low, because that would mean the boxes get very small near vertices, wasting memory.

Instead, use a larger convex object. Wikipedia lists the coordinates for an icosahedron here. Expanded, the twelve vertices are

 0.000000  1.000000  1.618034
 0.000000  1.000000 -1.618034
 0.000000 -1.000000 -1.618034
 0.000000 -1.000000  1.618034
 1.000000  1.618034  0.000000
 1.000000 -1.618034  0.000000
-1.000000  1.618034  0.000000
-1.000000 -1.618034  0.000000
 1.618034  0.000000  1.000000
 1.618034  0.000000 -1.000000
-1.618034  0.000000  1.000000
-1.618034  0.000000 -1.000000

This generates the 20 half-spaces, and if you use t=4 or t=8, you should have a good test case.

I have the illusion (?), that the root holds the hyperplanes of the whole polytope and then we divide the polytope to 8 boxes, with every one box containing a part of the polytope

No illusion, that is completely correct.

The root holds all the hyperplanes of the polytope, because the root holds the entire polytope and thus must contain all of the halfspaces too.

So, you split the root cell into 2^D boxes (D being the number of dimensions), in a way that at least some of these boxes have fewer hyperplanes.

You continue splitting each box into sub-boxes, until there are at most t hyperplanes in that box.

Moreover, a technical question:
How am I suppose to insert new data in an octree?

You don't need a generic insert routine, because you insert the nodes as you descend into the tree. Insert happens automatically, when you construct the tree.

In pseudocode:

construct(halfspace_array, halfspace_count,
          vertex_array, vertex_count,
          box_minimum, box_maximum):

    Count the number of halfspaces in halfspace_array
    that have planes intersecting the current box;
    save this count as halfspace_inside.
    Move these to the start of the array, and the others
    to the end of the array.

    If halfspace_inside < t:
        This is a leaf node, containing the halfspace_inside
        halfspaces. Store them, and return.
        (Note: this could be in some cases be an empty node,
         with no intersecting planes. Then you need to decide
         whether this is inside or outside -- just test the final
         halfspaces in halfspace_array to find out.)

    Count the number of vertices in vertex_array
    that lie within the current box. Move vertices outside
    the current box, or at a corner of the current box,
    to the end of the list; save the count of the others as
    vertex_inside.

    If there is a vertex near the center of the box,
    or anywhere inside the box but not in the center,
    pick that as the split point. Otherwise pick the center
    of the box.

    Call construct() for each of the new boxes thus generated.
    (Remember to update the box_minimum and/or
     box_maximum to reflect the box for each new sub-box.)

[in a binary tree] it has to go left or right, but in the octree?

Whenever you descend into a tree, you have a value you compare to the value on each node, to decide which way you go.

In a binary tree, the value of the node determines whether to go left or right. If you have (a), and the node has (d), then

a <  d     ->next[0] == ->left
a >= d     ->next[1] == ->right

In a quadtree, you compare two-component vectors -- like (row,column) coordinates for example --, and you have 2²=4 possible child nodes. If the vector you have is say (a,b), and the node has (d,e), then the node you pick is

a <  d && b <  e        ->next[0]
a >= d && b <  e        ->next[1]
a <  d && b >= e        ->next[2]
a >= d && b >= e        ->next[3]

You can think of the node having a coordinate point as its value; the quadtree picks the correct quadrant compared to the node point, based on the point you are holding.

For an octree, you have (a,b,c) and the node has (d,e,f), and you pick the node

a <  d && b <  e && c <  f     ->next[0]
a >= d && b <  e && c <  f     ->next[1]
a <  d && b >= e && c <  f     ->next[2]
a >= d && b >= e && c <  f     ->next[3]
a <  d && b <  e && c >= f     ->next[4]
a >= d && b <  e && c >= f     ->next[5]
a <  d && b >= e && c >= f     ->next[6]
a >= d && b >= e && c >= f     ->next[7]

Similarly, you can think of the node having a 3D point as its value, and the octree picking the correct octant compared to the node point, based on the point you are holding.

In any 2^D-tree, you compare D values. There are 2^D possible results for the comparisons, and each result combination refers to a specific child pointer. Thus, you have 2^D child pointers per node.

The way I do it is construct the child index i so that where (have_value_k >= node_value_k) is true, k = 0..D-1, bit k is set in i. This yields the exact same results as the tables above.

[std10093]

The pseudocode is actually the method described in post #23.

>Count the number of vertices in vertex_array
that lie within the current box.

I am going to count the vertices that lie on the faces of the polytope too, I guess you agree with this.

However, it is not obvious where vertices_inside will be used... :/