Thread: Quadtrees

  1. #31
    Ticked and off
    Join Date
    Oct 2011
    Location
    La-la land
    Posts
    1,728
    So, just out of curiosity, I wrote my own version. (It's not perfect; I need to review the edge cases, as I get spurious wrong results when the test point is exactly at the edge of the bounding box).

    It turns out that just halving the bounding box yields pretty good results. (I still suspect that when there is only one vertex within the bounding box, it is best to pick that vertex as the splitting point -- just in case that vertex participates in a large number of sides.)

    For the few test cases I did, ignoring the vertices completely (after building the half-spaces) yielded very good results for regular convex polyhedra; it also worked acceptably for a 100-vertex polyhedron with all vertices randomly on the unit sphere.

    In other words, you probably should drop the vertices as soon as the halfspaces are generated. The vertex information does not seem to help as much as I originally thought.

    It might be a better approach to save the bounding box of the polytope side, when the halfspaces are generated, instead. (The bounding box is just the minimum and maximum coordinates of the vertices that are at the plane of the halfspace, since the polytope is convex.)

    This does make the size of each halfspace definition much bigger (from D+1 doubles to 3D+1 doubles, where D is the number of dimensions):
    Code:
    typedef struct {
        point_t min;
        point_t max;
        point_t normal;
        double distance;
    } halfspace_t;
    However, it might allow a better logic of how to divide the box efficiently.

  2. #32
    SAMARAS std10093's Avatar
    Join Date
    Jan 2011
    Location
    Nice, France
    Posts
    2,694
    Thank you for the 700th post.

    Well, I want to first make it work and then improve it.

    Sorry for asking again, but now I am a bit unsure. Thanks for the great help so far and if you think you are bored or something, just do not answer (I won't get angry at all, since you have helped² me).

    I understood how to find the min and max box of a polytope, but now, you are talking about the min and max box of the halfspace. Maybe I am confused, because of the simplicity!

    An example will help me get it. Assume the vertices I had before, the halfspaces are:
    Code:
    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.57735 0.57735 0.57735 0.57735 // This is sqrt(1/3)
    So, now, which are the bounding boxes?
    Code - functions and small libraries I use


    It’s 2014 and I still use printf() for debugging.


    "Programs must be written for people to read, and only incidentally for machines to execute. " —Harold Abelson

  3. #33
    Ticked and off
    Join Date
    Oct 2011
    Location
    La-la land
    Posts
    1,728
    Quote Originally Posted by std10093 View Post
    I understood how to find the min and max box of a polytope, but now, you are talking about the min and max box of the halfspace.
    I phrased it badly. A halfspace is infinite; it splits the entire space into two halves, and therefore it cannot have a "bounding box".

    I shall try to explain what I tried to convey, but my tests indicate this idea didn't work that well either. At this point, considering memory and CPU use and code complexity, I believe it's best to just split the box into equal-sized sub-boxes when recursively generating the tree.

    When the tree structure is constructed, the way you split each box into 2D sub-boxes in D dimensions, determines how efficient/shallow or inefficient/deep your tree structure becomes.

    Since my original idea of using the vertices to decide how to split the boxes when recursively generating the octree did not work (as well as I'd thought), I was looking for an alternative approach.

    The idea I described last was to record the bounding box of the vertices that define each face, in each halfspace structure, and use these to decide how to split the boxes when the tree is constructed.

    Consider the tetrahedron
    Code:
    0 0 0
    1 0 0
    0 1 0
    0 0 1
    The faces, and the bounding box per face I was talking about, would be
    Code:
    Vertex  Vertex  Vertex   Halfspace  Bounding box
    0 0 0   1 0 0   0 1 0    0 0 -1 0   (0,0,0)-(1,1,0)
    0 0 0   1 0 0   0 0 1    0 -1 0 0   (0,0,0)-(1,0,1)
    0 0 0   0 1 0   0 0 1    -1 0 0 0   (0,0,0)-(0,1,1)
    1 0 0   0 1 0   0 0 1    √⅓√⅓√⅓√⅓   (0,0,0)-(1,1,1)
    These "bounding boxes" only tell us the bounding box of the original face that was used to construct the halfspace, i.e. the bounding box of the vertices that defined the face. They can be used to find regions of high contention (i.e. where it should be beneficial to split the space).

    However, the typical cases for convex polytopes, doing that didn't yield significantly better (shallower) trees. Definitely not better enough to offset the memory use.



    To be perfectly clear, my current structures and constructor functions look something like this:
    Code:
    typedef struct {
        double x;
        double y;
        double z;
    } point_t;
    
    typedef struct {
        double x;
        double y;
        double z;
        double d;
    } halfspace_t;
    
    typedef struct node  node_t;
    struct node {
        int          halfspaces; /* 0 */
        point_t      split;
        struct node *child[8];
    };
    
    typedef struct {
        int          halfspaces; /* >0 */
        halfspace_t  halfspace[HALFSPACES_MAX];
    } leaf_t;
    
    /* If the entire box is within the polyhedron */
    typedef struct {
        int          halfspaces; /* < 0 */
    } inside_leaf_t;
    
    typedef struct {
        point_t min;
        point_t max;
        node_t *tree;
    } convex_polyhedron_t;
    
    /* Recursive function to construct halfspace octrees.
     * Function may reorder the first nh halfspaces in h.
    */
    node_t *convex_polyhedron_node(const point_t min, const point_t max,
                                   halfspace_t *const h, const size_t nh);
    
    /* Construct a convex polyhedron based on halfspaces.
     * bbmin,bbmax define the bounding box for the polyhedron.
     * Return 0 if successful, errno otherwise.
    */
    int convex_polyhedron(convex_polyhedron_t *const p,
                          const point_t bbmin, const point_t bbmax,
                          halfspace_t *const h, const size_t nh);
    By the way, I'm using dot and circo from Graphviz to visualize the above convex_polyhedron_t structures. I have a C function which outputs the dot-language descriptions, which I can visualize using dot -Tx11 or circo -Tx11, or save to e.g. PDF (-Tpdf) or SVG (-Tsvg) files.

  4. #34
    SAMARAS std10093's Avatar
    Join Date
    Jan 2011
    Location
    Nice, France
    Posts
    2,694
    >By the way, I'm using dot and circo from Graphviz to visualize the above convex_polyhedron_t structures. I have a C function which outputs the dot-language descriptions, which I can visualize using dot -Tx11 or circo -Tx11, or save to e.g. PDF (-Tpdf) or SVG (-Tsvg) files.

    That's seems to be interesting. Very interesting! I have written the program in C++, but I do not think is a problem, is it? Moreover, I have no idea what this is, so would it be hard to try it? What should I do?



    As for the algorithm, I now have digest how the octree is going to be "automatically" be built, thanks to you, thus I am thinking various ways. So, let's say I split the original polytope somehow, the question I have, is why to divide based on the halfspaces?

    Why not to split, based on the vertices of the polytope? That way it seems easier to me, in the paper at least and more intuitive.

    One drawback I just thought, is that assuming we split the polytope by vertices, then in the leaves of the octree, we are going to have the vertices of the subpolytopes, but we need to generate the halfspaces, thus we have to pay the cost of the generation more than once (one per leaf).
    However, if we split based on the halfspaces, we pay only the cost of the split, since the halfspaces are ready.

    To make myself clear, I want to compare the following two approaches.
    First method - divide by vertices.
    • Split the original polytope, based on the vertices.
    • The leaves contain the vertices of the subpolytopes.
    • For every leaf, generate the halfspaces of the corresponding subpolytope.


    Second method - divide by halfspaces.
    • Generate the halfspaces of the original polytope.
    • Split the original polytope, based on the halfspaces.
    • The leaves contain the halfspaces of the subpolytopes.
    Code - functions and small libraries I use


    It’s 2014 and I still use printf() for debugging.


    "Programs must be written for people to read, and only incidentally for machines to execute. " —Harold Abelson

  5. #35
    Ticked and off
    Join Date
    Oct 2011
    Location
    La-la land
    Posts
    1,728
    Quote Originally Posted by std10093 View Post
    That's seems to be interesting. Very interesting! I have written the program in C++, but I do not think is a problem, is it? Moreover, I have no idea what this is, so would it be hard to try it? What should I do?
    No problem at all! Consider this small example program. It creates a small binary tree, and then spits out the definition in dot (language) to standard output. I'm using the pointers (node addresses) as the identifiers:
    Code:
    #include <stdlib.h>
    #include <string.h>
    #include <stdio.h>
    #include <errno.h>
    
    struct node {
        struct node *left;  /* Less than or equal to */
        struct node *right; /* Greater than */
        int          value;
    };
    
    struct node *tree_add(struct node **const rootptr, int value)
    {
        struct node *newnode, *curr, **pptr;
    
        /* Tree not specified? */
        if (!rootptr) {
            errno = EINVAL;
            return NULL;
        }
    
        /* Create node. */
        newnode = malloc(sizeof *newnode);
        if (!newnode) {
            errno = ENOMEM;
            return NULL;
        }
    
        /* Initialize node. */
        newnode->left = NULL;
        newnode->right = NULL;
        newnode->value = value;
    
        /* Descend into the tree. */
        pptr = rootptr;  /* Pointer to pointer used to get to curr. */
        curr = *rootptr; /* Current node. */
        while (curr)
            if (curr->value >= value) {
                pptr = &curr->left;
                curr = curr->left;
            } else {
                pptr = &curr->right;
                curr = curr->right;
            }
    
        /* Add node to where we dropped off the tree. */
        *pptr = newnode;
    
        /* Done! */
        return newnode;
    }
    
    /* Recursively describe the nodes of the tree in DOT.
    */
    void fdot_node(FILE *const out, const struct node *const node)
    {
        fprintf(out, "\t\"%p\" [ label=\"%d\", shape=oval ]\n", node, node->value);
    
        if (node->left) {
            fdot_node(out, node->left);
            fprintf(out, "\t\"%p\" -> \"%p\" [ taillabel=\"≤\" ]\n", node, node->left);
        }
    
        if (node->right) {
            fdot_node(out, node->right);
            fprintf(out, "\t\"%p\" -> \"%p\" [ taillabel=\">\" ]\n", node, node->right);
        }
    
        return;
    }
    
    /* Write the structure of the tree in DOT language.
     * Returns 0 if success, errno error code otherwise.
    */
    int fdot(FILE *const out, const struct node *const tree)
    {
        /* Invalid parameters? */
        if (!out || !tree)
            return errno = EINVAL;
    
        /* I/O error? Paranoid check, really.. */
        if (ferror(out))
            return errno = EIO;
    
        /* Start by describing the entire graph. */
        fprintf(out, "digraph {\n");
    
        /* Describe nodes recursively. */
        fdot_node(out, tree);
    
        /* Finish the graph. */
        fprintf(out, "}\n\n");
    
        /* Still everything OK in out? */
        if (ferror(out))
            return errno = EIO;
    
        /* Success! */
        return errno = 0;
    }
    
    int main(void)
    {
        struct node *tree = NULL;
    
        tree_add(&tree, 5);
        tree_add(&tree, 3);
        tree_add(&tree, 7);
        tree_add(&tree, 2);
        tree_add(&tree, 4);
        tree_add(&tree, 6);
        tree_add(&tree, 8);
        tree_add(&tree, 1);
        tree_add(&tree, 9);
    
        fdot(stdout, tree);
    
        return 0;
    }
    If you compile and run it, it generates something like
    Code:
    digraph {
    	"0x22b9010" [ label="5", shape=oval ]
    	"0x22b9030" [ label="3", shape=oval ]
    	"0x22b9070" [ label="2", shape=oval ]
    	"0x22b90f0" [ label="1", shape=oval ]
    	"0x22b9070" -> "0x22b90f0" [ taillabel="≤" ]
    	"0x22b9030" -> "0x22b9070" [ taillabel="≤" ]
    	"0x22b9090" [ label="4", shape=oval ]
    	"0x22b9030" -> "0x22b9090" [ taillabel=">" ]
    	"0x22b9010" -> "0x22b9030" [ taillabel="≤" ]
    	"0x22b9050" [ label="7", shape=oval ]
    	"0x22b90b0" [ label="6", shape=oval ]
    	"0x22b9050" -> "0x22b90b0" [ taillabel="≤" ]
    	"0x22b90d0" [ label="8", shape=oval ]
    	"0x22b9110" [ label="9", shape=oval ]
    	"0x22b90d0" -> "0x22b9110" [ taillabel=">" ]
    	"0x22b9050" -> "0x22b90d0" [ taillabel=">" ]
    	"0x22b9010" -> "0x22b9050" [ taillabel=">" ]
    }
    except that the pointer labels ("0x...") will vary, of course. If you save that as say example.dot, and then run
    Code:
    dot -Tpng example.dot > example.png
    the end result will look like this:
    Attachment 12845
    Note that dot does reorder the edges as it sees best; in particular, any edge starting from the left side from an oval is not necessarily the left child. That's why I added the tail labels to the arrows.

    Using
    Code:
    neato -Tpng example.dot > neato.png
    from the same data will generate an image like
    Attachment 12844

    If you generate SVG (-Tsvg) images, you can edit them in e.g. Inkscape, letting you fine-tune the images for e.g. your documentation.



    Quote Originally Posted by std10093 View Post
    So, let's say I split the original polytope somehow, the question I have, is why to divide based on the halfspaces? Why not to split, based on the vertices of the polytope?
    Consider a half-sphere in three or more dimensions. It is convex, and perfectly suitable for this algorithm. The flat side is obviously only one halfspace.

    If you split more than once, none of the vertices defining the flat face will be inside your boxes. However, some of the boxes will still intersect with the flat face, and thus must be limited by the flat face halfspace.

    Because of this, you cannot consider only those vertices that lie within current box, but must consider all vertices that take part in the faces that intersect with the current box. Because halfspace == face, it's just more efficient to use the halfspaces.

    If you want to use the vertices, you can, but then you must somehow keep all the necessary vertices along when you recursively generate the tree. You cannot just drop any vertex that happens to be outside the current box, because such vertices can still participate in the faces that intersect the current box.

    Quote Originally Posted by std10093 View Post
    The leaves contain the halfspaces of the subpolytopes.
    Note that in my method, I did not construct subpolytopes. My logic was more like:

    • Start with the halfspaces forming the original convex polytope,
      and with the bounding box of the original convex polytope.
    • If the current box intersects with not more than T halfspace planes (face hyperplanes), then this node in the tree is a leaf, and specifies these halfspaces.
      Otherwise:
      Split the current box into 2D sub-boxes, and repeat this step for each of them.

    Theoretically, it is possible that the current box is completely inside or completely outside the polytope. So, in my program I also checked for that, and used NULL for completely outside, and a special "all inside" leaf node for completely inside.

    This procedure means that within each box, we only consider those halfspaces of the original polytope that affect points within each box. Those halfspaces do not form a polytope, they are open to infinity in some directions. I don't know how to explain this well, but it is clear if you look at the nodes.

    Perhaps it would be best if I clean up my implementation and post it here?

    I'm pretty sure at this point that it won't affect you unduly, because you're certainly thinking about this at least as much as I have; but, it would allow us to look at my implementation and make it easier to discuss the choices I made, and why I made them. In particular, I'm not at all certain I fully followed the article, so relying on anything I've said or have done is not necessarily correct.

  6. #36
    SAMARAS std10093's Avatar
    Join Date
    Jan 2011
    Location
    Nice, France
    Posts
    2,694
    For the graphics, it seems that you assumed I have Linux and I still haven't installed them on my machine. Every time I wanted to do something on Linux, I used the machines of my university and I do not think they have this program installed.

    I had understood your thinking, but what I presented was my thinking. I believe, that with my thinking, every node of the octree will have something from the origin polytope, thus we won't waste any memory at all. However, I will need some space for the new subpolytopes, since they may have many in common with the original polytope, but they will have something new. Your idea, presented a way that would make use of only one array and somehow move data from one index to another. Do you mind explaining this one more time?

    Please, wait until I have something ready and then we post our implementations, because even though we have spent -pretty much- the same amount of thinking in the algorithm, you are far more experienced than me.

    // I really try to learn and you are a BIG help for me
    Code - functions and small libraries I use


    It’s 2014 and I still use printf() for debugging.


    "Programs must be written for people to read, and only incidentally for machines to execute. " —Harold Abelson

  7. #37
    Ticked and off
    Join Date
    Oct 2011
    Location
    La-la land
    Posts
    1,728
    Quote Originally Posted by std10093 View Post
    For the graphics, it seems that you assumed I have Linux
    No, Graphviz is available on all major OSes, including Windows and Mac OS X, not just Linux.

    Quote Originally Posted by std10093 View Post
    I had understood your thinking, but what I presented was my thinking. I believe, that with my thinking, every node of the octree will have something from the origin polytope, thus we won't waste any memory at all.
    Ah. I obviously didn't quite catch your idea.

    As to the memory use, my implementation uses the full array of halfspaces only during tree construction; after tree construction, the array is released. My tests indicate that this is actually quite frugal wrt. memory use.

    Quote Originally Posted by std10093 View Post
    Please, wait until I have something ready and then we post our implementations
    Will do.

    Quote Originally Posted by std10093 View Post
    Your idea, presented a way that would make use of only one array and somehow move data from one index to another. Do you mind explaining this one more time?
    I'd be delighted to, of course

    Let's start with a cube with vertices at ±1, ±1, ±1. The six half-spaces are
    Code:
    +1 0 0 +1     (the side towards +x)
    0 +1 0 +1     (the side towards +y)
    0 0 +1 +1     (the side towards +z)
    0 0 -1 +1     (the side towards -z)
    0 -1 0 +1     (the side towards -y)
    -1 0 0 +1     (the side towards -x)
    The bounding box is something like (ϵ-1,ϵ-1,ϵ-1)-(ϵ+1,ϵ+1,ϵ+1), where ϵ depends on how you test for inclusion in the bounding box.

    In my case, I use the rule xminx < xmax, in which case the bounding box values in C are
    Code:
        min.x = -1.0;
        min.y = -1.0;
        min.z = -1.0;
        max.x = nextafter(1.0, HUGE_VAL);
        max.y = nextafter(1.0, HUGE_VAL);
        max.z = nextafter(1.0, HUGE_VAL);
    i.e. max are the smallest values larger than the vertex coordinates. In practical terms, it means that in my logic, min is inclusive (inside the box), but max exclusive (outside).

    (If you wonder why I bother with this detail, not bothering with this was what caused the problems I mentioned I had with my code. One has to be unambiguous with the corner cases, or problems will ensue!)

    Okay, we're ready to construct the tree, with initial halfspace list h containing nh=6 halfspaces. Note that nh is local to each function call; if it is changed, the changes are only visible in any recursive calls that function makes, not in any parent calls.

    The initial call checks all halfspaces, and sees that all their planes intersect with the initial box (-1 to +1 in all coordinates). Therefore, it splits the space into 23 = 8 sub-boxes, and recursively calls each of them.

    The first sub-box contains (-1,-1,-1)-(0,0,0). Only threeNOTE 1 of the six halfspaces have planes that intersect with this sub-box, so the halfspaces in h get rearranged into
    Code:
    0 0 -1 +1     (the side towards -z)
    0 -1 0 +1     (the side towards -y)
    -1 0 0 +1     (the side towards -x)
    new nh = 3
    +1 0 0 +1     (the side towards +x)
    0 +1 0 +1     (the side towards +y)
    0 0 +1 +1     (the side towards +z)
    If the number of halfspaces per leaf is smaller than three, then this would get subdivided further, with the first recursive child getting box (-1,-1,-1)-(-0.5,-0.5,-0.5), and see only the new nh=3 halfspaces at the start of the array. It is free to rearrange those; it won't change anything the parent cares. (It does not matter in which order we have the halfspaces.)

    If the number of halfspaces per leaf is three or more, then these three halfspaces form a leaf node in the tree. The function allocates a suitable leaf node, copies the halfspace specifications, and returns to the caller. Note: you cannot just refer to the original halfspace array using e.g. a pointer, because the order of the halfspaces in the array changes.

    Back in the first call, we have only handled the first of the eight sub-boxes, and already the contents of the h array have been rearranged. The second sub-box handles (0,-1,-1)-(1,0,0), and rearranges the halfspace list left from the first sub-box,
    Code:
    0 0 -1 +1     (the side towards -z)
    0 -1 0 +1     (the side towards -y)
    -1 0 0 +1     (the side towards -x)
    +1 0 0 +1     (the side towards +x)
    0 +1 0 +1     (the side towards +y)
    0 0 +1 +1     (the side towards +z)
    nh = 6
    into
    Code:
    +1 0 0 +1     (the side towards +x)
    0 0 -1 +1     (the side towards -z)
    0 -1 0 +1     (the side towards -y)
    new nh = 3
    -1 0 0 +1     (the side towards -x)
    0 +1 0 +1     (the side towards +y)
    0 0 +1 +1     (the side towards +z)
    Remember, changes in nh are not visible to the parent call, only in any recursive calls the function might make.

    Because the first thing each call does, is rearrange the array, it is okay for any call to rearrange the array. Simply put, as long as the array each call sees still contains the exact same halfspaces, it is perfectly ok for any recursive call to reorder the halfspaces: we just do not rely at all on the order of the halfspaces.

    This repeats for the six other recursive sub-box calls.

    After the initial call has done the recursive calls for all eight sub-boxes, it is done, and returns. (My function calls return a pointer to the node they generate; it is a very simple and common approach. In this case the initial call has generated an inner node with eight child pointers.)

    Since each call that generated new leaf nodes copied the halfspaces (and inner nodes do not contain any halfspace definitions), the original halfspace array h is no longer needed, and can be freed.

    Overall, in addition to the memory required by the tree, and whatever stack space is temporarily needed by the recursive function calls, the memory overhead is just the original halfspace array h. For my 3D implementation, that is four doubles per face, and I consider it very frugal.


    1 Actually, this too depends slightly on how you handle the corner cases wrt. the halfspaces.
    In my case, I am assuming the end result is to create a fast polytope inclusion test. Therefore, I exclude all halfspaces that do not affect the inclusion test within the current box.
    If all eight corners of the current box are outside any half-space planes, the current box is completely outside the polytope (and I denote such leaves with NULL pointers).
    When all eight corners of the current box are either inside of or on the half-space plane of the halfspace, that halfspace is not needed to make the decision.
    Because I'm only interested in halfspaces where at least one corner of the current box is inside, and at least one corner outside the halfspace, I only get three halfspaces, not six.

  8. #38
    SAMARAS std10093's Avatar
    Join Date
    Jan 2011
    Location
    Nice, France
    Posts
    2,694
    We could set
    Code:
    max.x = nextafter(1.0, HUGE_VAL);
    as 1.0 + ε or 1.0 + ε², since we calculate the squared epsilon at some point (thus we do not pay the cost of the computation again). Of course, we could apply this to all maxes, not just the one of x.
    // Well actually I know see, that you imply this too, aren't you?

    So, you apply your logic and for the inner nodes, you only manipulate the original array. Then, for every leaf, you actually create a subarray, that contains some of the data of the main array, not new data, but only a part of the main array. I think I have an idea on this, but I will first have to give it a try..

    Finally, the only thing that seems to me a little obscure, is how do you divide the box and generate the sub-boxes? What's your logic behind this?
    Code - functions and small libraries I use


    It’s 2014 and I still use printf() for debugging.


    "Programs must be written for people to read, and only incidentally for machines to execute. " —Harold Abelson

  9. #39
    Ticked and off
    Join Date
    Oct 2011
    Location
    La-la land
    Posts
    1,728
    Quote Originally Posted by std10093 View Post
    We could set
    Code:
    max.x = nextafter(1.0, HUGE_VAL);
    Exactly. This is very important, if min is inclusive, but max exclusive.

    In other words, for the initial bounding box, each max coordinate must be greater than any of the corresponding vertex coordinates. I find the maximum vertex coordinate, and then use max.c = nextafter(vertex_max.c, HUGE_VAL); (but min.c = vertex_min.c).

    Quote Originally Posted by std10093 View Post
    for the inner nodes, you only manipulate the original array. Then, for every leaf, you actually create a subarray, that contains some of the data of the main array
    Yes, exactly.

    Quote Originally Posted by std10093 View Post
    Finally, the only thing that seems to me a little obscure, is how do you divide the box and generate the sub-boxes? What's your logic behind this?
    I didn't see any details on how to split the box in the article you linked to (page 580, SplitReduce function), so I investigated a couple of possibilities. I've described them to pretty good detail in this thread, but they both were duds.

    For now, the best option seems to simply split it in half.

    In practice: Let's say the current box has been deemed to be an inner node, and has to be split. The current bounding box is (min.x, min.y, min.z) - (max.x, max.y, max.z), max exclusive. I calculate their mean,
    Code:
    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);
    Note that I do store the split in the inner node, so that for point inclusion tests, it is easy to compare to and decide which child node to descend into.

    The child boxes will be
    Code:
    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, maxt.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)
    so that for point (x, y, z) inclusion test, the child to descend into is
    Code:
    = 0
    + 1 if x >= split.x
    + 2 if y >= split.y
    + 4 if z >= split.z


    The different methods of splitting the boxes affect only the split coordinates in each step. It must still be within the min and max, but a very good choice will create more efficient (shallower) trees.

    The third method I've been thinking about is to consider the coordinates at which each halfspace intersects the current box.

    If we ignore halfspaces coplanar to an axis, each halfspace opens to either positive or negative infinity along each axis. We can compute the maximum coordinate within the current box for each plane that opens to negative infinity, and the minimum coordinate within the current box for each plane that opens to positive infinity. If these two are sorted (in separate arrays, not mixed), comparing the two arrays tells us how the halfspaces "overlap" within the current box.

    If there is a vertex in the current box, and lots of faces (radiating from it), the maximum "overlap" should be pretty narrowly near the vertex. So, we could try splitting at the maximum "overlap" region. (Perhaps a sanity rule, so that the larger side of the split is never larger than say eight times the smaller side, might be in order too.)

    Does this work? I have no idea. I haven't worked out the math, but looking at some of the structure trees I've generated, there seems to be room for improvement in how the boxes are split. The two methods I've tried thus far, just did not improve from the always-split-in-half method.

  10. #40
    SAMARAS std10093's Avatar
    Join Date
    Jan 2011
    Location
    Nice, France
    Posts
    2,694
    Well, first I want to try this and make it work.

    You have a typo here
    Code:
    split.z = min.y + 0.5*(max.z - min.z);
    It needs to be z, not y.

    Yes, the article does not provide any information on how to split the boxes, because it maybe implies that we should just halve them.

    I am pretty sure I can start building the octree, but I have a problem in feeling the sub-boxes and I really want to learn, not just make it work.
    Quote Originally Posted by Nominal Animal View Post
    Let's start with a cube with vertices at ±1, ±1, ±1.
    The first sub-box contains (-1,-1,-1)-(0,0,0). Only threeNOTE 1 of the six halfspaces have planes that intersect with this sub-box, so the halfspaces in h get rearranged into
    Code:
    0 0 -1 +1     (the side towards -z)
    0 -1 0 +1     (the side towards -y)
    -1 0 0 +1     (the side towards -x)
    new nh = 3
    +1 0 0 +1     (the side towards +x)
    0 +1 0 +1     (the side towards +y)
    0 0 +1 +1     (the side towards +z)
    I think I can understand this.
    The way I visualize it, is like this:
    Quadtrees-unitcube-jpg
    So, the plane of x and y axes cut the cub in the middle. Then, because z is vertical to the page, 8 sub-cubes are there. One per quadrant, on the positive z axis, thus 4 sub-cubes and -likewise- 4 sub-cubes on the negative z axis, thus 8 in the whole 3D space.

    The sub-box containing (-1,-1,-1)-(0,0,0) is actually the sub-cube in the 3rd quadrant (left and down in the x-y plane in the image), which lies below the x-y plane.
    Quote Originally Posted by Nominal Animal View Post
    The second sub-box handles (0,-1,-1)-(1,0,0), and rearranges the halfspace list left from the first sub-box into
    Code:
    +1 0 0 +1     (the side towards +x)
    0 0 -1 +1     (the side towards -z)
    0 -1 0 +1     (the side towards -y)
    new nh = 3
    -1 0 0 +1     (the side towards -x)
    0 +1 0 +1     (the side towards +y)
    0 0 +1 +1     (the side towards +z)
    Now, with the same logic as above, I would say that this sub-box is actually the sub-cube in the 4th quadrant (right and down in the image), again lying below the x-y plane.

    However, I would expect the array to be rearranged into this:
    Code:
    -1 0 0 +1     (the side towards -x)
    0 0 -1 +1     (the side towards -z)
    0 +1 0 +1     (the side towards +y)
    new nh = 3
    +1 0 0 +1     (the side towards +x)
    0 -1 0 +1     (the side towards -y)
    0 0 +1 +1     (the side towards +z)
    So, where have I lost it? I googled a lot, but could not find something that would make my visualize the situation with the sub-boxes and the halfspaces.
    Code - functions and small libraries I use


    It’s 2014 and I still use printf() for debugging.


    "Programs must be written for people to read, and only incidentally for machines to execute. " —Harold Abelson

  11. #41
    Ticked and off
    Join Date
    Oct 2011
    Location
    La-la land
    Posts
    1,728
    Quote Originally Posted by std10093 View Post
    You have a typo
    Good catch. I can't edit the post anymore, so I'll have to leave it there.

    Quote Originally Posted by std10093 View Post
    The way I visualize it, is like this
    Yes, that's correct for the initial box.

    Quote Originally Posted by std10093 View Post
    The sub-box containing (-1,-1,-1)-(0,0,0) is actually the sub-cube in the 3rd quadrant (left and down in the x-y plane in the image), which lies below the x-y plane.
    There is no standard numbering for the quadrants (or octants, rather); you are completely free to devise your own numbering. The numbering I use is the simplest one in my opinion, as choosing the correct octant only requires D comparisons for D dimensions, addition, and bit shifts.

    Quote Originally Posted by std10093 View Post
    However, I would expect the array [for the (0,-1,-1)-(1,0,0) sub-box] to be rearranged into this:
    No. Why would you expect to see any other sides/halfspaces/planes than x=1 (+x), y=-1 (-y), and z=-1 (-z)?

    Just to be sure, I checked. My code gives me
    Code:
    0: -x -y -z    Halfspaces
                    0  0 -1 1   (z = -1 plane)
                    0 -1  0 1   (y = -1 plane)
                   -1  0  0 1   (x = -1 plane)
    
    1: +x -y -z    Halfspaces
                    0  0 -1 1   (z = -1 plane)
                    0 -1  0 1   (y = -1 plane)
                   +1  0  0 1   (x = +1 plane)
    
    2: -x +y -z    Halfspaces
                    0  0 -1 1   (z = -1 plane)
                   -1  0  0 1   (x = -1 plane)
                    0 +1  0 1   (y = +1 plane)
    
    3: +x +y -z    Halfspaces
                    0  0 -1 1   (z = -1 plane)
                   +1  0  0 1   (x = +1 plane)
                    0 +1  0 1   (y = +1 plane)
    
    4: -x -y +z    Halfspaces
                   -1  0  0 1   (x = -1 plane)
                    0 -1  0 1   (y = -1 plane)
                    0  0 +1 1   (z = +1 plane)
    
    5: +x -y +z    Halfspaces
                   +1  0  0 1   (x = +1 plane)
                    0 -1  0 1   (y = -1 plane)
                    0  0 +1 1   (z = +1 plane)
    
    6: -x +y +y    Halfspaces
                   -1  0  0 1   (x = -1 plane)
                    0 +1  0 1   (y = +1 plane)
                    0  0 +1 1   (z = +1 plane)
    
    7: +x +y +z    Halfspaces
                   +1  0  0 1   (x = +1 plane)
                    0 +1  0 1   (y = +1 plane)
                    0  0 +1 1   (z = +1 plane)
    where
    Code:
    -x means x = [ -1, 0 ),
    +x means x = [ 0, +1 ),
    -y means y = [ -1, 0 )
    +y means y = [ 0, +1 ),
    -z means z = [ -1, 0 ), and
    +z means z = [ 0, +1 ).
    and the sub-boxes are recursively called in the above order. I used a very slightly rotated cube (vertices not exactly at ±1) to make sure, then rounded the coordinates.

    You might also wish to draw the boxes in some orthographic projection. My preferred one is
    Attachment 12850
    where X is full right and half up, Y is full up and half left, and Z is full up; origin (0,0,0) is at the center of the image.

  12. #42
    SAMARAS std10093's Avatar
    Join Date
    Jan 2011
    Location
    Nice, France
    Posts
    2,694
    I think I understood, thanks.

    However, I spotted something. I had the impression, that the number of vertices is the same as the number of halfspaces. It happened before, with one tetrahedron we had and I asked and you replied affirmatively, but you probably meant for this example, that the number is the same.

    For the unit cube for instance, we have 8 vertices and 6 halfspaces.

    So, after some thinking, I would say, that the number of halfspaces is, at most, equal to the number of vertices. What do you say?
    Code - functions and small libraries I use


    It’s 2014 and I still use printf() for debugging.


    "Programs must be written for people to read, and only incidentally for machines to execute. " —Harold Abelson

  13. #43
    Ticked and off
    Join Date
    Oct 2011
    Location
    La-la land
    Posts
    1,728
    Quote Originally Posted by std10093 View Post
    I had the impression, that the number of vertices is the same as the number of halfspaces.
    Only for some specific polytopes.

    Quote Originally Posted by std10093 View Post
    I would say, that the number of halfspaces is, at most, equal to the number of vertices. What do you say?
    Nope. Consider an octahedron, with six vertices:
    Code:
     0  0 -1
    +1 +1  0
    -1 +1  0
    -1 -1  0
    +1 -1  0
     0  0 +1
    Yet, it has eight faces, and therefore is defined by eight halfspaces.

    (I am a bit worried I might have stated somewhere that the number of halfspaces is at most the number of vertices. If so, it's my error, and I'm sorry for the confusion. It is an easy error to make, because it is easy to forget that each face and halfspace is derived from a combination of vertices.)

    In general, the Euler characteristic, χ = vertices - edges + faces, is 0 for even-dimensional convex polytopes, and 2 for odd-dimensional convex polytopes.

    There are other interesting things polyhedral combinatorics tells about convex polytopes. For example, the Upper Bound Conjecture says that asymptotically, the maximum number of faces a convex polytope of d dimensions and n vertices is
    O(nd/2⌋)

    Interesting stuff.

  14. #44
    SAMARAS std10093's Avatar
    Join Date
    Jan 2011
    Location
    Nice, France
    Posts
    2,694
    Really, very very interesting material.

    The splitReduce algorithm is dealing with convex polytopes.

    Is the octahedron a convex polytope? It has 6 vertices and 8 faces. It lies into the category of 3D, thus d is 3.

    The upper bound conjecture, says that the max number of faces is O(n^(⌊d/2⌋)), where d dimensions and n vertices.
    d = 3, thus d/2 = 1.5 and with the floor function this yields ⌊d/2⌋ = 1.
    n = 6, thus faces = 6^1 = 6. But the faces, are obviously 8.
    The only thing that can justify this, is that the octahedron is not convex, but it sure seems like a convex polytope to me.

    I would like to emphasize that the algorithm focuses only on convex polytopes. The question I have, is which is the maximum number of halfspaces, that n vertices can generate? Is the answer to this, the Upper Bound Conjecture?

    A quick Greek lesson.
    hedra means face.
    octo means eight, thus octo + hedra yields an octahedron (it has eight faces).
    Code - functions and small libraries I use


    It’s 2014 and I still use printf() for debugging.


    "Programs must be written for people to read, and only incidentally for machines to execute. " —Harold Abelson

  15. #45
    Ticked and off
    Join Date
    Oct 2011
    Location
    La-la land
    Posts
    1,728
    Quote Originally Posted by std10093 View Post
    Is the octahedron a convex polytope?
    Yes, because it can be described as the intersection of halfspaces. (Or, using the other definition, all lines between any two points inside the polytope are completely inside the polytope, so the polytope is convex.)

    You forgot how to read O() notation. Since it is asymptotic, it does not (necessarily) apply to small numbers of dimensions. Since it is just the order of, it might have a very large constant (that is not included in the notation).

    The maximum number of halfspaces in a convex polytope of D dimensions and N vertices is given by fD-1:
    Code:
    Odd D:
                  ⎛(D-1)/2 ⎛N-D-1+i⎞⎞
        faces ≤ 2 ⎜   ∑    ⎜       ⎟⎟
                  ⎝  i=0   ⎝   i   ⎠⎠
    
    Even D:
                  ⎛ D/2 ⎛N-D-1+i⎞⎞   ⎛N-D/2-1⎞
        faces ≤ 2 ⎜  ∑  ⎜       ⎟⎟ - ⎜       ⎟
                  ⎝ i=0 ⎝   i   ⎠⎠   ⎝  D/2  ⎠
    Assuming I calculated right, the maximum number of faces in a D-dimensional convex polytope is
    D = 3: faces ≤ 2 × vertices - 4
    D = 4: faces ≤ 1/2 × vertices2 - 3/2 × vertices
    D = 5: facesvertices2 - 7 × vertices + 12
    D = 6: faces ≤ 1/6 × vertices3 - 3/2 × vertices2 + 10/3 × vertices
    D = 7: faces ≤ 1/3 × vertices3 - vertices2 + 74/3 × vertices - 40
    D = 8: faces ≤ 1/24 × vertices4 - 3/4 × vertices3 + 107/24 × vertices2 - 35/4 vertices
    D = 9: faces ≤ 1/12 × vertices4 - 13/6 × vertices3 + 251/12 × vertices2 - 533/6 × vertices + 140
    and so on.

Popular pages Recent additions subscribe to a feed