# Finding a bounding volume (line segments) from a kDop definition.

I have a problem. I'm trying to recover a bounding volume (actually line segments that form the bounding volume) from a kDop definition (in a 3D space). (its to draw the kDop on screen)

In my kDOP structure i have the Min/Max values calculated for each axis. (And well, i know the axis used)

I tried before coming here 3 ways to resolve it :

1) Intersecting 3 neighboring planes to find a point and then connecting them (Sometimes the intersection is out of the k-DOP, i didn't manage to find a way to figure the correct planes to use)

2) Starting from the lines of the AABB (DOP6) and then intersecting the lines with each of the k-DOP planes, but when a line is completely behind a plane (i.e. not visible), we have to add two or more new lines and i don't know which is to be drawed and which not. :P (and i get too many lines drawed)

3) Brute force method. Intersect each plane with each other plane to get a set of lines and then intersect these lines with each plane to find the correct line segment (or discard the line) (too slow :()

What's the best way to go? Can you give some explanation on it?

Thanks

EDIT : About the slowness. Well yes, i need to do it for 100-120 objects under 16,667 ms (n^3 method takes 0,96 ms (and i've optimized pretty much everything))

EDIT 2: Answering the questions from TonyK :

1) (1,0,0)(0,1,0)(0,0,1) // AABB (Axis Aligned Bounding Box)

(1,1,1)(1,-1,1)(1,1,-1)(1,-1,-1) // Corners

(1,1,0)(1,0,1),(0,1,1),(1,-1,0),(1,0,-1),(0,1,-1) // Edges

These are the axes used. Are fixed and I can happen to use a subset of those. (only AABB (DOP6) or AABB + Corners (DOP14) or AABB+Corners+Edges (DOP26)). Yes, the only thing changing is the Min/Max values. (The normals for each plane are prestored (26 normals in the worst case (DOP26)))

2) My point of view is the camera of the simulation and "can" change. (the line segment obtained will be multiplied by World/View/Projection matrix to get the final position on the Viewport, but that's a GPU detail)

3) Mhn...well, there are Axes that do not intersect each other.

4) What do you mean by "treating the bounding plane as opaque"?

EDIT 3:

Regarding the other question, No, I treat them as wireframe so the lines under the kDOP surface will appear.

-
Just to translate it into purely math. language. You know the N half-spaces whose intersection is a convex polytope. You want to find the edges of this polytope. The trivial $O(N^3)$ algorithm (find the intersection line of each pair of planes and see what segment other half spaces cut on this line) is, apparently, too slow. How to do it faster? – fedja Feb 1 '10 at 14:50
For those like me who have no idea what a kDOP is, the wikipedia page at en.wikipedia.org/wiki/Bounding_volume has a small explanation. – Mariano Suárez-Alvarez Feb 1 '10 at 15:34
Now we're talking milliseconds, perhaps you can fill in some of the details. How many axes are there, and can you give a simple description of them (e.g. symmetry axes of the cube)? Are the axes fixed in space, with only their min/max values changing? Is your point of view at infinity (and is it fixed)? Do all the axes intersect each other? And you are treating the bounding planes as opaque, right? – TonyK Feb 1 '10 at 17:32
In other words, your axes are the symmetry axes of the cube :-) (and my question about whether they all intersect was a stupid one). By the 'bounding planes', I meant the surface of the k-DOP: if you treat them as opaque, then you only draw lines that are not obscured by the k-DOP itself (so a cube would typically require nine lines instead of twelve). – TonyK Feb 1 '10 at 17:58
I thought a little.... Can't i simply precalculate the vertexes of the various edges manually using 3 planes as I have a limited (ALWAYS The same) number of planes/vertex to calculate? Someone has an image with the simmetry axes of the cube? So i can take a look. This way it should be quite faster (I just calculate the right vertex) :P – feal87 Feb 1 '10 at 20:05

OK, it seems like your $N=52$ then and the normals to the planes are fixed. Then you are right that $N^3$ method takes 1ms time. Indeed, to find the parametric coefficients of the intersection of any 2 planes once you prestore the line directions and the inverse matrices takes just 6 multiplications and 3 additions, to find the cut takes 52 cycles of 6 multiplications, 5 additions, one division, and two comparisons each (keeping the max and the min), and doing all 1326 pairs that requires 1326*(9+52*14)=977262 elementary operations, which puts your processor speed at 10^9 multiplications per second, which is about right. There is an $O(N^2)$ algorithm if you know a point inside the polytope, which, once you are talking about kDop, you, probably do, but we need to count operations very carefully to see if it beats the $N^3$ one.

The key problem is that in the worst case when every plane matters, you may have about 3*52=156 edges which is just one tenth of the maximal possible amount and if you spend the same 9+52*14 operations on each edge, you will be just on the edge of what you need. Still, it is worth trying. It works under the assumption that no 4 planes intersect at one point. To make sure that it is the case, just perturb your Mins and Maxes by some random numbers that are much bigger than the precision with which you can compute the points but much less than what is visible to the eye (say, $10^{-4}$ times the object size, also make sure that you decrease Min's and increase Max's) but do all comparisons with much higher precision (like $10^{-12}$ times the object size) where the object size is just the maximal span in your directions, which you can find in linear time.

Start with finding one edge. To this end, just take any point inside the polytope, compute the vertices of the dual polytope with respect to that point, choose any direction, and take the furthest vertex of the dual polytope from the origin. Now, you know that this plane certainly matters. It will take you just linear time, so I'll not even count that. Once you figured out one plane, you can find an edge on that plane in the time 50*(9+52*14)=36850 elementary operations just as in the $N^3$ method. Keep track on which planes give you the boundary points.

Now, at each moment, keep track of edges you already have, and the vertices you have with their planes. For each vertex also keep the track of whether you know each of the adjacent edges (according to the pair of planes) completely. In the beginning, you have the list of 2 vertices with one known and two unknown edges.

Take the first vertex in the list with at least one unknown adjacent edge. You know the pair of planes for this edge, so you can find the new edge in time (9+26*15)=737, where you have 26 instead of 52 because you can find whether you'll need Max or Min in each direction in one comparison. Once you have found the edge, see if its other end is one of the already found vertices and mark the corresponding edge at that vertex as known (about 156*2=312 operations). If not, add its other vertex to the end of the list. This gives you 1049 operations per edge.

Thus, your total is 36850+156*1049=200494 operations, which is 5 times better than my count for the $N^3$ algorithm but still a bit slower than you need.

There is a room for improvement. The most obvious one is that you may want to keep the list of already found vertices for each pair of planes (the update costs 3 operations and reduces 312 to just 6. Thus, you are down to 36850+156*153=154318, which is within the limits, but just barely. Another thing is that you can try to reduce the time you need to find the first edge. but I do not see how to do it at the moment.

All right, first thing first. The main part of the cure is to construct the admissibility table that, for each pair of directions, lists all directions that may form a vertex with the given two. Surprisingly, it is very short. I'll keep the directions in the 0,$\pm 1$ format like you (not normalized). I assume that you keep an array of these directions p of length 26. The only thing I want you to make sure about the order is that p[2k] and p[2k+1] are opposite for each k=0,1,...,12. This will spare a few operations in the edge tracing algorithm.

First, let us exclude all degenerate pairs. That is pretty easy. The pair (A,B) is degenerate iff A+B is a multiple of some other direction in the list. In practical terms, it means that all non-zero entries of $A+B$ are equal in absolute value. The first thing is to create the table of pairs of directions and precompute some information for them. To make it working for both algorithms, I suggest that for each pair (i,j) of indices, you create the structure Q[i,j] with the following fields

Boolean nondegeneracy; true if the pair is non-degenerate, false if it is degenerate. You should get 288 non-degenerate pairs if the order of i,j matters (144 with i

vector normal; normal=$p[i]\times p[j]$ is the direction orthogonal to both p[i] and p[j]

iprojector, jprojector; those are just found as $v=p[j]\times\operatorname{normal}$; iprojector$=v/(v\cdot p[i])$ and similarly for jprojector. This will allow you to find a point X with given scalar products $a=p[i]\cdot X$ and $b=p[j]\cdot X$ as a times iprojector plus b times jprojector (6 multiplications and 3 additions)

an array goodindices of structures. Each structure element goodindices[m] will have 3 fields: integer k such that p[i],p[j],p[k] may make a vertex, and 3 scalar products iscal=iprojector$\cdot$p[k],jscal=jprojector$\cdot$p[k],and nscal=normal$\cdot$p[k]. To construct it, you take your pair a=p[i],b=p[j] of directions and include the structure with the first field k (once you know k, you can compute the other three fields) into goodindices if the following two conditions hold:

1)the determinant det=det(a,b,c) is not 0

2) for every index m different from i,j,k at least one of the determinants det(a,b,d), det(a,d,c), and det(d,b,c) where d=p[m] has sign opposite to the sine of det, i.e., if, say, det > 0, one of them is strictly less than 0 (since all coordinates are integer, compare to -0.5). This checks that there is no direction in the cone formed by a,b,c.

Funnily enough, the total number of indices in goodindices is never greater than 6.

Also, let Max[i] be the maximal scalar product of your object points with p[i] (I prefer to keep all 26 directions and Max array to keeping 13 and Min and Max arrays just to avoid any casework anywhere)

Now, when the table is ready, it is very easy to run the N^3 algorithm. For each pair i,j with i < j do the following.

If Q[i,j].nondegeneracy=false, just skip the pair altogether

Determine tmin as the maximum over all structures S in Q[i,j].goodindices such that S.nscal<0 of the ratios (M[S.k]-M[i]S.iscal-M[j]S.jscal)/S.nscal .

Determine tmax as the minimum over all structures S in Q[i,j].goodindices such that S.nscal>0 of the ratios (M[S.k]-M[i]S.iscal-M[j]S.jscal)/S.nscal.

If tmin < tmax, put X=M[i]Q[i,j].iprojector+M[j]Q[i,j].jprojector.; compute the points X+tmin Q[i,j].normal and X+tmax Q[i,j].normal and draw the edge between them (or store it for later drawing).

You should optimize here updating tmin and tmax when looking at each k and stopping immediately if tmax gets less than or equal to tmin. This cycle over k is the most time-consuming part when defining the edge, so how you write it determines the real speed of this algorithm. Optimize as much as possible here to make sure that you do not do any extraction from array twice, etc. Each little extra operation counts and memory is not a problem. That's why I suggested to prestore the scalar products as well, though it might seem a bit silly when we are talking of cycle of length 6 within a cycle of length 144

That's all. It should run fast enough to be acceptable already and it is robust in the sense that you do not need to perturb your maxes, check for paralellness anywhere, etc., etc.

The other algorithm is still about 4 times faster but I should check that it has no stability issues before posting it. Try this first and tell me the running time.

-
Thanks a lot for the detailed and precise explanation!!! It should be easy to implement having such a good guide. (I'll do it anyway even if I can find another solution just as an educational project) A note, is it possible to do as I said in the last comment of the question and simply precalculate each vertex manually preknowing which planes to intersect (as i don't change EVER the planes) simply by doing a 3 planes intersect and linking the points one to one (it is really tedious to write every combination, but should VERY fast) What do you think? – feal87 Feb 1 '10 at 20:20
N is 26, I think, not 52. So fedja's method should be fast enough. – TonyK Feb 1 '10 at 20:22
(yes its 26. 13 min and 13 max with the normals inverted) – feal87 Feb 1 '10 at 20:24
The whole point is that despite the directions of the planes stay fixed, the vertices are not always formed by the intersections of the same triples (in a striking contrast with $\mathbb R^2$ where you just go counterclockwise) though you are right that there are some triples that have no chance to form vertex (the ones that form a cone inside which you have another direction), so you do not know the combinatorial type in advance. – fedja Feb 1 '10 at 20:42
But, of course, you can do the precomputation to determine for each pair of directions which other directions have any chance to make a vertex with them, which may shave some extra time off in both algorithms. Just do not do it manually: the machine can easily determine if a given vector has positive coefficients with respect to 3 other vectors and construct the table for you. – fedja Feb 1 '10 at 20:50

may I know the order of your planes? when I change the order of the planes, the lines that are drawn changes as well

-
when you create the admissibility table, the order of normal for that table must correspond to the order of normal for the KDOP that is used to compute the min/max values? – user54555 Jul 1 '14 at 19:36