BSP split plane determination

Published October 09, 2018
Advertisement

It's been a while since my last blog entry, but the problem I posed on my earlier blog entry still persists - how to efficiently choose a good split plane for an n-vector data structure. To summarize the structure, geographic points are stored as n-vectors (unit vectors) in a binary tree. Branch nodes of this tree define a plane that bisects/halves the unit sphere - one child of the branch contains points that are below the plane, the other child contains points that are above. Like all trees, leaf nodes are turned into branch nodes as they are filled beyond a threshold. That is the basic idea of it.

 

Split plane determination must occur when a leaf becomes a branch or when a branch becomes unbalanced. In either case, the method called to determine the split plane is agnostic as to why it was called - it simply receives a bunch of points that it must determine a good split for. My initial implementation was naive:

 


bestsplit = ...
bestscore = infinite;
for (i=0; i<numpoints; i++) {
  for (j=i+1; j<numpoints; j++) {
    split = (points[i]-points[j]).normalize();
    score = 0;
    for (k=0; k<numpoints; k++) {
      dot = split.dot(points[k]);
      // accumulate some scoring heuristics, such as is each point above/below the split, how far above/below, etc.
      score += ...;
    }
    if ( bestscore > score ) {
      bestscore = score;
      bestsplit = split;
    }
  }
}

So basically - for all points, determine the point between it and every other point, normalize it, and consider this difference as the normal to the split plane, then test it to see how it "scores" using some heuristics. Probably you can see how this will perform miserably for any significant number of points. The big-O notation is something like n^3 which doesn't include the expensive normalization that is called n^2 times. Several other methods were tested, such as fixing the normal of the split plane to be perpendicular to the x, y, or z axis, but these also proved too expensive and/or also had test cases where the split determination was unsatisfactory.

 

Heuristics

Enter calculus. If we can represent each of the heuristics as a mathematical function, we can determine when the function reaches what is called a "critical point". Specifically we are interested in the critical point that is the global maximum or minimum, depending on which heuristic it is for. So far we have three of these.

1. Similar Distance

We don't want a split plane where, for example, all points above are nearby and all points below are far away. The points should be as evenly distributed as possible on either side. Given that the dot product of a plane normal and any other point is negative when the point is below the plane and positive when the point is above, and the absolute value of the dot product increases as distance to the plane increases, the sum of all dot products for a good split will be at or close to zero. If we let \(P\) be the array of points, \(N\) be the number of points in the array, and \(S\) the split plane, the following function will add all dot products:

\(SumOfDots = \displaystyle\sum_{i=1}^{N} P_i \cdot S\)

The summation here is not really part of a mathematical function, at least not one we can perform meaningful calculus on, since the calculation must be done in the code. We don't know ahead of time how many points there will be or what their values are, so the function should be agnostic in this regard. As written we cannot use it without inordinate complication, but consider that it is really doing this:

\(SumOfDots = \displaystyle\sum_{i=1}^{N} P_{ix} * S_x + P_{iy} * S_y + P_{iz} * S_z\)

This summation expanded will look something like:

\(SumOfDots = P_{1x} * S_x + P_{1y} * S_y + P_{1z} * S_z + P_{2x} * S_x + P_{2y} * S_y + P_{2z} * S_z + P_{3x} * S_x + P_{3y} * S_y + P_{3z} * S_z + \cdots\)

We can then rewrite this as:

\(SumOfDots = S_x*(P_{1x} + P_{2x} + P_{3x} + \cdots) + S_y*(P_{1y} + P_{2y} + P_{3y} + \cdots) + S_z*(P_{1z} + P_{2z} + P_{3z} + \cdots)\)

Or similarly:

\(SumOfDots = S_x*(\displaystyle\sum_{i=1}^{N} P_{ix}) + S_y*(\displaystyle\sum_{i=1}^{N} P_{iy}) + S_z*(\displaystyle\sum_{i=1}^{N} P_{iz})\)

As far as the mathematical function is concerned, the sums are constants and we can replace them with single characters to appear concise:

\(A = \displaystyle\sum_{i=1}^{N} P_{ix}\)

\(B = \displaystyle\sum_{i=1}^{N} P_{iy}\)

\(C = \displaystyle\sum_{i=1}^{N} P_{iz}\)

We can pre-calculate these in the code as such:


double A = 0, B = 0, C = 0;
for (i=0; i<N; i++) {
  A += points[i].x;
  B += points[i].y;
  C += points[i].z;
}

We will rewrite the function with these constants:

\(SumOfDots = S_x*A + S_y*B + S_z*C\)

This is great so far. We are interested in when this function reaches zero. To make it simpler, we can square it which makes the negative values positive, and then we become interested in when this function reaches a global minimum:

\(SquaredSumOfDots = (S_x*A + S_y*B + S_z*C)^2\)

So again, when this function reaches zero it means that the points on either side of the split plane - as denoted by \(S\) - are spread apart evenly. This does not mean that \(S\) is a good split plane overall, as the points could all lie on the plane, or some other undesirable condition occurs. For that we have other heuristics.

As a final note, since the vector \(S\) represents a unit normal to a plane, any determination of it must be constrained to space of the unit sphere:

\(S_x^2 + S_y^2 + S_z^2 = 1\)

 

2. Large Distance

Practically speaking, the points will not be random and will originate from a grid of points or combination of grids. But random or not, if the points form a shape that is not equilateral - in other words if they form a rectangular shape instead of a square or an ellipse instead of a circle - the larger axis of the shape should be split so that child areas are not even less equilateral. To achieve this we can emphasize that the sum of the absolute value of all the dot products is large, meaning that the points are farther away from the split plane. To do this, we want to find the global maximum of a function that calculates the sum:

\(SumOfAbsoluteDots=\displaystyle\sum_{i=1}^{N} |P_i \cdot S|\)

Unfortunately there is no way, that I know of, to handle absolute values and still determine critical points, so we need to rewrite this function without the absolute value operator. Really all we are interested in is when this function reaches a maximum, so we can replace it with a square:

\(SumOfSquaredDots=\displaystyle\sum_{i=1}^{N} (P_i \cdot S)^2\)

Not unlike the previous heuristic function, we need to rewrite this so that \(S\) is not contained in the summation, and we can extract the constant values. If we skip some of the expanding and reducing the squared dot product we will arrive at this step:

\(SumOfSquaredDots=(S_x^2*\displaystyle\sum_{i=1}^{N} P_{ix}^2) + (S_y^2*\displaystyle\sum_{i=1}^{N} P_{iy}^2) + (S_z^2*\displaystyle\sum_{i=1}^{N} P_{iz}^2) + (S_x * S_y * 2 * \displaystyle\sum_{i=1}^{N} P_{ix}*P_{iy})+ (S_x * S_z * 2 * \displaystyle\sum_{i=1}^{N} P_{ix}*P_{iz})+ (S_y * S_z * 2 * \displaystyle\sum_{i=1}^{N} P_{iy}*P_{iz})\)

Again we will create some named constants to be concise:

\(D = \displaystyle\sum_{i=1}^{N} P_{ix}^2\)

\(E = \displaystyle\sum_{i=1}^{N} P_{iy}^2\)

\(F = \displaystyle\sum_{i=1}^{N} P_{iz}^2\)

\(G = 2 * \displaystyle\sum_{i=1}^{N} P_{ix}*P_{iy}\)

\(H = 2 * \displaystyle\sum_{i=1}^{N} P_{ix}*P_{iz}\)

\(I = 2 * \displaystyle\sum_{i=1}^{N} P_{iy}*P_{iz}\)

As before, these can be pre-calculated:


double D = 0, E = 0, F = 0, G = 0, H = 0, I = 0;
for (i=0; i<N; i++) {
  D += points[i].x * points[i].x;
  E += points[i].y * points[i].y;
  F += points[i].z * points[i].z;
  G += points[i].x * points[i].y;
  H += points[i].x * points[i].z;
  I += points[i].y * points[i].z;
}
G *= 2.0;
H *= 2.0;
I *= 2.0;

And then the function becomes:

\(SumOfSquaredDots=(S_x^2*D) + (S_y^2*E) + (S_z^2*F) + (S_x * S_y * G)+ (S_x * S_z * H) + (S_y * S_z * I)\)

Again when this function maximizes, the points are farthest away from the split plane, which is what we want.

3. Similar number of points

A good split plane will also have the same number of points on either side. We can again use the dot product since it is negative for points below the plane and positive for points above. But we cannot simply sum the dot products themselves, since a large difference for a point on one side will cancel out several smaller differences on the other. To account for this we normalize the distance to be either +1 or -1:

\(SumOfNormalizedDots=\displaystyle\sum_{i=1}^{N} \frac{P_i \cdot S}{\sqrt{(P_i \cdot S)^2}}\)

Expanding this becomes:

\(SumOfNormalizedDots=\displaystyle\sum_{i=1}^{N} \frac{P_{ix} * S_x + P_{iy} * S_y + P_{iz} * S_z}{\sqrt{(S_x^2*P_{ix}^2) + (S_y^2*P_{iy}^2) + (S_z^2*P_{iz}^2) + (S_x*S_y*2*P_{ix}*P_{iy})+ (S_x*S_z*2*P_{ix}*P_{iz})+ (S_y*S_z*2*P_{iy}*P_{iz})}}\)

Unfortunately there is no way to reduce this function so that we can extract all \(S\) references out of the summation, and as with the previous heuristics, put all \(P\) references into constants that we pre-calculate and use to simplify the function. So at present, this heuristic is not able to be used. I am still working on it. The problem is that each iteration of the sum is dependent on the values of \(S\) in such a way that it cannot be extracted.

Putting it all together

What we ultimately want to do is combine the heuristic functions into one function, then use this one function find the critical point - either the global minimum or global maximum depending on how we combine them. The issue is that for \(SquaredSumOfDots\) we want the global minimum and for \(SumOfSquaredDots\) we want the global maximum. We can account for this by inverting the former so that we instead want the global maximum for it. The combined function then becomes:

\(Combined = SumOfSquaredDots - SquaredSumOfDots\)

Applying the terms from each function, we get:

\(Combined = (S_x^2*D) + (S_y^2*E) + (S_z^2*F) + (S_x * S_y * G)+ (S_x * S_z * H) + (S_y * S_z * I) - (S_x*A + S_y*B + S_z*C)^2\)

Expanding the square on the right side and combining some terms will give us:

\(Combined = (S_x^2*(D-A^2)) + (S_y^2*(E - B^2)) + (S_z^2*(F - C^2)) + (S_x * S_y * (G - (2 * A * B)))+ (S_x * S_z * (H - (2 * A * C))) + (S_y * S_z * (I - (2 * B * C)))\)

Again we can combine/pre-calculate some constants to simplify:

\(J = D-A^2\)

\(K = E-B^2\)

\(L = F-C^2\)

\(M = G - (2 * A * B)\)

\(N = H - (2 * A * C)\)

\(O = I - (2 * B * C)\)


double J = D - (A * A);
double K = E - (B * B);
double L = F - (C * C);
double M = G - (2 * A * B);
double N = H - (2 * A * C);
double O = I - (2 * B * C);

And then apply these to the function:

\(Combined = (S_x^2*J) + (S_y^2*K) + (S_z^2*L) + (S_x * S_y * M)+ (S_x * S_z * N) + (S_y * S_z * O)\)

Because we are lazy, we then plug this into a tool that does the calculus for us. And this is where I am currently blocked on this issue, as I have found no program that can perform this computation. I've actually purchased Wolfram Mathematica and attempted it with the following command:

Maximize[{((x^2)*j) + ((y^2)*k) + ((z^2)*l) + (x*y*m) + (x*z*n) + (y*z*o), ((x^2) + (y^2) + (z^2)) == 1}, {x, y, z}]

After 5 or 6 days this had not finished and I had to restart the computer it was running on. I assumed that if it took that long, it would not complete in a reasonable amount of time. I will update this blog entry if I make any progress on this problem as well as the 3rd heuristic function.

 

Further Optimizations

While I have not gotten this far yet, it may be ultimately necessary to emphasize (or de-emphasize) one heuristic over the others in order to further optimize the split determination. This could be done by simply multiplying each heuristic function by a scalar value to either increase or decrease it's emphases on the final result. The reason I haven't researched this yet is because if I cannot find the global maximum without these extra variables, I certainly cannot do it with them.

0 likes 6 comments

Comments

lawnjelly

Blimey! Are you still working on this? I thought you'd have had it solved that afternoon... hehe :)

I still have an inkling feeling you are barking up the wrong tree(!) with a BSP tree, for exactly the kind of reasons you are finding, difficulty in building / maintaining etc.

As well as my original suggestions from your last blog post I'd like to also suggest Delaunay triangulation as another useful structure for prestoring the nearest neighbours, and allowing fast recursion out. You could also use Delaunay triangles for rough visualization. 

Also 'gridding' up a sphere would be useful in a number of your cases, if only for construction of other data structures, but maybe for fast location too. So something that can tessellate to a sphere would be good.. Blender has an 'icosphere' made up of triangles so that might be worth looking into.

October 10, 2018 07:28 AM
JoeJ

To find a good split plane i would do this:

Calculate average line that fits the point set.

Sort points along line to find split which divides in two equal sized sets (or accept unbalanced tree and use median or center of mass).

 

Average line was discussed here: https://www.gamedev.net/forums/topic/698600-how-close-are-pointa-forming-a-line/?tab=comments#comment-5389006

(To make my algorithm work in 3D using band 3 of spherical harmonics should work, but something like least squares regression should work as well)

 

However, i would prefer octree over BSP to avoid those worries.

 

October 10, 2018 09:00 AM
lawnjelly

Here's a nice article on some solutions:

https://en.wikipedia.org/wiki/Nearest_neighbor_search

The one I suggested storing the proximity information on the nodes is called in the article 'proximity graph methods'. BSP (presumably arbitrary plane) is one of their suggested space partitioning methods, although I'm not convinced myself :) , I too would be more inclined for a regular partitioning of some kind, but am willing to be converted! Perhaps BSP might end up with slices like Delaunay triangulation anyway.

October 10, 2018 10:12 AM
jorgander
On 10/10/2018 at 3:28 AM, lawnjelly said:

Blimey! Are you still working on this? I thought you'd have had it solved that afternoon... hehe :)

I still have an inkling feeling you are barking up the wrong tree(!) with a BSP tree, for exactly the kind of reasons you are finding, difficulty in building / maintaining etc.

As well as my original suggestions from your last blog post I'd like to also suggest Delaunay triangulation as another useful structure for prestoring the nearest neighbours, and allowing fast recursion out. You could also use Delaunay triangles for rough visualization. 

Also 'gridding' up a sphere would be useful in a number of your cases, if only for construction of other data structures, but maybe for fast location too. So something that can tessellate to a sphere would be good.. Blender has an 'icosphere' made up of triangles so that might be worth looking into.

One of the reasons why I want to stick with a BSP tree is that queries are fast and simple... one test per level of hierarchy using just a dot product. Also since the data points are represented by unit vectors, and querying compares these against another unit vector using the dot product (e.g. to find the closest N points to another arbitrary point), the implementation translates well to graphics hardware, specifically a pixel shader. I've already got a shader implemented for calculating and displaying land/coastlines - rendering these at 2560x1440 on the CPU takes a few minutes but on the GPU takes a tenth of a second. So I'm hoping the difficult in building/maintaining pays off down the road. But we'll see.

 

Yeah this is my hobby project, so it takes a back seat to other things.

October 11, 2018 08:16 PM
lawnjelly

Maybe you can post a screenshot of the existing version? :)

The end use also doesn't seem to be clear (to me at least lol), you want a fast structure for building interpolated data at points on the earth, but is this ultimately for visualization (like google earth)? Are you writing the interpolation to texture maps, or intending to do the calculation for each screen pixel individually (more like a ray tracing approach)? 

October 12, 2018 07:34 AM
jorgander
On 10/12/2018 at 3:34 AM, lawnjelly said:

Maybe you can post a screenshot of the existing version? :)

The end use also doesn't seem to be clear (to me at least lol), you want a fast structure for building interpolated data at points on the earth, but is this ultimately for visualization (like google earth)? Are you writing the interpolation to texture maps, or intending to do the calculation for each screen pixel individually (more like a ray tracing approach)? 

With some help I've been able to work through the math part and am implementing it now. Apparently there is a another way to determine the global maximum of quadratic form mathematical functions, and the one I have now is a quadratic. The details are in this mathematica stackexhange post I created. If it works it will be quite significant, as it means that except for the loop to build the constant values, the time to determine split planes is constant regardless of the number of points. But there are potential problems with any solution:

  • For the constant (summation) values and ensuing math, I am using the {{long double}} data type. I don't think this will overflow for any practical uses of the implementation (it would have to mean billions or trillions of points, I am guessing), but I suppose it is possible.
  • If I ever am able to combine the 3rd heuristic function, it might change the form of the function from quadratic and I would have to revisit the math part.

To answer your questions:

Yes, the the end use is visualization of animated weather data over a set of parameters. For example, rainfall between now and 7 days from now, for a particular geographic region, etc. Again, the idea is that if I can render weather data at one point in time fast enough (say, less than a second), animating over time will be feasible. It won't run at the FPS of video games, but fast enough to look pretty.

As the data structure will be pushed up to GPU memory and used by a pixel shader, the calculation will be per-pixel which will create images/animation at the appropriate resolution for the screen it's running on. To compose images/animation for other resolutions, the implementation could be modified to instead use a compute shader, which is not that different and in fact it's what I though of first - instead of writing to the frame buffer, we would instead write to a texture. Simple enough.

October 14, 2018 05:24 PM
You must log in to join the conversation.
Don't have a GameDev.net account? Sign up!
Profile
Author
Advertisement
Advertisement