How would YOU explain to someone the concept of the quadtrees?
PS - I know about google too :)
Printable View
How would YOU explain to someone the concept of the quadtrees?
PS - I know about google too :)
O_o
A "Quadtree" is really just an abstract data structure.
How would you explain a tree with any other fixed number of children?
You probably wouldn't explain such a tree. You'd probably explain the reason such an abstraction exists by explaining some of the real fixtures.
So, what flavor of "Quadtree" are you discussing?
Soma
I am extremely interested in this paper and I want to implement its algorithm in page 580, but I have to study a lot first.
About the practical part of this, is there a library extremely easy to be imported which uses quadtrees? I have heard of Boost, but I have never used it.
The way I read the description of the SplitReduce function, you actually need a 2D-tree for D-dimensional data. (It seems the paper is written using two-dimensional terms, while occasionally reminding the reader that the algorithm is not limited to just two dimensions.)
Simply put, it describes how to divide a D-dimensional cube recursively, until each sub-D-cube contains a part of the original D-polytope surface that can be described using at most T half-spaces.
I personally would start with something like
While you can omit the midpoint from the data structure, it is extremely useful near the vertices of the original D-polytope. When the vertex is within the current D-cube, set the midpoints to the vertex coordinates; this should avoid exessive recursion near vertices that are part of many faces, I believe. (Consider a polyhedron approximating a wheel, with the hub being a single vertex. That vertex may be part of thousands of faces.)Code:#define D 3 /* Dimensions */
#define T 8 /* Half-spaces */
typedef struct {
double x[D]; /* Unit normal vector */
double d; /* Minimum distance to origin */
} halfspace_t;
typedef struct node node_t;
struct node {
size_t halfspaces; /* zero */
double midpoint[D];
struct node *child[1<<D];
};
typedef struct {
size_t halfspaces; /* 1..T, inclusive */
halfspace_t halfspace[T];
} leaf_t;
typedef struct {
double minimum[D]; /* Bounding box */
double maximum[D]; /* Bounding box */
node_t *tree;
} polytope_t;
Convex D-polytopes with up to T faces can be descibed using a single leaf_t, since they can be described with the same number of half-spaces. (In fact, each face plane is simply converted to a half-space, pointing towards the inside of the convex D-polytope.)
A half-space test is trivial. If n is the unit normal for the half-space, pointing towards the included half-space, d is the minimum distance from origin to the half-space, and the dividing plane itself is considered part of the half-space, then point p is in the half-space if and only if
n · p ≥ d
Descending leaves is similarly trivial. You compare each coordinate against the midpoint, and if equal or greater, set the corresponding bit of the child pointer index to one:
The only tricky thing is to write a function that calculates the number of half-space planes within a given D-cube. This is needed when generating the node tree, to find out when the current D-cube contains a simple enough part of the original D-polytope to be described using at most T half-spaces.Code:leaf_t *descend(node_t *node, const double x[D])
{
while (node && !node->halfspaces) {
size_t i = 0, d = D;
while (d-->0)
i |= (x[d] >= node->midpoint[d]) ? 1<<d : 0;
node = node->child[i];
}
return (leaf_t *)node;
}
For a convex polytope it is simple, because you can describe the [D]-polytope using half-spaces (that match the face planes). Calculate z=n·p-d for each half-space (n,d), and each corner p. If there is at least one corner where z<0, and at least one corner where z≥0, then that half-space is within the D-cube.
For other D-polytopes, each half-space is limited by a (flat) polygon, so the above test is not sufficient. One option is to construct the polygon from the intersection of the plane and the D-cube, and check whether the two polygons (which are both flat on the half-space plane) do intersect. However, I'm sure there are easier/better/more efficient ways to do that.
So, I can not use a library? I have to build my own quadtree?
By the way, I haven't yet attend the Computational Geometry class, so almost everything seem difficult from me... I will start googling...
Technically, a quadtree is just a structure that references 4 other structures of the same type.Quote:
Originally Posted by std10093
In practice, they are frequently used for partitioning a 2D area.
Imagine a canvas, draw a line down the middle from top to bottom, and another from the middle left to the right. Viola, you have 4 quadrants, in which you may rinse and repeat.
What do you want to use quadtrees for? What they look like in code is going to depend on their application immensely, so any generic library is going to just look something like:Quote:
Originally Posted by std10093
Code:class quadtree {
quadtree* quadrants [4];
void deepen () {
int i = 4;
while (i -- > 0)
quadrants [i] = new quadtree;
}
};
I have understood the quadrant concept :)
Well, as I said in link #3, I want to implement the algorithm in page 580 :)
Regarding #4, now it seems that I am starting to understand what Nomimal says. However, I am stuck on the first step of the algorithm, with the ε-approximation. Should I construct this approximation from the given polytope P? And if so, how I am supposed to do this? Can not find an algorithm for this purpose.
Moreover, I am not really sure how the input P would look like, since I have not an input sample.
Finally, it seems that I need to build my own quadtree. Should I start from a binary tree (for example this one, which is in C but imagine the same concept in C++) and expand it? Should I augment it with more functions??
//So many questions :/
I think it's just a quirk the authors had to employ to be mathematically correct. Since the original polytope is only described using half-spaces (planes matching the polytope faces, with the "inner" half-space being towards the insides of the original polytope), I guess it's technically an ε-approximation.Quote:
Originally Posted by std10093
The way I see it, the data describing the polytope in your program is the ε-approximation of the "actual" polytope.
As to the data structures in your C or C++ program, at minimum you need the half-spaces matching each face of the polytope. The half-spaces match the face planes; you need the unit normal vector of the face (pointing "out" from the polytope), and the minimum distance from the plane (if extended to infinity) to origin.
If the polytope is convex, then those face planes intersect at the polytope edges.
If the polytope is not convex, you may also need the vertices (or edges) defining each face.
Are you implementing the algorithm in 2D or 3D? In 2D, the simplest polytope is a triangle, for example (0,0),(1,0),(0,1).
In 3D, the simplest polytope is the tetrahedron, for example (0,0,0),(1,0,0),(0,1,0),(0,0,1).
The above link will show you how to get the three (2D) or four (3D) parameters you need for each face plane (and both triangles and tetrahedra are always convex, so the face planes should suffice for these), but if you want, I can give you the exact values as soon as you show your data structure you wish to hold them in.
Why not? The quad (2D) or octree (3D) you need is basically the same as your link, except that you have four or eight pointers instead of two. Usually algorithms like this one need their own implementation anyway, so you won't be able to use any existing tree functions; you'll always have to tweak them a bit. Not much.Quote:
Originally Posted by std10093
Personally, I'd write the hard stuff first, though, and write the tree-handling code when you know what kind of data structures and details you need (in the nodes). It's not like you'll need a lot of tree-related code!
For a 3D octree, consider this:
When the polytope tree is generated, you don't need any fancy manipulation functions. Aside from freeing an entire tree (as shown in the unused helper above), you don't even need to delete any nodes, ever!Code:#include <stdlib.h>
#include <math.h>
/* Point in 3D.
*/
typedef struct {
double x;
double y;
double z;
} vec3_t;
/* Plane or half-space in 3D.
* An array or half-spaces in 3D is terminated with .d == -HUGE_VAL.
*/
typedef struct {
vec3_t n;
double d;
} plane3_t;
/* Types of different nodes.
*/
typedef enum {
OUTSIDE = 0, /* Not used for convex polytopes */
INSIDE = 1, /* Not needed for convex polytopes (except boxes) */
NODE,
LEAF,
} nodetype_t;
/* Common part of all node types.
*/
struct node {
nodetype_t type;
};
/* NODE-type node.
*/
struct node_node {
nodetype_t type; /* = NODE */
vec3_t mid;
struct node *next[8];
};
/* LEAF-type node.
*/
struct node_leaf {
nodetype_t type; /* = LEAF */
plane3_t *planes; /* Array of planes */
};
/* Polytope.
*/
typedef struct {
vec3_t min;
vec3_t max;
struct node *tree;
} polytope_t;
/* Helper function: discard tree.
*/
void node_free(struct node *tree)
{
if (tree) {
if (tree->type == NODE) {
struct node_node *const n = (struct node_node *)tree;
node_free(n->next[0]);
node_free(n->next[1]);
node_free(n->next[2]);
node_free(n->next[3]);
node_free(n->next[4]);
node_free(n->next[5]);
node_free(n->next[6]);
node_free(n->next[7]);
}
free(tree);
}
}
/* Return 1 if v is inside polytope p.
*/
int inside_polytope(const vec3_t v, const polytope_t *const p)
{
const struct node *n;
/* No polytope? */
if (!p)
return 0;
/* Outside the polytope bounding box? */
if (v.x < p->min.x || v.y < p->min.y || v.z < p->min.z ||
v.x > p->max.x || v.y > p->max.y || v.z > p->max.z)
return 0;
/* Find the leaf containing v. */
n = p->tree;
while (1) {
/* NULL leaf pointer indicates OUTSIDE. */
if (!n)
return 0;
/* Entire box outside the polytope? */
if (n->type == OUTSIDE)
return 0;
/* Entire box inside the polytope? */
if (n->type == INSIDE)
return 0;
/* Leaf node? */
if (n->type == LEAF) {
const plane3_t *list = ((const struct node_leaf *)n)->planes;
/* No planes means INSIDE. (Paranoid check.) */
if (!list)
return 1;
/* List is terminated with .d == -HUGE_VAL.
* If v is outside any of the half-spaces defined in the list,
* return 0 (immediately). v is inside the polytope only
* if it is inside all the half-spaces relevant to this box.
*/
while ((*list).d > -HUGE_VAL)
if ((*list).n.x*v.x + (*list).n.y*v.y + (*list).n.z*v.z > (*list).d)
return 0;
else
list++;
/* Not outside any of the half-spaces, therefore INSIDE. */
return 1;
}
/* Non-leaf node? */
if (n->type == NODE) {
const struct node_node *const nn = (const struct node_node *)n;
n = nn->next[ 1*(v.x >= nn->mid.x)
+ 2*(v.y >= nn->mid.y)
+ 4*(v.z >= nn->mid.z) ];
continue;
}
/* This code is never reached, unless you add new node types. */
return 0;
}
}
There are two real problems you should concentrate on solving, instead:
First, how to efficiently subdivide the box. You can always just halve it, but near vertices, that may cause a VERY deep hierarchy. (My suggestion is to always subdivide the box so that there is at most one vertex; and if there is a vertex, the vertex is at a corner of the box.)
Second, how to efficiently determine which faces are relevant to the current box. (I've actually described one possible solution for convex polytopes already in this thread.)
Everything else should be very straightforward. In particular, seeing that you have a binary tree structure well in hand (I didn't actually review your code, though!), you shouldn't have any issues with the tree stuff at all.
Don't let the term "quadtree" or "octree" scare you. If you are familiar with pointers, they're very tame beasts. In this particular instance, they even behave very, very nicely -- no tree manipulations needed, aside from building a tree and tearing a tree down.
I am very confident with pointers and binary trees :D
The code is going to be in C++.
The paper assumes that the polytope is convex, so let's head that way. As a matter of fact, here is the sentence:
"Given a convex polytope P in R^d, presented as the intersection of halfspaces..."
Can one suggest how the input would look like?
I will go for a 2D approach now and then, jumping into 3D would not be that hard I guess, since at that time, I will have felt the algorithm.
Now, it's 3 o'clock at the night, so tomorrow, I should get started...right? :D
> (I didn't actually review your code, though!)
The code of the binary tree was given to us by a very very good² professor and I have used it many times, so take as granted it is ok. I am going to augment the binary tree right now :)
I'd recommend going straight for the 3D case, since there are many simplifications you can do in 2D that you cannot do in 3D, and they're easier to discover by going from the complex to simple.Quote:
Originally Posted by std10093
The obvious choice is to define the polytope using the vertices, as they are enough to define the faces for a convex polytope.Quote:
Originally Posted by std10093
You need the unit normal of each face, pointing away from the polytope, and the minimum distance from each face plane (extended to infinity) to origin.
If three points p1, p2, and p3 are non-collinear vertices of the same face in 3D, then
v = (p2 - p1) × (p3 -p1)If n points outwards from the convex polyhedron, then
n = v / sqrt(v · v)
d = n p1
p · n > dif and only if point p is outside the polyhedron (by virtue of the point p being outside the face of the convex polyhedron).
Given all half-spaces of a polyhedron, point p is outside the polyhedron if it is outside any of the half-spaces.
For example, here is a program that reads in any number of 3D vertices from a file, and outputs both the (read) vertices, and the half-spaces formed by the faces:
The logic used in the above program is very simple.Code:#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <errno.h>
#include <math.h>
typedef struct {
double x;
double y;
double z;
} point_t;
typedef struct {
double x;
double y;
double z;
double d;
} halfspace_t;
/* Read one or more 3D points from the specified stream,
* until unreadable input is encountered.
*/
size_t read_points(point_t **const listptr,
size_t *const sizeptr,
FILE *const in)
{
point_t *list, p;
size_t size, used;
/* Invalid parameters? */
if (!listptr || !sizeptr || !in) {
errno = EINVAL;
return (size_t)0;
}
/* Error in input? */
if (ferror(in)) {
errno = EIO;
return (size_t)0;
}
/* Already encountered end of stream? */
if (feof(in)) {
errno = 0;
return (size_t)0;
}
/* If list is NULL, size is always 0. */
list = *listptr;
size = (list) ? *sizeptr : 0;
used = 0;
/* Read loop. */
while (fscanf(in, " %lf %lf %lf", &p.x, &p.y, &p.z) == 3) {
/* Grow the list if necessary. */
if (used >= size) {
size = (used | 127) + 129;
list = realloc(list, size * sizeof *list);
if (!list) {
errno = ENOMEM;
return (size_t)0;
}
*listptr = list;
*sizeptr = size;
}
/* Add to list. */
list[used++] = p;
}
/* Done. */
errno = 0;
return used;
}
/* 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 - p3.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_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;
/* Assume origin is inside the convex polyhedron. */
if (curr.d < 0.0) {
curr.x = -curr.x;
curr.y = -curr.y;
curr.z = -curr.z;
curr.d = -curr.d;
}
/* 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++;
if (d < -epsilon)
n_in++;
}
/* There must be at least one point inside,
* and no points outside, for this to be valid. */
if (n_out > 0 || n_in < 1)
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;
}
int main(int argc, char *argv[])
{
point_t *vertex_list = NULL;
size_t vertex_count;
size_t vertex_maxcount = 0;
halfspace_t *face_list = NULL;
size_t face_count;
size_t face_maxcount = 0;
size_t i;
double eps;
char dummy;
if (argc < 2 || argc > 3 || !strcmp(argv[1], "-h") || !strcmp(argv[1], "--help")) {
fprintf(stderr, "\n");
fprintf(stderr, "Usage: %s [ -h | --help ]\n", argv[0]);
fprintf(stderr, " %s EPSILON [ VERTEX-FILE ]\n", argv[0]);
fprintf(stderr, "\n");
fprintf(stderr, "This program reads four or more 3D points from VERTEX-FILE.\n");
fprintf(stderr, "These points should define a convex polyhedron, with origin\n");
fprintf(stderr, "inside the polyhedron.\n");
fprintf(stderr, "\n");
fprintf(stderr, "If VERTEX-FILE is - or omitted, then the points are read\n");
fprintf(stderr, "from standard input.\n");
fprintf(stderr, "\n");
fprintf(stderr, "Epsilon defines the desired precision. Everything between\n");
fprintf(stderr, "-EPSILON and +EPSILON, inclusive, is considered zero.\n");
fprintf(stderr, "\n");
fprintf(stderr, "The program outputs the half-spaces that define the polyhedron,\n");
fprintf(stderr, "with one half-space per line. The four components on each line\n");
fprintf(stderr, "are the unit normal vector for the half-space (x y z), and\n");
fprintf(stderr, "the minimum distance (d) from the infinite plane to origin.\n");
fprintf(stderr, "If (hx,hy,hz,hd) are those four parameters, and\n");
fprintf(stderr, "\thx*x + hy*y + hz*z > hd\n");
fprintf(stderr, "then point (x,y,z) is outside the half-space.\n");
fprintf(stderr, "\n");
return 1;
}
if (sscanf(argv[1], " %lf %c", &eps, &dummy) != 1 || eps < 0.0) {
fprintf(stderr, "%s: Epsilon must be a nonnegative number.\n", argv[1]);
return 1;
}
/* Read the vertices. */
if (argc > 2 && strcmp(argv[2], "-")) {
FILE *in;
in = fopen(argv[2], "rb");
if (!in) {
fprintf(stderr, "%s: %s.\n", argv[2], strerror(errno));
return 1;
}
vertex_count = read_points(&vertex_list, &vertex_maxcount, in);
if (!vertex_count && errno) {
fprintf(stderr, "%s: %s.\n", argv[2], strerror(errno));
return 1;
}
if (ferror(in) || fclose(in)) {
fprintf(stderr, "%s: %s.\n", argv[2], strerror(EIO));
return 1;
}
if (vertex_count < 3) {
fprintf(stderr, "%s: Need at least three vertices.\n", argv[2]);
return 1;
}
} else {
vertex_count = read_points(&vertex_list, &vertex_maxcount, stdin);
if (!vertex_count && errno) {
fprintf(stderr, "standard input: %s.\n", strerror(errno));
return 1;
}
if (ferror(stdin)) {
fprintf(stderr, "standard input: %s.\n", strerror(EIO));
return 1;
}
if (vertex_count < 3) {
fprintf(stderr, "standard input: Need at least three vertices.\n");
return 1;
}
}
/* Construct faces based on the vertices. */
face_count = halfspaces(&face_list, &face_maxcount, vertex_list, vertex_count, eps);
if (!face_count && errno) {
fprintf(stderr, "%s.\n", strerror(errno));
return 1;
}
/* Print vertices. */
printf("%lu vertices: x y z\n", (unsigned long)vertex_count);
for (i = 0; i < vertex_count; i++)
printf("\t%12.6f %12.6f %12.6f\n",
vertex_list[i].x, vertex_list[i].y, vertex_list[i].z);
/* Print half-spaces. */
printf("%lu half-spaces: x y z d\n", (unsigned long)face_count);
for (i = 0; i < face_count; i++)
printf("\t%9.6f %9.6f %9.6f %12.6f\n",
face_list[i].x, face_list[i].y, face_list[i].z, face_list[i].d);
/* Discard vertex list. */
free(vertex_list); vertex_list = NULL;
vertex_count = vertex_maxcount = 0;
/* Discard face list. */
free(face_list); face_list = NULL;
face_count = face_maxcount = 0;
return 0;
}
It constructs the half-space defined by each unique triplet of points, using the formula I showed above. The triple loop over all vertices is such that 0 <= i1 < i2 < i3 < vertices. Because the program makes the added assumption that origin is within the convex polyhedron, the face normal vector must point away from origin; if not, we reverse the half-space. That way this triple loop suffices.
If a face is defined by more than three vertices, the above would generate duplicate half-spaces. However, an extra loop verifies that the candidate half-space curr is not already known. (Because there are fewer faces than vertices, it is better to do this check earlier rather than later.)
If there is at least one vertex "inside" the candidate half-space, and none "outside", then that half-space is one that defines the convex polytope. (This omits "flat" extrusions from the polytope. You could accept all half-spaces that have no vertices "outside", to include those "flat" extrusions, too. I didn't, because I don't think convex polytopes should have such features.)
The epsilon parameter to the program is important to detect collinear or nearly-collinear points. For double-precision floating point numbers (that have 53 bits of precision) and three dimensions, I suggest something like one thirty-millionth of the maximum coordinate difference (because 3*(30 million)2 is close to but under 252), unless you know your precision. If you know your coordinates are defined to three three decimal places, then something like 0.0001 is a good choice.
The above code and logic should be fairly straightforward, barring any brainfarts on my part. If you see any, or need clarification on something, just say so.
In the general case, there are many better algorithms that can be used to construct the faces based on vertices alone. The most well-known is probably Delaunay triangulation. It generates a triangle mesh of the convex hull of the specified points. Since some of the triangles may be co-planar -- for example, given a cube you'll get two co-planar triangles per face -- you need to ignore duplicate half-spaces (like I do in the above program, although again, there are more efficient methods too). Of course, you also need to make sure the face normal points away from the polytope. (For this, the best option is to find (or decide, like I did) one point that is known to be inside the convex polytope. Then, all the face normals must point away from that point.)
Hope this helps.
Thanks for all your answers Nomimal and yes, it helps². (brainfarts made me laugh :D ). I remind that we investigate only the convex polytopes.
Well, the algorithm is recursive, thus I guess in every moment of the recursion the ε-approximation will be ok, as you mentioned.Quote:
Originally Posted by Nominal Animal
What does this code do? :)Quote:
Originally Posted by Nominal Animal
I agree. By setting the z parameter to z, I jump into 2D easily.Quote:
Originally Posted by Nominal Animal
Where the methods are in the link planes you provided before.Quote:
Originally Posted by Nominal Animal
Since we calculate the n pointing outwards, we do not need to check whether n points outwards or not, right? If not, let me know.Quote:
Originally Posted by Nominal Animal
It picks the correct octant to descend into.Quote:
Originally Posted by std10093
Because the code is 3D, the tree is an octree (23=8), where each node has eight children. (If you halve a cube along each axis, you get eight subcubes.)
Instead of just halving each edge at the midpoint, I store the cut coordinate in the mid member. It'll be halfway the cube in many cases, but near vertices, it can be used to limit the depth of the tree.
The subcube to descend into is
Instead of writing out the eight (24) conditions, I used the C rules instead: (int)(boolean logic) evaluates to 1 if it is true, and 0 if false.Code:0 if v.x < nn->mid.x && v.y < nn->mid.y && v.z < nn->mid.z
1 if v.x >= nn->mid.x && v.y < nn->mid.y && v.z < nn->mid.z
2 if v.x < nn->mid.x && v.y >= nn->mid.y && v.z < nn->mid.z
3 if v.x >= nn->mid.x && v.y >= nn->mid.y && v.z < nn->mid.z
4 if v.x < nn->mid.x && v.y < nn->mid.y && v.z >= nn->mid.z
5 if v.x >= nn->mid.x && v.y < nn->mid.y && v.z >= nn->mid.z
6 if v.x < nn->mid.x && v.y >= nn->mid.y && v.z >= nn->mid.z
7 if v.x >= nn->mid.x && v.y >= nn->mid.y && v.z >= nn->mid.z
Obviously -- I mean, you said that n point outwards, so n point outwards!Quote:
Originally Posted by std10093
But yes, we only need to check n when we initially generate the half-spaces from the vertices, never afterwards.
The issue when generating the half-spaces from the vertices is that the direction of n (and n') in
n' = (p2 - p1) × (p3 - p1)does depend on the order of the points p1, p2, and p3 around the face -- it'll point outwards if the points are in counterclockwise order, and inwards if in clockwise order.
n = n / || n' ||
d = n · p1
Now, the triple loop 0 <= i1, i2, i3 < N does N3 iterations. I used 0 <= i1 < i2 < i3 < N, which only does N3/6-N2/2+N/3 iterations, or less than one sixth, but only iterates over each triplet once. Because of that, I had to verify the direction of my n vectors. (If I had done the full triple loop, I would have gotten all possible combinations; then I wouldn't have had to check the orientation. But cutting the number of iterations to under one-sixth was more sensible to me.)
As long as each of your half-spaces is such that there is at least one vertex inside, and no vertices outside, they have the correct orientation.
So, I am playing around with the code in post #12.
I input
0 0 0
1 0 0
0 1 0
0 0 1
but the halfspaces are zero, because of this piece of code, which makes the loop to continue always.
Code:/* There must be at least one point inside,
* and no points outside, for this to be valid. */
if (n_out > 0 || n_in < 1)
continue;
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).
The corresponding halfspaces() function isCode: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
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; 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;
}
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.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;
}
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.
I agree with all you said. I am playing around with the code and check this
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:-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)
Consider a cube (-1,-1,-1)-(+1,+1,+1):Quote:
Originally Posted by std10093
The six half-spaces areCode:-1 -1 -1
-1 -1 +1
-1 +1 -1
-1 +1 +1
+1 -1 -1
+1 -1 +1
+1 +1 -1
+1 +1 +1
(x,y,z) of the half-space is the unit normal vector of the face.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
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,
The half-spaces are nowCode: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 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.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)
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) = 0The right side is exactly equivalent to our use of d, i.e. we use
⇔ n·p = n·p0
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.
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 :D
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.
In this thread, I have used:Quote:
Originally Posted by std10093
- 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).
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.)Quote:
Originally Posted by std10093
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?
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. ;)
Exactly.Quote:
Originally Posted by std10093
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.
You know that's my spirit too ;)Quote:
Originally Posted by Nominal Animal
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?
Yes.Quote:
Originally Posted by std10093
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).
Yes, exactly.Quote:
Originally Posted by std10093
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
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.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 */);
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
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.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 */);
(Darn it, I don't understand why my posts get this long.)
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 :)
the bounding box is (0,0,0) - (1,1,1).Quote:
Originally Posted by std10093
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.
Right; sorry for speeding ahead.Quote:
Originally Posted by std10093
Note, you don't need the bounding box if you are generating the halfspace array based on the vertices.
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 :D // 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
, thus you check if it has to go left or right, but in the octree?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 */
}
Yup, and you end up with a tree with one node only.Quote:
Originally Posted by std10093
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
This generates the 20 half-spaces, and if you use t=4 or t=8, you should have a good test case.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
No illusion, that is completely correct.Quote:
Originally Posted by std10093
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.
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.Quote:
Originally Posted by std10093
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.)
Whenever you descend into a tree, you have a value you compare to the value on each node, to decide which way you go.Quote:
Originally Posted by std10093
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
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 isCode:a < d ->next[0] == ->left
a >= d ->next[1] == ->right
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.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]
For an octree, you have (a,b,c) and the node has (d,e,f), and you pick the node
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.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]
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.
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... :/
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):
However, it might allow a better logic of how to divide the box efficiently.Code:typedef struct {
point_t min;
point_t max;
point_t normal;
double distance;
} halfspace_t;
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:
So, now, which are the bounding boxes?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)
I phrased it badly. A halfspace is infinite; it splits the entire space into two halves, and therefore it cannot have a "bounding box".Quote:
Originally Posted by std10093
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
The faces, and the bounding box per face I was talking about, would beCode:0 0 0
1 0 0
0 1 0
0 0 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).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)
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:
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.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.
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.
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:Quote:
Originally Posted by std10093
If you compile and run it, it generates something likeCode:#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;
}
except that the pointer labels ("0x...") will vary, of course. If you save that as say example.dot, and then runCode: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=">" ]
}
the end result will look like this:Code:dot -Tpng example.dot > example.png
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
from the same data will generate an image likeCode:neato -Tpng example.dot > neato.png
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. :)
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.Quote:
Originally Posted by std10093
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.
Note that in my method, I did not construct subpolytopes. My logic was more like:Quote:
Originally Posted by std10093
- 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.
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 :D
No, Graphviz is available on all major OSes, including Windows and Mac OS X, not just Linux.Quote:
Originally Posted by std10093
Ah. I obviously didn't quite catch your idea.Quote:
Originally Posted by std10093
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.
Will do.Quote:
Originally Posted by std10093
I'd be delighted to, of course :)Quote:
Originally Posted by std10093
Let's start with a cube with vertices at ±1, ±1, ±1. The six half-spaces are
The bounding box is something like (ϵ-1,ϵ-1,ϵ-1)-(ϵ+1,ϵ+1,ϵ+1), where ϵ depends on how you test for inclusion in the bounding box.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)
In my case, I use the rule xmin ≤ x < xmax, in which case the bounding box values in C are
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).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);
(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
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.)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 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,
intoCode: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
Remember, changes in nh are not visible to the parent call, only in any recursive calls the function might make.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)
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.
We could set
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.Code:max.x = nextafter(1.0, HUGE_VAL);
// 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?
Exactly. This is very important, if min is inclusive, but max exclusive.Quote:
Originally Posted by std10093
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).
Yes, exactly.Quote:
Originally Posted by std10093
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.Quote:
Originally Posted by std10093
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,
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.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);
The child boxes will be
so that for point (x, y, z) inclusion test, the child to descend into isCode: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)
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.
Well, first I want to try this and make it work.
You have a typo here
It needs to be z, not y. :)Code:split.z = min.y + 0.5*(max.z - min.z);
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. ;)
I think I can understand this.Quote:
Originally Posted by Nominal Animal
The way I visualize it, is like this:
Attachment 12848
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.
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.Quote:
Originally Posted by Nominal Animal
However, I would expect the array to be rearranged into this:
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:-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)
Good catch. I can't edit the post anymore, so I'll have to leave it there.Quote:
Originally Posted by std10093
Yes, that's correct for the initial box.Quote:
Originally Posted by std10093
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
No. Why would you expect to see any other sides/halfspaces/planes than x=1 (+x), y=-1 (-y), and z=-1 (-z)?Quote:
Originally Posted by std10093
Just to be sure, I checked. My code gives me
whereCode: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)
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.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 ).
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.
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?
Only for some specific polytopes.Quote:
Originally Posted by std10093
Nope. Consider an octahedron, with six vertices:Quote:
Originally Posted by std10093
Yet, it has eight faces, and therefore is defined by eight halfspaces.Code:0 0 -1
+1 +1 0
-1 +1 0
-1 -1 0
+1 -1 0
0 0 +1
(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(n⌊d/2⌋)
Interesting stuff.
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).
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.)Quote:
Originally Posted by std10093
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:
Assuming I calculated right, the maximum number of faces in a D-dimensional convex polytope isCode: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 ⎠
D = 3: faces ≤ 2 × vertices - 4and so on.
D = 4: faces ≤ 1/2 × vertices2 - 3/2 × vertices
D = 5: faces ≤ vertices2 - 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
Oh I see.. Well, I do also assume that your calculations are correct, but I would like to take a shot for the D = 3 case.
The equation is this, but I do not seem to get what k represents.
EDIT:
How do we call in English the d - i and k -i that lie into a parenthesis? You can see it in the link. I know how to calculate it.
The combination one can answer. Yes, that's true.
Assume we have the combination of n and k. In Greek, we can say n per k, can we say that in English too?
fi describes the number of i-dimensional features in the convex polytope, so k-1 is the number of dimensions in the feature we are interested in; for halfspaces, k-1 = D-1.Quote:
Originally Posted by std10093
Quote:
Originally Posted by std10093
is the binomial coefficient n, k. Wikipedia says it is often read as "n choose k" in English.Code:⎛n⎞
⎝k⎠
Maybe I am tired from the code. In regards with post #41, in order to find the intersecting halfspaces with the current box of the node, isn't enough to do something like this:
?Code:if ( halfspace.x >= node.box.minBox.x && halfspace.y >= node.box.minBox.y && halfspace.z >= node.box.minBox.z
&& halfspace.x < node.box.maxBox.x && halfspace.y < node.box.maxBox.y && halfspace.z < node.box.maxBox.z )
[IRRELEVANT]
PS: Remember the architecture assignment I had to do and we figured out a way to replace division by multiplication and shifting, unroll the loop, etc.?
Well, it seemed that I had no completely get our task and I was graded badly²! However, I learnt a lot from all this experience and the knowledge is far more precious than a grade, in my case (still my GPA is good).
For the history:
- I thought that every range would get 20 numbers, but that was not a must after all, with the input tests the PhD'es had, having 22 in some and 19 in others. As a result, when the code sensed that a range did not have 20 elements, it was programmed to assume that we are out of range.
// My confusion was justified, because of a tip the professor gave inside the class.- When the data was out of range, the code was able to find this situation, but did not break the loop, which was a secondary requirement.
// I knew that, but it was not something important I guess, the first one was the real penalty.
Thanks again for the great things you thought me at that project and sorry for messing up with the grade.
[/IRRELEVANT]
No, it is not enough.Quote:
Originally Posted by std10093
The vector (halfspace.x, halfspace.y, halfspace.z) is only normal (perpendicular) to the halfspace plane; it does not by itself define the halfspace.
Consider a cube not at origin; say, at (2,3,4)-(3,4,5) instead. So, it'll have vertices at
but halfspacesCode:2 3 4
3 3 4
2 4 4
3 4 4
2 3 5
3 3 4
2 4 5
3 4 5
The initial bounding box is obviously x=[2,3], y=[3,4], z=[4,5], and none of the halfspace normal vectors are in it. But the same logic applies as to the previous cube we've discussed in this thread: all six halfspaces intersect the initial bounding box, and each of the eight subcubes will have three intersecting halfspaces.Code:-1 0 0 -2
0 -1 0 -3
0 0 -1 -4
0 0 +1 5
0 +1 0 4
+1 0 0 3
So, what you do have to do, is check whether each halfspace is inside or outside the vertex of the current box, i.e. at the eight corners of the current box:
(Of course, since floating point numbers are not exact, I personally use q<-ϵ, -ϵ<=q<=ϵ, and q>ϵ, instead.)Code:q = halfspace.x * x + halfspace.y * y + halfspace.z * z - halfspace.d
q < 0: (x, y, z) is inside the halfspace
q = 0: (x, y, z) is at the plane of the halfspace
q > 0: (x, y, z) is outside the halfspace
The WinMIPS64 assembly one, this thread? The task basically being to compute a five-bin histogram of 100 inputs, each input being between 0 and 4999 inclusive (thus each bin being 1000 integers "wide").Quote:
Originally Posted by std10093
(I think the most optimized -- quite insane -- version I sent you ran in 720 cycles per 100 inputs, using 12 bits per bin (five bins) in a 64-bit register; or quoting myself in a private message, "for any multiple of 8 Length, it runs in 18+6.75*Length cycles".)
So, your computation of the histogram was correct, but the conclusions your code made based on the histogram were incorrect?Quote:
Originally Posted by std10093
Hm. If so, I would have graded differently. (I place much more emphasis on getting the logic and computation right, than having the exact expected result.)
I do recall we discussed input range checking at some point. To me, as the lecturer/grader, detecting invalid range would have been even more important than fast execution.Quote:
Originally Posted by std10093
In any case, I do sincirely hope none of my advice caused you to veer off from the original assignment. Maybe, if we hadn't played with the different techniques and methods on how to implement a histogram calculation efficiently in 64-bit MIPS assembly, you might have concentrated more on the actual assignment, and not on the other stuff we discussed.. :(
About the halfspaces, I was also thinking that this could not be enough, because of the minimum distance to origin of the halfspace! However, it is 6 in the morning here in Greece and I was playing football and going out all the days, so I am tired I guess, but that's how things go. Thanks for pushing me a bit. ;)
[EDIT]
What are x, y and z?Code:q = halfspace.x * x + halfspace.y * y + halfspace.z * z - halfspace.d
q < 0: (x, y, z) is inside the halfspace
q = 0: (x, y, z) is at the plane of the halfspace
q > 0: (x, y, z) is outside the halfspace
We have used the information from the halfspace, what about the information of the bounding box (which has min and max coords)?
[/EDIT]
Well, the code I sent was of 820 cycles, but still the team that took 10/10 had 1100 cycles....
The edition of 720 was too insane! haha -but fun to read.
>So, your computation of the histogram was correct, but the conclusions your code made based on the histogram were incorrect?
Not really. The professor had something in course and I thought he meant that we should take as granted that every range would get exactly 20 elements.
>I place much more emphasis on getting the logic and computation right, than having the exact expected result.
Damn right! However, after discussing with others, we concluded that they really just run the code and checked the results. They did not even have a look at our documentations!
>I do recall we discussed input range checking at some point. To me, as the lecturer/grader, detecting invalid range would have been even more important than fast execution.
Me too, but I can accept this penalty, the other one was really strict.
>In any case, I do sincirely hope none of my advice caused you to veer off from the original assignment. Maybe, if we hadn't played with the different techniques and methods on how to implement a histogram calculation efficiently in 64-bit MIPS assembly, you might have concentrated more on the actual assignment, and not on the other stuff we discussed..
My friend, I want to do more than just getting a good grade. I want to learn. I want to be a scientist in our field, just like you! If the trade-off we had to dealt with was what you thought me and my grade, I was for sure going for the knowledge you gifted me.
However, I am pretty sure, that it was my fault that I lost this requirements. ;)
The coordinates you are going to test, so one corner of the bounding box in this case.Quote:
Originally Posted by std10093
Remember, if you have halfspace (halfspace.x, halfspace.y, halfspace.z, halfspace.d), and point (point.x, point.y, point.z), thenIf you check the corners of the current bounding box, and at least one corner is inside, and one corner outside, then the plane of the halfspace must intersect the current bounding box, right?Code:halfspace.x * point.x + halfspace.y * point.y + halfspace.z * point.z < halfspace.d: point is inside halfspace
halfspace.x * point.x + halfspace.y * point.y + halfspace.z * point.z = halfspace.d: point is on the plane of the halfspace
halfspace.x * point.x + halfspace.y * point.y + halfspace.z * point.z > halfspace.d: point is outside the halfspace
If none of the corners of the current bounding box is outside the halfspace, then the entire box is inside that halfspace, and that halfspace is not relevant to testing whether any point in the current box is inside the convex polytope.
If all the corners of the current bounding box are outside, for any halfspace, then the entire box is outside the convex polytope.
(This is because it is enough for a point to be outside any halfspace for it to be outside the convex polytope. Conversely, a point must be inside (or on the plane of) all halfspaces for it to be inside the convex polytope.)
Pity that thread is closed, I would have liked to post that code. (I do think I have the file still.) The technique itself is pretty much insane -- constructing a five-bin histogram in a single 64-bit register, then using several of those registers in parallel to unroll a loop. It's one of the those things that are possible, but unless it is write-only code (and the speed is absolutely required), it is not sane (maintainable, readable) -- even a sane description of how it works is probably longer than the code itself.Quote:
Originally Posted by std10093
That's a pity -- but really, that's how this world works.Quote:
Originally Posted by std10093
Usually, it is better to answer the question that the asker believes they asked, instead of answering either the question they literally asked, or the question they should have asked. (The second gives you a reputation as a difficult smart-alec, and many people cannot handle "a student" doing the third to them. Especially professors who have already fallen off the tip of their field, and need crutches to maintain the illusion in their head. Beware of anyone who thinks that they've learned enough when they get their first full professorship; they're the worst.)
In other words, you need to learn to read what the asker attempted to ask. I treat it as a puzzle; I don't have the social skills to determine it directly.
I did no such thing! :mad:Quote:
Originally Posted by std10093
I don't think knowledge can be given, or really even shared. A friend, tutor, teacher, mentor can help you discover it, but that's it.
(The best ones can help you discover stuff they don't know themselves; I've had the privilege to meet a few of those. I'm definitely not one of those, as I lack the social talent they have.)
I merely described some of the techniques possible, and gave you a bit of feedback on your implementation. In particular, I didn't show you my version until you had submitted yours (you didn't want me to, either). If I were the lecturer on that course, I'd hope the TAs were there to give you the assistance I did here. Answering theoretical questions and assessing particular details will help the students apply the course content in their assignments; heck, that's why TAs are needed!
(Of course, many professors use TAs to grade assignments, and don't care a bit whether the TAs help the students at all. Most students are too timid to complain about bad TAs to their lecturer, even if a TA spews utter incorrect bull instead of explaining the correct solution. The damage one bad apple can do in such a situation is surprisingly widespread, so stuff like that should be caught and avoided, even if handled quietly.)
In your case, I suspect the TAs did most of the grading, too.
Actually, I've had something very similar happen to me on a C course, well over a decade ago. (It was a basic sort command assignment; I discovered that I could significantly cut the wall clock time used if I read the entries into a tree, instead of using an array or a list and sorting when all data has been read.) I talked to the lecturer, specifically not to adjust the grade, but to ask if he understood why I made the choices I had.
I was much happier after the short chat, when I was sure the lecturer had understood me (and my work) correctly.
Perhaps you could do the same? It's best to mention first that you don't want them to adjust your grade, you just want to discuss your solution. (At least here, it's very common for students to try and whine a higher grade, even cry a little bit. Some are even proud of how good they are at it. Best make sure you're not one of those, in my opinion.)
Yes I remember when the point is inside a halfspace, because we had discussed it when we were up determining where a point lies inside a polytope.
Now, the question is, if a halfspace intersects a bounding box? I had in mind, that for a bounding box, it is enough to have a 3D point describing its min boundary and one for its max.
And I can't figure it out.
>I did no such thing!
>I don't think knowledge can be given, or really even shared. A friend, tutor, teacher, mentor can help you discover it, but that's it.
I phrase it badly, but that's what I meant. You show me the way in finding the path that leads to this knowledge. ;)
>In your case, I suspect the TAs did most of the grading, too.
True, but every professor has 200-300 students at my university, so I pretty much understand this.
> Most students are too timid to complain about bad TAs to their lecturer
To be fair, there are really many students that complain about TAs, because they do not actually do their homework. This is really bad and TAs are gotten accused unfairly. I saw this, while I had my scholarship in Switzerland and the university was a private one.
I did talk to him, but he was not positive.
Yes.Quote:
Originally Posted by std10093
You still have to test each of the eight corners of the box,
to know whether the halfspace plane intersects the box or not.Code:min.x min.y min.z
max.x min.y min.z
min.x max.y min.z
max.x max.y min.z
min.x min.y max.z
max.x min.y max.z
min.x max.y max.z
max.x max.y max.z
If you were only to test the two points (first and last above), you'd only know if the halfspace intersected the diagonal of the box. Yet, there are many possible planes that do intersect the box but do not intersect the diagonal.
Pity.Quote:
Originally Posted by std10093
Well, you did learn a very valuable lesson, though: always make sure you understand how your task will be evaluated, and direct your efforts accordingly.
It's a good rule to go by even if you do the evaluation yourself.
Oh you meant the corners of the bounding box! No idea why I could not get this earlier.
So, I have my function for testing if a halfspace intersects a bounding box like this:
Is it correct?Code:**
* Increase counters if needed.
* @param q - compare it with epsilon
* @param in_at - counter for points inside,
* or at the plane of the halfspace
* @param out - counter for points outside
* the halfspace
*/
void Guadtree::incCounters(const double& q, size_t& in_at, size_t& out)
{
/*
* q < 0: (x, y, z) is inside the halfspace.
* q = 0: (x, y, z) is at the plane of the halfspace.
* q > 0: (x, y, z) is outside the halfspace.
* Use epsilon, because of double variables.
*/
if(q < -epsilon || (-epsilon <= q && q <= epsilon))
in_at++;
else
out++;
}
/**
* Test if a halfspace intersects a bounding box.
* @param halfspace - the halfspace to be tested
* @param minBox - the min of the bounding box
* @param maxBox - the max of the bounding box
* @return - answer of test
*/
bool Guadtree::halfspaceIntersectsBox(const Halfspace& halfspace, const Point3d*& minBox, const Point3d*& maxBox)
{
size_t in_at = 0; // Counter for corners inside, or at the plane, of the halfspace
size_t out = 0; // Counter for corners outside of the halfspace
const Point3d halfspaceCoords = halfspace.getCoords();
const double halfspaceD = halfspace.getD();
/* Check all the corners of the bounding box, eight in number.
* The corners are:
* min.x min.y min.z
* max.x min.y min.z
* min.x max.y min.z
* max.x max.y min.z
* min.x min.y max.z
* max.x min.y max.z
* min.x max.y max.z
* max.x max.y max.z
*/
double q = halfspaceCoords * *minBox - halfspaceD;
incCounters(q, in_at, out);
q = halfspaceCoords.x() * maxBox->x() + halfspaceCoords.y() * minBox->y()
+ halfspaceCoords.z() * minBox->z() - halfspaceD;
incCounters(q, in_at, out);
/* If at least one corner is inside(or at the plane), and one corner outside,
* then the plane of the halfspace intersects the bounding box.
*/
if(in_at >= 1 && out == 1)
return true;
q = halfspaceCoords.x() * minBox->x() + halfspaceCoords.y() * maxBox->y()
+ halfspaceCoords.z() * minBox->z() - halfspaceD;
incCounters(q, in_at, out);
if(in_at >= 1 && out == 1)
return true;
q = halfspaceCoords.x() * maxBox->x() + halfspaceCoords.y() * maxBox->y()
+ halfspaceCoords.z() * minBox->z() - halfspaceD;
incCounters(q, in_at, out);
if(in_at >= 1 && out == 1)
return true;
q = halfspaceCoords.x() * minBox->x() + halfspaceCoords.y() * minBox->y()
+ halfspaceCoords.z() * maxBox->z() - halfspaceD;
incCounters(q, in_at, out);
if(in_at >= 1 && out == 1)
return true;
q = halfspaceCoords.x() * maxBox->x() + halfspaceCoords.y() * minBox->y()
+ halfspaceCoords.z() * maxBox->z() - halfspaceD;
incCounters(q, in_at, out);
if(in_at >= 1 && out == 1)
return true;
q = halfspaceCoords.x() * minBox->x() + halfspaceCoords.y() * maxBox->y()
+ halfspaceCoords.z() * maxBox->z() - halfspaceD;
incCounters(q, in_at, out);
if(in_at >= 1 && out == 1)
return true;
q = halfspaceCoords * *maxBox - halfspaceD;
incCounters(q, in_at, out);
if(in_at >= 1 && out == 1)
return true;
}
What is your opinion about this code? It is C++.
>Well, you did learn a very valuable lesson, though: always make sure you understand how your task will be evaluated, and direct your efforts accordingly.
>It's a good rule to go by even if you do the evaluation yourself.
We learn and get better constantly! :D
//At 7.30 I am leaving for vacation, for a week and I decided not to take a computer/laptop with me. I have to work at my aunt's restaurant too.
I am also very concerned on how to print the octree! What do you suggest? What should it be my approach?
// I am a travelling in two hours, so I must go to bed (of course, I will eat before going to bed :)). See you in one week! :)
A couple of notes:Quote:
Originally Posted by std10093
- It's Quad, not Guad. :)
- To handle the cases efficiently, you probably should differentiate all three cases: inside, on the plane, and outside.
You see, if one or more of the corners of the bounding box are on the plane, but none are inside or outside, then that halfspace is irrelevant for the current box.- You check the results too early.
You need to check all corners first, and then determine whether the halfspace intersects the bounding box. You don't know the result before checking all eight corners.
As a comparison, here is the relevant bit of how I'd do it:
Note the differences to your code. The loop is over each halfspace.Code:halfspace_t *const h; /* Array of halfspaces */
int nh; /* Number of halfspaces in the array */
int i = 0;
/* Loop over known halfspaces: */
while (i < nh) {
const double d[8] = { min.x*h[i].x + min.y*h[i].y + min.z*h[i].z - h[i].d,
max.x*h[i].x + min.y*h[i].y + min.z*h[i].z - h[i].d,
min.x*h[i].x + max.y*h[i].y + min.z*h[i].z - h[i].d,
max.x*h[i].x + max.y*h[i].y + min.z*h[i].z - h[i].d,
min.x*h[i].x + min.y*h[i].y + max.z*h[i].z - h[i].d,
max.x*h[i].x + min.y*h[i].y + max.z*h[i].z - h[i].d,
min.x*h[i].x + max.y*h[i].y + max.z*h[i].z - h[i].d,
max.x*h[i].x + max.y*h[i].y + max.z*h[i].z - h[i].d };
const int corners_inside = (d[0] < -eps)
+ (d[1] < -eps)
+ (d[2] < -eps)
+ (d[3] < -eps)
+ (d[4] < -eps)
+ (d[5] < -eps)
+ (d[6] < -eps)
+ (d[7] < -eps);
const int corners_outside = (d[0] > eps)
+ (d[1] > eps)
+ (d[2] > eps)
+ (d[3] > eps)
+ (d[4] > eps)
+ (d[5] > eps)
+ (d[6] > eps)
+ (d[7] > eps);
if (corners_outside >= 8) {
/*
* This sub-box is completely outside the polytope.
*
* I use NULL pointer to declare entire sub-box "outside",
* i.e. NODE_OUTSIDE_POLYTOPE.
*/
return NULL;
}
if (corners_outside == 0 || corners_inside == 0) {
/*
* Either all corners are inside the halfspace, or
* the plane of halfspace intersects a corner only.
*
* Since no part of the current sub-box is outside the polytope,
* halfspace i is irrelevant to the current sub-box.
* So, "remove" the halfspace (by moving it to the end of array).
*/
swap_halfspaces(&h[i], &h[--nh]);
continue;
}
/* Keep this halfspace. */
i++;
}
/* If there are no halfspaces left,
* this entire box is inside the polytope.
*/
if (nh < 1) {
/*
* Return a special "fully inside the polytope" node.
*/
return NODE_INSIDE_POLYTOPE;
}
The corner tests exclude the plane of the halfspace itself. (Even if we needed it, we could just count it based on there being 8 corners, and every corner falling into one of the categories; i.e. corners_at_plane = 8 - corners_inside - corners_outside.)
The determination whether a halfspace intersects the current bounding box can only be made after all eight corners have been tested. (However, if a halfspace shows that the current bounding box is completely outside the polytope, i.e. all eight corners outside, we know that the current bounding box is completely outside the polytope.)
A halfspace is only relevant to the polytope inclusion test, if it intersects the current bounding box. (If the halfspace plane intersects corners only, it is irrelevant, because it does not affect the test for any point inside the current bounding box.)
Only those halfspaces where one or more corners are inside, and one or more corners are outside, are needed for the inclusion test. My code keeps those at the start of the h array, and moves the others to the end of the list.
(We've already discussed how this allows you to use a single array of halfspaces with a recursive function.)
Just remember that in C, (true) evaluates to 1, and it is quite okay to count the number of true values using summation like I do above. In C++, converting bool to int always yields 0 or 1 for false and true, respectively, so the same works in C++, too.
If you want to make it explicit (for other programmers), you can write it as ((condition) ? 1 : 0) instead. Given sufficient optimization flags, a C compiler should generate the exact same code anyway.. but some do not.
Have a nice vacation!Quote:
Originally Posted by std10093
Tried to bypass some calculations, but I searched and finally I agree with you, so let's calculate data for all the corners.
Basically, I want a function that tells me if a halfspace intersects with the current sub-box or not, thus, I did actually use pretty much your code, well, except from the corners_outside >= 8, that I wrote it as corners_outside == 8, since we have exactly 8 corners, so the > check is redundant.
And here it is, hope this is ok.
Now I am going to write the print function and actually see, if my code works, because it executes successfully, but I am worried because, I had very few segmentation faults and relevant errors. Maybe I have some some logical errors!Code:/**
* Test if a halfspace intersects a bounding box.
* @param halfspace - the halfspace to be tested
* @param minBox - the min of the bounding box
* @param maxBox - the max of the bounding box
* @return - result of test
*/
bool Quadtree::halfspaceIntersectsBox(const Halfspace& halfspace, const Point3d*& minBox, const Point3d*& maxBox)
{
const Point3d h = halfspace.getCoords(); // Coordinates of halfspace
const double hD = halfspace.getD(); // Distance of halfspace
/* Check all the corners of the bounding box, eight in number.
* The corners are:
* min.x min.y min.z
* max.x min.y min.z
* min.x max.y min.z
* max.x max.y min.z
* min.x min.y max.z
* max.x min.y max.z
* min.x max.y max.z
* max.x max.y max.z
*/
const double d[8] = { h * *minBox - hD,
h.x() * maxBox->x() + h.y() * minBox->y()
+ h.z() * minBox->z() - hD,
h.x() * minBox->x() + h.y() * maxBox->y()
+ h.z() * minBox->z() - hD,
h.x() * maxBox->x() + h.y() * maxBox->y()
+ h.z() * minBox->z() - hD,
h.x() * minBox->x() + h.y() * minBox->y()
+ h.z() * maxBox->z() - hD,
h.x() * maxBox->x() + h.y() * minBox->y()
+ h.z() * maxBox->z() - hD,
h.x() * minBox->x() + h.y() * maxBox->y()
+ h.z() * maxBox->z() - hD,
h * *maxBox - hD };
const int corners_inside = (d[0] < -epsilon)
+ (d[1] < -epsilon)
+ (d[2] < -epsilon)
+ (d[3] < -epsilon)
+ (d[4] < -epsilon)
+ (d[5] < -epsilon)
+ (d[6] < -epsilon)
+ (d[7] < -epsilon);
const int corners_outside = (d[0] > epsilon)
+ (d[1] > epsilon)
+ (d[2] > epsilon)
+ (d[3] > epsilon)
+ (d[4] > epsilon)
+ (d[5] > epsilon)
+ (d[6] > epsilon)
+ (d[7] > epsilon);
/*
* A halfspace is only relevant to the polytope inclusion test,
* if it intersects the current bounding box. (If the halfspace
* plane intersects corners only, it is irrelevant, because it
* does not affect the test for any point inside the current
* bounding box). Only those halfspaces where one or more corners
* are inside, and one or more corners are outside, are needed
* for the inclusion test.
*/
if(corners_inside == 0 || corners_outside == 0 || corners_outside == 8)
return false;
return true;
}
EDIT:
Here is the print function and now, let's start testing and... debugging of course! :D
Thanks for the wish for the vacations. :)Code:/**
* Print the tree.
*/
void Quadtree::print()
{
print(root);
}
/**
* Print the current node.
* @param node - the node to be printed.
*/
void Quadtree::print(node*& node)
{
if(node)
{
print(node->child[0]);
std::cout << node->split << std::endl;
print(node->child[1]);
std::cout << node->split << std::endl;
print(node->child[2]);
std::cout << node->split << std::endl;
print(node->child[3]);
std::cout << node->split << std::endl;
print(node->child[4]);
std::cout << node->split << std::endl;
print(node->child[5]);
std::cout << node->split << std::endl;
print(node->child[6]);
std::cout << node->split << std::endl;
print(node->child[7]);
std::cout << node->split << std::endl;
}
}
So, lets see:Quote:
Originally Posted by std10093
- corners_inside == 0
None of the corners are inside the box.
Most of the box is therefore outside the halfspace; one to four corners, up to four edges, and possibly one face of the box may be on the plane of the halfspace, though.- corners_outside == 0
All corners are either inside or on the plane of the halfspace.
There are no points inside this box that are outside this halfspace.- corners_outside == 8
All eight corners are outside the halfspace, so the entire box is outside the halfspace.
Therefore, this box is completely outside the polytope.
The latter two are clearly okay, I'm only worried about the first one. It really depends on what do you use this test for?
(corners_inside == 0) can be true if one or more corners are on the plane of the halfspace. This means that even if corners_inside==0, the halfspace may be coplanar with one side of the box. Points exactly on the halfspace plane are usually considered inside.
Now that I think about it, I think the line 44 in the code I listed in post #56 is incorrect; it should only test for (corners_outside==0).
Consider the case where the halfspace is coplanar to a side of the box. Those four corners are on the plane (neither inside nor outside), and the other four are outside. So, corners_inside==0 and corners_outside==4.
I think that kind of a halfspace should be considered to intersect the box.
The same logic applies even if the halfspace plane intersects only two points (coplanar to an edge of the box) or one point (intersects at that corner only).
I'm not absolutely certain yet, though. It needs testing, something I think I neglected to do for the code I showed in #56. In particular, I think that the points that happen to sit exactly on the plane should be considered "inside" -- even if every other point within the box is considered outside. Hmm.. needs more thinking.
Yes, I do agree that this box is relevant, since it has its 4 corners on the plane. I am going to think of it more however.
>It really depends on what do you use this test for?
When we are building the tree, we want to split a node, only if the number of the intersecting halfspaces with the box of this node is larger than t (which is the desired query time, which we get from input).
So, for every halfspace of the polytope, when we are in the root, we check if the current halfspace intersects the current box and the answer to this question comes from the Quadtree::halfspaceIntersectsBox().
By the way, for the print of the tree, I modified the function into this:
However, I am not sure yet if this is correct. I am going to test it after my match tonight.Code:/**
* Print the tree.
*/
void Quadtree::print()
{
std::cout << "Octree" << std::endl;
print(root);
}
/**
* Print the current node.
* @param node - the node to be printed.
*/
void Quadtree::print(node*& node)
{
if(node)
{
print(node->child[0]);
std::cout << node->split << std::endl;
print(node->child[1]);
print(node->child[2]);
print(node->child[3]);
print(node->child[4]);
print(node->child[5]);
print(node->child[6]);
print(node->child[7]);
}
}
Yep. These corner cases (no pun intended) are tricky; they need hard thought to get right.Quote:
Originally Posted by std10093
It's not like it's the first time I've gotten the corner cases wrong myself, even in this thread ;)
For more or less the same purpose, then.Quote:
Originally Posted by std10093
The reason I asked is that I do my test when deciding which halfspaces are relevant to this box, and which are not. The decision whether to split the node or not is only done later, after I've checked all the halfspaces. I only need the number of relevant halfspaces for that decision.
The difference is very subtle. If a halfspace intersects the box so that all points in the box are still inside the box, the halfspace is not interesting. If a halfspace intersects the box so that the rest of the box is outside the halfspace, the halfspace is interesting. So, it's not a plain intersects-or-not test.
Naming and commenting your functions in a way that makes it clear to the programmer what it does, is very important. I should know, because I just made that error in post #56 myself!
Aaaand since I cannot edit post #56 anymore, here's the fixed version of that code snippet:
If anyone else is wondering, my habit of using ">= 8" instead of "== 8" is perhaps silly, but it brings me comfort. I've seen too many bugs caused by using "==" when ">=" is the correct choice (for some rare or exceptional input), that going well and truly the other way gives me the warm fuzzy feeling. :biggrin:Code:halfspace_t *const h; /* Array of halfspaces */
int nh; /* Number of halfspaces in the array */
int i = 0;
/* Loop over known halfspaces: */
while (i < nh) {
const double d[8] = { min.x*h[i].x + min.y*h[i].y + min.z*h[i].z - h[i].d,
max.x*h[i].x + min.y*h[i].y + min.z*h[i].z - h[i].d,
min.x*h[i].x + max.y*h[i].y + min.z*h[i].z - h[i].d,
max.x*h[i].x + max.y*h[i].y + min.z*h[i].z - h[i].d,
min.x*h[i].x + min.y*h[i].y + max.z*h[i].z - h[i].d,
max.x*h[i].x + min.y*h[i].y + max.z*h[i].z - h[i].d,
min.x*h[i].x + max.y*h[i].y + max.z*h[i].z - h[i].d,
max.x*h[i].x + max.y*h[i].y + max.z*h[i].z - h[i].d };
const int corners_inside = (d[0] < -eps)
+ (d[1] < -eps)
+ (d[2] < -eps)
+ (d[3] < -eps)
+ (d[4] < -eps)
+ (d[5] < -eps)
+ (d[6] < -eps)
+ (d[7] < -eps);
const int corners_outside = (d[0] > eps)
+ (d[1] > eps)
+ (d[2] > eps)
+ (d[3] > eps)
+ (d[4] > eps)
+ (d[5] > eps)
+ (d[6] > eps)
+ (d[7] > eps);
if (corners_outside >= 8) {
/*
* This sub-box is completely outside the polytope.
*
* I use NULL pointer to declare entire sub-box "outside",
* i.e. NODE_OUTSIDE_POLYTOPE.
*/
return NULL;
}
if (corners_outside == 0) {
/*
* None of the corners are outside the halfspace, so
* all points inside the box are inside this halfspace.
* Therefore, this halfspace is irrelevant to this box.
* So, "remove" the halfspace (by moving it to the end of array).
*/
swap_halfspaces(&h[i], &h[--nh]);
continue;
}
/* Keep this halfspace. */
i++;
}
/* If there are no halfspaces left,
* this entire box is inside the polytope.
*/
if (nh < 1) {
/*
* Return a special "fully inside the polytope" node.
*/
return NODE_INSIDE_POLYTOPE;
}
In other words, it's just a mannerism I have. If you have a sharp eye, you'll notice many others in my code snippets.. but they all have a story behind them, and I like 'em.
So, we use for the same exact reason. I just have a function inside the loop you have.
>If a halfspace intersects the box so that all points in the box are still inside the box, the halfspace is not interesting.
I suspect, you did not want to write the italics and you have made a typo, because I do not get it.
Obviously I thought that - as you said before - A halfspace is only relevant to the polytope inclusion test, if it intersects the current bounding box.
So, I am going to rename the function to ... hmm..I can't think of anything short and illustrative now..
D'oh. Make the latter "box" "halfspace": "If a halfspace intersects the box so that all points in the box are still inside the halfspace, the halfspace is not interesting."Quote:
Originally Posted by std10093
Yes, but not all halfspaces that intersect the current bounding box are relevant to the polytope inclusion test.Quote:
Originally Posted by std10093
If the halfspace plane is such that it coincides with a corner, edge, or face of the box, but otherwise the box is inside the halfspace, then that halfspace is irrelevant to the current box: There are no points that test "outside that halfspace" within the current box.
This is exactly the case (corners_outside == 0).
(It takes quite a bit of thought to grasp correctly. Like I said, I messed this detail up earlier. It is possible I still haven't thought it out right, so I recommend you don't take my word for it.)
Some variation of "halfspace relevant" to box? :)Quote:
Originally Posted by std10093
By the way, I personally firmly believe that stuff like that -- that the code uses names and comments that drive the reader and future programmers towards the correct concepts, and perhaps even their intuition towards the correct assumptions -- is an important part of that which makes great code great.
The more I work with code, the more weight I place on the long-term viability and maintainability of it. It is certainly the direction I wish to develop my own skills.
I think I have the code working..
However, I have as test case the unit cube. In your code, does it produce a swallow or normal-depth tree? Because in mine, it goes on and on! And I think that the reason is my logic and not some logical bug.
[EDIT]
Same for the tetrahedron.
[/EDIT]
You have an array of the halfspaces and you re-arrange them, thus you only use one array and the cost of re-arrangement.
I have my node as this right now:
Where the vector holds the indexes of the array(the same array as yours) which allows for getting away form the re-arrangement cost, but requires some space. And if you think that our tree might get big, then we will have some space multiplied by the number of nodes!Code:struct Node {
Point3d split; // Split point
BoundingBox* box; // Pointer to bounding box. Free it after construction
// -1 means that all halfspaces intersect the box (for example, in the root)
std::vector<int> relevantHalfspaces; // Indexes of halfspaces that intersect the current box. The indexes are to the halfspaces vector of the
// polytope
Node* child[1<<Dimension]; // Pointers to children
};
Now, I am thinking of it, I have to destroy the vectors, in all inner leaves.
What do you say about my approach? Are you happy from yours? :)
What are you using as the maximum number of halfspaces per leaf node? I only test with 3 or larger, as smaller limits generate unnaturally deep trees; the intersection of two or three halfspaces are quantized to least-significant-unit sized bounding boxes... I definitely consider 3 the smallest sane limit to test.Quote:
Originally Posted by std10093
For a cube, I get one inner node with eight leaf nodes containing three halfspaces per leaf. (This is with perturbing the vertices very slightly, as otherwise my algorithm detects that three of the sides are coplanar with the box, and keeps only the three halfspaces, all in the root node.)
For the (0,0,0),(1,0,0),(0,1,0),(0,0,1) tetrahedron, I get two inner nodes; one with split point (0.5000000000000002,0.5000000000000002,0.500000000 0000002) with four leaf nodes completely outside the polytope, three leaf nodes with three halfspaces, and one inner node child. The inner node child has split point (0.2500000000000001,0.2500000000000001,0.250000000 0000001), with four leaf nodes with three halfspaces, three leaf nodes with two, and one with one halfspace.
The reason for the split point coordinates not being exactly half and a quarter is the nextafter() logic we discussed above: in my implementation, the upper limit is exclusive, i.e. nextafter(maxCoord, HUGE_VAL); due to rounding, the mean of the bounds is just above half.
Why would you store the bounding box in the node structure? It is only needed during construction! For the bounding box, I actually use six variables in the recursive function call, as it is easier to construct the bounding box using separate variables for the next recursive call.Quote:
Originally Posted by std10093
Hmm. You need either the halfspaces, or the child nodes, never both.Quote:
Originally Posted by std10093
To be exact, there are four different types of nodes you need to handle:
- Inner node -- has split point coordinates, and 2dimensions child nodes
- Outside leaf node -- leaf node where all points in the bounding box are outside the polytope (I use NULL pointer for this)
- Inside leaf node -- leaf node where all points in the bounding box are inside the polytope (I use an easily-detected normal leaf node for this)
- Normal leaf node -- contains or refers to halfspaces which tell which points in the bounding box are inside the polytope, and which are outside
The root node, or the structure containing the root node to the tree, must also have the bounding box of the polytope itself, so that you don't need to test points you know for sure are outside the polytope.
Aside from the superfluous bounding box in the node structure, I worry that having indexes to the halfspaces (4 or 8 bytes per index) rather than the halfspace itself (16 or 32 bytes per halfspace) is going to be a slowdown.Quote:
Originally Posted by std10093
The reason is cache behaviour. When you test for point inclusion in the polytope, you first need to descend through the tree. This is unavoidable, and a cost we just have to accept.
However, in my case, with the halfspaces contained within the leaf, that data will be cached, and immediately available to the CPU, to determine whether the point is inside the polytope or not.
In your case, the halfspaces are rarely stored near each other, so the CPU has to cache parts of the halfspace array too.
Remember: this implementation makes sense only if lots of halfspaces are needed to define the convex polytope, so we cannot assume all of the halfspace array will be cached from point inclusion test to another.
Obviously, only comparing the real-world results from both approaches would tell for sure.. but this is something I did consider when designing my node structures. So, I'm confident (that the indirect method is measurably slower for real-world use cases), but I don't know it for a fact.
As to my own algorithm implementation, I'm quite happy with it. I'm certain there is room for improvement, but it's efficient, clean, and easy-to-read as is. Before posting it, I would do a final pass of checking the comments and testing it with some random polytopes, but that's it.
Improving the split point logic from just halving would be interesting to do research on, but aside from a pure research topic, is not worth the effort.
(I suspect that first finding a case where the tree structure generated by the halving split point logic is clearly poor, would also provide enough information to devise a way to easily avoid these poor cases. I believe they occur only in cases where a small set of vertices defines many faces/halfspaces. For regular polytopes, the halving logic works really well.)
Methods of extending this to nonconvex polytopes are also very interesting, and perhaps more relevant to real-world use cases.
I see your points.
Except this one.
For this I am getting this:Quote:
For the (0,0,0),(1,0,0),(0,1,0),(0,0,1) tetrahedron, I get two inner nodes; one with split point (0.5000000000000002,0.5000000000000002,0.500000000 0000002) with four leaf nodes completely outside the polytope, three leaf nodes with three halfspaces, and one inner node child. The inner node child has split point (0.2500000000000001,0.2500000000000001,0.250000000 0000001), with four leaf nodes with three halfspaces, three leaf nodes with two, and one with one halfspace.
V = 4, epsilon = 1e-05, t = 3 <-- The goal is at most 3 halfspaces per leaf
Vertices
0,0,0
1,0,0
0,1,0
0,0,1
Halfspaces
-0 -0 -1 -0
0 -1 0 0
-1 -0 -0 -0
0.577350269189626 0.577350269189626 0.577350269189626 0.577350269189626
and then the octree, creates another level with eight children, as expected and the 6th child (i.e. the pre-last), creates more levels and goes one forever.
The 6th child has this data:
0.250000000050000,0.000000000000000,0.000000000000 000 //split point
Bounding box is: ( 0.000000000000000, 0.000000000000000, 0.000000000000000 ) - ( 0.500000000100000, 1.000000000100000, 1.000000000100000 )
and it intersects all halfspaces.
And then the 6th child of the upper 6th child has this:
0.125000000075000,0.000000000000000,0.000000000000 000
Bounding box is: ( 0.000000000000000, 0.000000000000000, 0.000000000000000 ) - ( 0.250000000150000, 1.000000000200000, 1.000000000200000 )
and it intersects all halfspaces too!
And the icosahedron and the unit cube had the same problem. In the 6th child, they would create children again and again, forever!
Vertices and halfspaces are as expected. I can even see from the split coordinates that even the initial bounding box should be correct (or at least the same as I have).
Both that split point, and that bounding box, strike me as very odd.Quote:
Originally Posted by std10093
If your original bounding box is [0,1) for all coordinates, then that bounding box is obviously wrong, right? The first-level children should have either [0,0.5) or [0.5,1) for each coordinate for the bounding box, with split point coordinates being either 0.25 or 0.75, never zero.
If you print out the first level child nodes -- their split point coordinates and bounding box coordinates, what do you get? Draw a 3D image (using e.g. the isometric perspective I mentioned above), it helps me a lot to visualize the situation correctly.
That's what I spotted too and that's why I asked you.
Here is where I am confused: When computing the bounding box of the 7th (for example) child, you have this formula: (split.x, split.y, split.z) - (max.x, max.y, max.z).
I guess, split point is from the 7th child and max from the parent, correct?
On the same child, we have to calculate the split point by this formula: split.x = min.x + 0.5*(max.x - min.x); (same for y and z).
The min and max come from the parent, correct?
I am testing, but something goes wrong..
The split point is calculated by this formula:
and the max box by these formulas: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);
*/
and there seem to be no typos..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, max.y, split.z)
4: (min.x, min.y, split.z) - (split.x, split.y, max.z)
5: (split.x, min.y, split.z) - (max.x, split.y, max.z)
6: (min.x, split.y, split.z) - (split.x, max.y, max.z)
7: (split.x, split.y, split.z) - (max.x, max.y, max.z)
Incorrect. Both are from parent.Quote:
Originally Posted by std10093
Think about it for a minute. The parent is a box, with split point somewhere inside it; the common corner for all the child boxes. Where are the edges for the eight child boxes?
Well, I have the code with "both are from parent", because of the thought you just said, but it will just not work. (giving the [0,05) and [0.5,1) you mentioned before.Quote:
Originally Posted by Nominal Animal
Consider the following. I'll use the notation 1+ to mean the smallest number greater than 1.
Let's assume current bounding box is x=[0,1+), y=[0,1+), and z=[0,1+). We discover that there are more relevant halfspaces that can be referred to in a single node, so the current node needs to be split.
We split the current bounding box at the mean of bounding box coordinates (halving the bounding box along each axis), so the current split point will be (0.5+, 0.5+, 0.5+).
So, we havemin = ( 0, 0, 0 )for the current box, right?
max = ( 1+, 1+, 1+ )
split = ( 0.5+, 0.5+, 0.5+ )
The bounding boxes for each child will be:
just like in your post #67; I only added some spaces to make it easier to read.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, max.y, split.z)
4: ( min.x, min.y, split.z) - (split.x, split.y, max.z)
5: (split.x, min.y, split.z) - ( max.x, split.y, max.z)
6: ( min.x, split.y, split.z) - (split.x, max.y, max.z)
7: (split.x, split.y, split.z) - ( max.x, max.y, max.z)
I recommend you modify your recursive function to just output the bounding box it receives, and not recurse any further. The above do produce
I definitely do not see your (0.0,0.0,0.0)-(0.5+,1.0+,1.0+) in the above list. (Which would be very surprising anyway, as it would fill the upper half of the box.)Code:0: ( 0.0, 0.0, 0.0 ) - ( 0.5+, 0.5+, 0.5+ )
1: ( 0.5+, 0.0, 0.0 ) - ( 1.0+, 0.5+, 0.5+ )
2: ( 0.0, 0.5+, 0.0 ) - ( 0.5+, 1.0+, 0.5+ )
3: ( 0.5+, 0.5+, 0.0 ) - ( 1.0+, 1.0+, 0.5+ )
4: ( 0.0, 0.0, 0.5+) - ( 0.5+, 0.5+, 1.0+ )
5: ( 0.5+, 0.0, 0.5+) - ( 1.0+, 0.5+, 1.0+ )
6: ( 0.0, 0.5+, 0.5+) - ( 0.5+, 1.0+, 1.0+ )
7: ( 0.5+, 0.5+, 0.5+) - ( 1.0+, 1.0+, 1.0+ )
The bounding box for the child 6 will be (0.0,0.5+,0.5+)-(0.5+,1.0+,1.0+), and if it is subdivided further (it shouldn't), it's split point will be (0.25+,0.75+,0.75+).
Don't get frustrated. I know this works.
I was modifying things yesterday and now I have the same boxes as you, but on the first child, the tree starts to get really deep. Is that normal? I think not! I have to investigate! I left the tree grow only one level and I saw that 1st, 2nd, 3rd and 4th child intersect all the halfspaces, thus need to be split.
Here is my octree, grown up until one level only.
If I have understood well, the split point is being computed from the current box and the current box from the parent's data.Code:Octree
0.250002500050000,0.250002500050000,0.250002500050000
Bounding box is: ( 0.000000000000000, 0.000000000000000, 0.000000000000000 ) - ( 0.500005000100000, 0.500005000100000, 0.500005000100000 )
/* That's the root */
0.500005000000000,0.500005000000000,0.500005000000000
Bounding box is: ( 0.000000000000000, 0.000000000000000, 0.000000000000000 ) - ( 1.000010000000000, 1.000010000000000, 1.000010000000000 )
0.750007500050000,0.250002500050000,0.250002500050000
Bounding box is: ( 0.500005000000000, 0.000000000000000, 0.000000000000000 ) - ( 1.000010000100000, 0.500005000100000, 0.500005000100000 )
0.250002500050000,0.750007500050000,0.250002500050000
Bounding box is: ( 0.000000000000000, 0.500005000000000, 0.000000000000000 ) - ( 0.500005000100000, 1.000010000100000, 0.500005000100000 )
0.750007500050000,0.750007500050000,0.250002500050000
Bounding box is: ( 0.500005000000000, 0.500005000000000, 0.000000000000000 ) - ( 1.000010000100000, 1.000010000100000, 0.500005000100000 )
0.250002500050000,0.250002500050000,0.750007500050000
Bounding box is: ( 0.000000000000000, 0.000000000000000, 0.500005000000000 ) - ( 0.500005000100000, 0.500005000100000, 1.000010000100000 )
0.750007500050000,0.250002500050000,0.750007500050000
Bounding box is: ( 0.500005000000000, 0.000000000000000, 0.500005000000000 ) - ( 1.000010000100000, 0.500005000100000, 1.000010000100000 )
0.250002500050000,0.750007500050000,0.750007500050000
Bounding box is: ( 0.000000000000000, 0.500005000000000, 0.500005000000000 ) - ( 0.500005000100000, 1.000010000100000, 1.000010000100000 )
0.750007500050000,0.750007500050000,0.750007500050000
Bounding box is: ( 0.500005000000000, 0.500005000000000, 0.500005000000000 ) - ( 1.000010000100000, 1.000010000100000, 1.000010000100000 )
EDIT: Wait, I found a bug.. No, nothing changed, have to look harder!
In the function that tests if a halfspace is relevant to box, if I modify it to this:
the tree becomes a tree of one level only, but this is wrong, as you have discussed...Code:if(cornersInside == 0 || cornersOutside == 8 || cornersOutside == 0)
then halfspace is relevant
If you mean that the root of your octree has bounding box
( 0.0000000000, 0.0000000000, 0.0000000000 ) - ( 1.0000100001, 1.0000100001, 1.0000100001 )and the bounding boxes for its childs are
( 0.0000000000, 0.0000000000, 0.0000000000 ) - ( 0.5000050001, 0.5000050001, 0.5000050001 )with the split point in the center of each bounding box, it looks right to me.
( 0.5000050000, 0.0000000000, 0.0000000000 ) - ( 1.0000100001, 0.5000050001, 0.5000050001 )
( 0.0000000000, 0.5000050000, 0.0000000000 ) - ( 0.5000050001, 1.0000100001, 0.5000050001 )
( 0.5000050000, 0.5000050000, 0.0000000000 ) - ( 1.0000100001, 1.0000100001, 0.5000050001 )
( 0.0000000000, 0.0000000000, 0.5000050000 ) - ( 0.5000050001, 0.5000050001, 1.0000100001 )
( 0.5000050000, 0.0000000000, 0.5000050000 ) - ( 1.0000100001, 0.5000050001, 1.0000100001 )
( 0.0000000000, 0.5000050000, 0.5000050000 ) - ( 0.5000050001, 1.0000100001, 1.0000100001 )
( 0.5000050000, 0.5000050000, 0.5000050000 ) - ( 1.0000100001, 1.0000100001, 1.0000100001 )
Hmm.. You do not enlarge the bounding box in the recursive call, do you? My argument about the necessity of enlarging the bounding box as the upper limit is exclusive and not inclusive applies only to the initial bounding box. In other words, if
xMIN = MIN( x0, x1, ..., xn-1 )with xi, yi, zi being the exact coordinates for vertex i, then the initial bounding box is
yMIN = MIN( y0, y1, ..., yn-1 )
zMIN = MIN( z0, z1, ..., zn-1 )
xMAX = MAX( x0, x1, ..., xn-1 )
yMAX = MAX( y0, y1, ..., yn-1 )
zMAX = MAX( z0, z1, ..., zn-1 )
In the recursive calls, the bounding box coordinates are treated as exact, no adjustment with nextafter(). The split point is calculated as the center (mean) of the bounding box,Code:boundingbox.min.x = xMIN;
boundingbox.min.y = yMIN;
boundingbox.min.z = zMIN;
boundingbox.max.x = nextafter(xMAX, HUGE_VAL);
boundingbox.max.y = nextafter(yMAX, HUGE_VAL);
boundingbox.max.z = nextafter(zMAX, HUGE_VAL);
which is numerically more stable than just halving their sum.Code:split.x = boundingbox.min.x + 0.5 * (boundingbox.max.x - boundingbox.min.x);
split.y = boundingbox.min.y + 0.5 * (boundingbox.max.y - boundingbox.min.y);
split.z = boundingbox.min.z + 0.5 * (boundingbox.max.z - boundingbox.min.z);
Because our lower limit of bounding box is inclusive, but the upper limit exclusive, the coordinate range min,split does not include the split point, but split,max does.
That also means that when descending the tree (in the polytope point inclusion function), we use
to determine which child (0..7) to descend into.Code:child = 1 * (testpoint.x >= split.x)
+ 2 * (testpoint.y >= split.y)
+ 4 * (testpoint.z >= split.z);
Every unique point describable using floating-point coordinates is inside exactly one of the child nodes. The split point of the parent, will be inside the last child.
If we did not do the initial bounding box upper limit adjustment, we would exclude those points of the polytope, that had the exact maximum coordinates (maximum coordinates from any vertex) -- again, because the upper limit is exclusive.
While that slice is infinitesimally thin to mathematicians (0.999... being equal to 1), it is measurable for us programmers; and that is exactly what nextafter() was designed to help with. I've seen how ignoring this kind of seemingly irrelevant details cause nastily visible seams and discontinuities in raytracing for example. Me being me, I find it important to do these details right, if you do it at all.
You make your own choices. As I've demonstrated, I am often wrong anyway.
Exactly.Quote:
Originally Posted by std10093
Ah-ha!Quote:
Originally Posted by std10093
What does your code do when there are no relevant halfspaces? I think you end up recursing forever in that case!
For the above tests, I have either expanded the polytope bounding box slightly, or perturbed the vertex coordinates slightly, so that my code does not find that three of the halfspaces are coplanar with the polytope bounding box -- for both the (0,0,0)-(1,0,0) cube, and the (0,0,0),(1,0,0),(0,1,0),(0,0,1) tetrahedron.
You see, in both cases the X=0, Y=0, and Z=0 faces are coplanar to the halfspace planes, and can therefore be ignored.
I am assuming that your code dismisses those halfspaces too, at least in the recursive case. You end up calling your recursive function for at least one of the sub-boxes, with no relevant halfspaces (because the logic has culled the halfspaces coplanar to the bounding box out).
If you change the test logic, you retain those otherwise culled halfspaces, and therefore avoid the bug.
That is absolutely correct, if the cube or tetrahedron coordinates are exact integers.Quote:
Originally Posted by std10093
The algorithm notes that three of the cube or tetrahedron faces/halfspace planes are coplanar with the bounding box (lower limits), and therefore irrelevant. This yields only three (cube) or one (tetrahedron) relevant halfspaces for the initial bounding box (entire polytope), which is obviously storeable in one node; thus yielding a "tree" with exactly one node.
This is obviously desirable for real life code, as it speeds polytope inclusion tests up, and produces shallower/faster trees.
To avoid this effect, I either enlarge the bounding box, in which case all halfspaces are inside the initial bounding box, not coplanar with any part of it, or perturb the vertex coordinates by a fraction (typically by a few hundredths).
When I test my algorithm, I test all three cases, because all three are possible in real-life data. (Well, except enlarging the initial bounding box, as I calculate that directly from the vertex coordinates.)
>What does your code do when there are no relevant halfspaces? I think you end up recursing forever in that case!
I am afraid that is not the case... I am going to continue searching and if I found nothing I am going to compare it with your implementation and hope that everything will work ok. :)
The tree is becoming huge and I can not find the bug yet, so I just left the tree grow one level, thus 9 nodes in total, 1 root and 8 children and tested the inclusion test.
I have
Test point = 1.100000000000000,0.500000000000000,0.700000000000 000
and I am getting the answer yes. I descend to the 5th child, were the first three halfspaces are relevant
Shouldn't I receive no? Since the x of the inclusion test is out of the unit cube. I have the polytope of the tetrahedron.Code:-0.000000000000000,-0.000000000000000,-1.000000000000000 <- that's the coordinates
and -0.000000000000000 <- that's d
0.000000000000000,-1.000000000000000,0.000000000000000
and 0.000000000000000
-1.000000000000000,-0.000000000000000,-0.000000000000000
and -0.000000000000000
Your test point seems to be outside the polytope bounding box, so it should not even descend into the tree.Quote:
Originally Posted by std10093
Should I post my code yet?
I've been playing with implementing a better vertices-to-halfspaces algorithm than the naive O(n4) one currently used; especially one that calculates the convex hull of the given points (and handles degenerate cases well). I can distract myself with that, if you don't want to see my code (yet).
I think it's time to get a look at your code.
I had implemented GRAHAM-SCAN algorithm for finding convex hulls, that runs in O(nlogn).
Recall that the index of the child for descending into the tree is:Quote:
Originally Posted by Nominal Animal
so if we compute the index, it is correct to be five, since the Test point is 1.1, 0.5, 0.7Code:child = 1 * (testpoint.x >= split.x)
+ 2 * (testpoint.y >= split.y)
+ 4 * (testpoint.z >= split.z);
But, we should probably first check if the test point is inside the bounding box.
This seems to do the trick:
Code:/**
* Checks if a point lies inside the bounding box of the polytope.
* @param box - the bounding box of the polytope (of the root)
* @param p - test point
* @return - result
*/
bool Quadtree::insidePolytopeBox(const BoundingBox* box, const Point3d& p)
{
for(int i = 0 ; i < D ; ++i) // Check all coordinates
if(box->getMinBox()[i] > p[i] || box->getMaxBox()[i] < p[i])
return false;;
return true;
}
No; you are forgetting that the halfspaces retained for inclusion tests in each node (or its child nodes) are only relevant to points within the bounding box of that node.Quote:
Originally Posted by std10093
Remember, we specifically store only the information relevant to that bounding box. So, descending into the tree with a point outside the original bounding box yields undefined results.
In practice, it should not matter, as every direction from inside the halfspace should be bounded by a halfspace.
Our implementations don't do that, though, because of the way we define "interesting halfspace": we ignore halfspaces that are perpendicular to the X, Y, or Z axes, on the lower-coordinate side of the polytope (the inside==0, outside==0 case; i.e. one or more corners of the box are on the halfspace plane, but no corner is clearly inside or outside the halfspace). This does not mean there should be any possibility of infinite recursion or stuff like that; it just means the tree constructed is only valid for points inside the original bounding box.
Take the (0,0,0),(1,0,0),(0,1,0),(0,0,1) tetrahedron as an example. If you test point (-1,-1,-1), ignoring the original bounding box, you get result inside. That is because the three axis-perpendicular faces are coplanar with the bounding box, and found "uninteresting" -- as they don't do any further limiting to point inclusion compared to the bounding box alone. The tree is really just one leaf with the diagonal halfspace, as the bounding box coincides with the other sides. So, negative coordinates will always yield "inside", because they are "inside" if you only consider the diagonal halfspace.
If you were to consider the halfspaces that are coplanar with a bounding box corner, edge, or side interesting, then the initial bounding box does not matter, and it would be correct to descend into child five in the case you showed.
Ok, but having like the one I posted above ( insidePolytopeBox ), has a really low cost, so I think I am gonna use it.
I got rid of the BoundingBox and I am calculating the boxes locally.. The problem is that my tree gets very very big and that is because I find many interesting halfspaces, until I have created many levels...
The any-dimensional code is way too spaghetti for me to post it; I cringe every time I look at it.Quote:
Originally Posted by std10093
Besides, the algorithm does not truly scale to "any" number of dimensions D, because each inner node will need to contain 2D child pointers. On 32-bit systems, a single 28-dimensional inner node takes a gigabyte!
The relationship between the number of dimensions D, and leaf nodes having up to t halfspaces, is such that on 64-bit systems where pointers and doubles are the same size, inner and leaf nodes are about the same size ift ≃ 2D / (D + 1)For 20-dimensional data, an inner node is the same size as a leaf node with about 50,000 halfspaces.
In other words, I strongly believe this algorithm is suitable for small number of dimensions only.
Also, I haven't had time to implement a nice convex hull algorithm, either; partly because the methods get pretty hairy (to fully grok) when there are more than three spatial dimensions. (I personally would need to take the time to develop the two and three dimensional cases in parallel first, then compare the two to develop the generic case.)
I think I'll just clean up my 3D convex hull case, post it here, and call it done.