How would YOU explain to someone the concept of the quadtrees?
PS - I know about google too
How would YOU explain to someone the concept of the quadtrees?
PS - I know about google too
Code - functions and small libraries I use
It’s 2014 and I still use printf() for debugging.
"Programs must be written for people to read, and only incidentally for machines to execute. " —Harold Abelson
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
“Salem Was Wrong!” -- Pedant Necromancer
“Four isn't random!” -- Gibbering Mouther
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.
Code - functions and small libraries I use
It’s 2014 and I still use printf() for debugging.
"Programs must be written for people to read, and only incidentally for machines to execute. " —Harold Abelson
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...
Code - functions and small libraries I use
It’s 2014 and I still use printf() for debugging.
"Programs must be written for people to read, and only incidentally for machines to execute. " —Harold Abelson
Technically, a quadtree is just a structure that references 4 other structures of the same type.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:
Code:class quadtree { quadtree* quadrants [4]; void deepen () { int i = 4; while (i -- > 0) quadrants [i] = new quadtree; } };
Last edited by Yarin; 07-14-2013 at 11:01 AM.
I have understood the quadrant concept
Well, as I said in link #3, I want to implement the algorithm in page 580
Code - functions and small libraries I use
It’s 2014 and I still use printf() for debugging.
"Programs must be written for people to read, and only incidentally for machines to execute. " —Harold Abelson
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 :/
Code - functions and small libraries I use
It’s 2014 and I still use printf() for debugging.
"Programs must be written for people to read, and only incidentally for machines to execute. " —Harold Abelson
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.
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.
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
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?
Code - functions and small libraries I use
It’s 2014 and I still use printf() for debugging.
"Programs must be written for people to read, and only incidentally for machines to execute. " —Harold Abelson
> (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
Code - functions and small libraries I use
It’s 2014 and I still use printf() for debugging.
"Programs must be written for people to read, and only incidentally for machines to execute. " —Harold Abelson
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.
The obvious choice is to define the polytope using the vertices, as they are enough to define the faces for a convex polytope.
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.
Last edited by Nominal Animal; 07-21-2013 at 09:34 AM.
Thanks for all your answers Nomimal and yes, it helps². (brainfarts made me laugh
Well, the algorithm is recursive, thus I guess in every moment of the recursion the ε-approximation will be ok, as you mentioned.
What does this code do?
I agree. By setting the z parameter to z, I jump into 2D easily.
Where the methods are in the link planes you provided before.
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.
Code - functions and small libraries I use
It’s 2014 and I still use printf() for debugging.
"Programs must be written for people to read, and only incidentally for machines to execute. " —Harold Abelson
It picks the correct octant to descend into.What does this code do?
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!
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;
Code - functions and small libraries I use
It’s 2014 and I still use printf() for debugging.
"Programs must be written for people to read, and only incidentally for machines to execute. " —Harold Abelson
