Like Tree1Likes

Quadtrees

This is a discussion on Quadtrees within the Tech Board forums, part of the Community Boards category; Good catch, std10093. I'll investigate. Edit 1: There was a typo when defining v3 . The last component should be ...

  1. #16
    Registered User
    Join Date
    Oct 2011
    Posts
    890
    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).
    Code:
    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
    Code:
    /* 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 N3/6-N2/2+N/3 iterations (instead of N3).
    Code:
    /* 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.
    Last edited by Nominal Animal; 07-23-2013 at 12:15 PM.

  2. #17
    SAMARAS std10093's Avatar
    Join Date
    Jan 2011
    Location
    Nice, France
    Posts
    2,694
    I agree with all you said. I am playing around with the code and check this
    Code:
    -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..
    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. #18
    Registered User
    Join Date
    Oct 2011
    Posts
    890
    Quote Originally Posted by std10093 View Post
    I do not understand how we represent the halfspace.
    Consider a cube (-1,-1,-1)-(+1,+1,+1):
    Code:
    -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
    Code:
    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,
    Code:
    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
    Code:
    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 p0. Then, a point p is on the plane if and only if
    n·(p - p0) = 0
    n·p = n·p0
    The right side is exactly equivalent to our use of d, i.e. we use
    d = n·p0
    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.
    Last edited by Nominal Animal; 07-24-2013 at 10:48 AM. Reason: Fixed sign of n·p0

  4. #19
    SAMARAS std10093's Avatar
    Join Date
    Jan 2011
    Location
    Nice, France
    Posts
    2,694
    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
    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. #20
    SAMARAS std10093's Avatar
    Join Date
    Jan 2011
    Location
    Nice, France
    Posts
    2,694
    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.
    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

  6. #21
    Registered User
    Join Date
    Oct 2011
    Posts
    890
    Quote Originally Posted by std10093 View Post
    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(ni pi, 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).

    Quote Originally Posted by std10093 View Post
    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?

  7. #22
    SAMARAS std10093's Avatar
    Join Date
    Jan 2011
    Location
    Nice, France
    Posts
    2,694
    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.
    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

  8. #23
    Registered User
    Join Date
    Oct 2011
    Posts
    890
    Quote Originally Posted by std10093 View Post
    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 NPARENT 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 2D 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 NPARENT - 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.

  9. #24
    SAMARAS std10093's Avatar
    Join Date
    Jan 2011
    Location
    Nice, France
    Posts
    2,694
    Quote Originally Posted by Nominal Animal View Post
    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?
    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

  10. #25
    Registered User
    Join Date
    Oct 2011
    Posts
    890
    Quote Originally Posted by std10093 View Post
    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).

    Quote Originally Posted by std10093 View Post
    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
    Code:
    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
    Code:
    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.)

  11. #26
    SAMARAS std10093's Avatar
    Join Date
    Jan 2011
    Location
    Nice, France
    Posts
    2,694
    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
    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

  12. #27
    Registered User
    Join Date
    Oct 2011
    Posts
    890
    Quote Originally Posted by std10093 View Post
    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.

    Quote Originally Posted by std10093 View Post
    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.

  13. #28
    SAMARAS std10093's Avatar
    Join Date
    Jan 2011
    Location
    Nice, France
    Posts
    2,694
    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
    Code:
    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?
    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

  14. #29
    Registered User
    Join Date
    Oct 2011
    Posts
    890
    Quote Originally Posted by std10093 View Post
    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
    Code:
     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.

    Quote Originally Posted by std10093 View Post
    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 2D 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.

    Quote Originally Posted by std10093 View Post
    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:
    Code:
    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.)
    Quote Originally Posted by std10093 View Post
    [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
    Code:
    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 22=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
    Code:
    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
    Code:
    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 2D-tree, you compare D values. There are 2D possible results for the comparisons, and each result combination refers to a specific child pointer. Thus, you have 2D child pointers per node.

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

  15. #30
    SAMARAS std10093's Avatar
    Join Date
    Jan 2011
    Location
    Nice, France
    Posts
    2,694
    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... :/
    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

Page 2 of 6 FirstFirst 123456 LastLast
Popular pages Recent additions subscribe to a feed

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21