Bug: celltree nodes sometimes contain incorrect bounding planes

[Huite]

I ran into some issues with my numba port of this cell tree implementation. It would not locate points in a cell, while the point was clearly located inside of the cell (by visual inspection). I checked, and the C++ implementation here has the same issue.

I did some digging, and it turns out to be a small bug in the construction of the tree here for Lmax and Rmin: https://github.com/NOAA-ORR-ERD/cell_tree2d/blob/5d72c595249c198986ea9d6e97b3f69be4111de6/src/cell_tree2d.cpp#L305-L306

It can happen that the values of Lmax and Rmin are not monotonically increasing with the buckets. This can happen for example when a grid contains large and small cells. See this visual example:

The centroid of 1 is left of the centroid of 68. So if two buckets contain one cell each, 1 is in the first, and 68 is in the second. In this case, I had four buckets in total during tree construction. Plane was 2, so in the code quoted above, the Lmax of bucket 1 is taken. But Lmax of bucket 0 is to the right!

This is easily fixed by taking the maximum Lmax value of all buckets left of the split, and the minimum Rmin value of all buckets right of the split.

Something like this:

double Lmax = std::numeric_limits<double>::min();
double Rmin = std::numeric_limits<double>::max();

for (unsigned int b = 0; b < plane; b++) {
  buk = buks.at(b);
  if (buk.Lmax < Lmax) {
      Lmax = buk.Lmax;
  }
}

for (unsigned int b = plane; b < buks.size(); b++) {
  buk = buks.at(b);
  if (buk.Rmin < Rmin) {
      Rmin= buk.Rmin;
  }
}

root.Lmax = Lmax;
root.Rmin = Rmin;

Cross-ref: https://github.com/Deltares/numba_celltree/issues/2

I've put some example data which triggers the problem there as well.

[jay-hennen]

Thank you very much for the issue report, and for locating a solution! I'll be working on getting this fixed shortly.

And a huge thank you for extending the functionality to line intersection and Numba. It's heartening to see the code I wrote all those years ago being improved upon.

As an aside, when I wrote this I did it as a re-implementation of the algorithm directly from the paper. References to the published code and the VTKCellTreeLocator were slim to none; though I was aware of both at the time, I could not understand them until I had pretty much done it myself. It's also worth pointing out that the first release was only intended for use on triangular (exclusive) grids, and the polygonal support was added years later.

I bring this up because at the time I tried a lot of awful triangular grids and failed to find a (valid) one that would cause a problem like this. If you CAN do this incidentally, I would be extremely interested in seeing such a topology.

[Huite]

For completeness sake, I included the (voronoi) grid that caused this problem for me (also posted in the cross ref'ed issue):

```python def voronoi_grid(): vertices = np.array( [ [-0.08166667, 0.877], [-0.06966667, 0.88], [-0.07566667, 0.88233333], [-0.06266667, 0.88166667], [-0.04966667, 0.88333333], [0.01066667, 0.888], [-0.06, 0.88633333], [-0.04566667, 0.88866667], [-0.08933333, 0.87566667], [-0.03733333, 0.88633333], [0.016, 0.89366667], [0.021, 0.893], [-0.006, 0.89066667], [0.002, 0.89433333], [-0.03466667, 0.89566667], [-0.06633333, 0.90066667], [-0.02366667, 0.89133333], [-0.01633333, 0.89866667], [-0.061, 0.906], [-0.08333333, 0.90066667], [-0.053, 0.90266667], [-0.07733333, 0.90366667], [0.00733333, 0.9], [0.005, 0.90866667], [-0.065, 0.91466667], [0.028, 0.91333333], [0.02266667, 0.90933333], [-0.044, 0.90666667], [0.01466667, 0.91466667], [-0.074, 0.92266667], [-0.02066667, 0.915], [-0.04633333, 0.91766667], [-0.05933333, 0.92233333], [0.00466667, 0.92733333], [-0.03566667, 0.92333333], [-0.009, 0.91966667], [-0.06466667, 0.92633333], [-0.00566667, 0.92966667], [-0.04166667, 0.935], [-0.00933333, 0.93833333], [-0.017, 0.94466667], [-0.05466667, 0.94966667], [-0.031, 0.943], [-0.081, 0.955], [-0.07866667, 0.95966667], [-0.04066667, 0.95133333], [-0.03333333, 0.96], [-0.04, 0.96633333], [-0.04866667, 0.96866667], [-0.09366667, 0.96833333], [-0.08966667, 0.971], [-0.05833333, 0.96233333], [-0.067, 0.965], [-0.09133333, 0.97733333], [-0.10433333, 0.98233333], [-0.08066667, 0.97633333], [-0.08466667, 0.98233333], [-0.06966667, 0.977], [-0.099, 0.98366667], [-0.052, 0.98266667], [-0.06266667, 0.986], [-0.09933333, 0.988], [-0.09766667, 0.991], [-0.048, 0.98666667], [-0.049, 0.995], [-0.06666667, 0.996], [-0.08933333, 0.99], [-0.03966667, 0.996], [-0.10533333, 0.99666667], [-0.057, 1.001], [-0.03666667, 1.002], [-0.08866667, 0.99866667], [-0.09666667, 1.00433333], [-0.07, 1.00166667], [-0.07866667, 1.00333333], [-0.082, 1.01366667], [-0.06333333, 1.016], [-0.075, 1.01733333], [-0.05766667, 1.02], [-0.09033808, 0.87758482], [-0.08933333, 0.872], [-0.07665982, 0.88451826], [-0.06037555, 0.88914993], [-0.04945385, 0.8935359], [-0.05186486, 0.89585586], [-0.06711625, 0.89816133], [-0.08303333, 0.89856667], [-0.08733333, 0.89866667], [-0.07763333, 0.90576667], [-0.069, 0.91466667], [-0.07354667, 0.92017333], [-0.07596154, 0.92419231], [-0.06489888, 0.92877154], [-0.04946667, 0.9376], [-0.0564, 0.9462], [-0.05843218, 0.95301379], [-0.0569668, 0.95720885], [-0.0807439, 0.95269512], [-0.087, 0.955], [-0.0825, 0.9635], [-0.08511261, 0.97315766], [-0.0857451, 0.97001961], [-0.09366667, 0.965], [-0.097, 0.96833333], [-0.09467296, 0.97918868], [-0.09796907, 0.98134708], [-0.10433333, 0.98], [-0.10609231, 0.98333846], [-0.10323333, 0.9893], [-0.10605436, 0.99549499], [-0.10495967, 0.9978624], [-0.100392, 1.003656], [-0.09764359, 1.00791538], [-0.08851538, 1.01544359], [-0.075, 1.021], [-0.05776437, 1.02166092], [-0.052, 1.02], [-0.06204241, 1.01286489], [-0.06533333, 1.00633333], [-0.057, 1.005], [-0.03666667, 1.005], [-0.031, 1.002], [-0.03533333, 0.99166667], [-0.04286036, 0.9829955], [-0.04837387, 0.97759009], [-0.04772973, 0.97428829], [-0.03788839, 0.97224585], [-0.02424138, 0.96389655], [-0.01386667, 0.95093333], [-0.00304, 0.94305333], [0.00546667, 0.9276], [0.00620513, 0.92702564], [0.01379977, 0.92203527], [0.02983333, 0.91516667], [0.03234359, 0.91085128], [0.02514201, 0.9033925], [0.02297315, 0.89581879], [0.02360177, 0.89002655], [0.01121622, 0.8847027], [-0.0065451, 0.88358039], [-0.02319655, 0.88334138], [-0.03472165, 0.88045704], [-0.04998995, 0.8802621], [-0.0595, 0.8785], [-0.06966667, 0.876], [-0.08069072, 0.87480412], ] ) faces = np.array( [ [13, 5, 10, 22, -1, -1, -1, -1, -1], [27, 14, 16, 17, 30, 34, 31, -1, -1], [17, 12, 13, 22, 23, 35, 30, -1, -1], [24, 18, 20, 27, 31, 32, -1, -1, -1], [34, 30, 35, 37, 39, 40, 42, 38, -1], [66, 56, 55, 57, 60, 65, 73, 74, 71], [60, 59, 63, 64, 69, 65, -1, -1, -1], [80, 8, 79, -1, -1, -1, -1, -1, -1], [79, 8, 0, 2, 81, -1, -1, -1, -1], [2, 1, 3, 6, 82, 81, -1, -1, -1], [6, 4, 7, 83, 82, -1, -1, -1, -1], [83, 7, 9, 14, 27, 20, 84, -1, -1], [15, 85, 84, 20, 18, -1, -1, -1, -1], [86, 85, 15, 21, 19, -1, -1, -1, -1], [87, 86, 19, -1, -1, -1, -1, -1, -1], [87, 19, 21, 88, -1, -1, -1, -1, -1], [88, 21, 15, 18, 24, 89, -1, -1, -1], [89, 24, 32, 36, 29, 90, -1, -1, -1], [90, 29, 91, -1, -1, -1, -1, -1, -1], [91, 29, 36, 92, -1, -1, -1, -1, -1], [36, 32, 31, 34, 38, 93, 92, -1, -1], [93, 38, 42, 45, 41, 94, -1, -1, -1], [94, 41, 95, -1, -1, -1, -1, -1, -1], [96, 95, 41, 45, 46, 47, 48, 51, -1], [43, 97, 96, 51, 52, 44, -1, -1, -1], [97, 43, 98, -1, -1, -1, -1, -1, -1], [98, 43, 44, 99, -1, -1, -1, -1, -1], [99, 44, 52, 57, 55, 100, -1, -1, -1], [50, 101, 100, 55, 56, 53, -1, -1, -1], [49, 102, 101, 50, -1, -1, -1, -1, -1], [102, 49, 103, -1, -1, -1, -1, -1, -1], [103, 49, 50, 53, 104, -1, -1, -1, -1], [58, 105, 104, 53, 56, 66, 62, 61, -1], [106, 105, 58, 54, -1, -1, -1, -1, -1], [106, 54, 107, -1, -1, -1, -1, -1, -1], [107, 54, 58, 61, 108, -1, -1, -1, -1], [108, 61, 62, 68, 109, -1, -1, -1, -1], [109, 68, 110, -1, -1, -1, -1, -1, -1], [68, 62, 66, 71, 72, 111, 110, -1, -1], [111, 72, 112, -1, -1, -1, -1, -1, -1], [72, 71, 74, 75, 113, 112, -1, -1, -1], [113, 75, 77, 114, -1, -1, -1, -1, -1], [77, 76, 78, 115, 114, -1, -1, -1, -1], [78, 116, 115, -1, -1, -1, -1, -1, -1], [76, 117, 116, 78, -1, -1, -1, -1, -1], [74, 73, 118, 117, 76, 77, 75, -1, -1], [73, 65, 69, 119, 118, -1, -1, -1, -1], [64, 67, 70, 120, 119, 69, -1, -1, -1], [70, 121, 120, -1, -1, -1, -1, -1, -1], [67, 122, 121, 70, -1, -1, -1, -1, -1], [63, 123, 122, 67, 64, -1, -1, -1, -1], [124, 123, 63, 59, -1, -1, -1, -1, -1], [52, 51, 48, 125, 124, 59, 60, 57, -1], [48, 47, 126, 125, -1, -1, -1, -1, -1], [46, 127, 126, 47, -1, -1, -1, -1, -1], [45, 42, 40, 128, 127, 46, -1, -1, -1], [39, 129, 128, 40, -1, -1, -1, -1, -1], [37, 33, 130, 129, 39, -1, -1, -1, -1], [131, 130, 33, -1, -1, -1, -1, -1, -1], [35, 23, 28, 132, 131, 33, 37, -1, -1], [28, 26, 25, 133, 132, -1, -1, -1, -1], [134, 133, 25, -1, -1, -1, -1, -1, -1], [135, 134, 25, 26, -1, -1, -1, -1, -1], [22, 10, 11, 136, 135, 26, 28, 23, -1], [137, 136, 11, -1, -1, -1, -1, -1, -1], [5, 138, 137, 11, 10, -1, -1, -1, -1], [139, 138, 5, 13, 12, -1, -1, -1, -1], [140, 139, 12, 17, 16, -1, -1, -1, -1], [9, 141, 140, 16, 14, -1, -1, -1, -1], [4, 142, 141, 9, 7, -1, -1, -1, -1], [3, 143, 142, 4, 6, -1, -1, -1, -1], [144, 143, 3, 1, -1, -1, -1, -1, -1], [0, 145, 144, 1, 2, -1, -1, -1, -1], [80, 145, 0, 8, -1, -1, -1, -1, -1], ] ) return vertices, faces ```

This grid causes a problem here:

vertices, faces = voronoi_grid()
tree = CellTree2d(vertices, faces, -1)
tree.locate_points([[-0.02, 0.90]])  # should return 1, but returns -1

Interestingly, I haven't run into issues with triangular grids at all as far as I'm aware of. Only when taking a look at these voronoi grids did I observe odd gaps when rasterizing them.

I think this is due to the limited number of neighbors of triangular cells. My gut feeling is that this causes less abrupt shifts between cell sizes during construction of a node. If you were looking to trigger it for triangles, you want to set the number of buckets during construction high, and (I think) the number of leaves as well. That should increase the probability that a triangle in one of the buckets exceeds Lmax or Rmin of the buckets at either side of plane. But it's difficult to purposefully engineer such a failure; if the cells are split earlier in construction, it won't show up.

I've added a validation method just in case: it checks for every node whether its children (or the bounding boxes for leaf nodes) fall inside of its bounding box (which can be constructed using the total bounding box of the tree, and then progressively taking Lmax or Rmin as you traverse the tree).

I recognize the understanding part: only after finishing the porting (very fortunate that it's public domain!) did I have a decent grasp of the algorithm! In terms of further developments, I've also implemented barycentric interpolation, and I'll probably extend to 3D at some point. I'm a little conflicted about numba part currently: I had to do quite some twiddling to match C++ performance and the nature of the code is that I won't get my "casual" Python colleagues interested anyway. Distribution is easy (no static binaries), but does require an initial JIT compile (which can be cached fortunately). On the other hand, I can disable JIT entirely and seamlessly debug in Python!

And I still want to see if the benchmarks from the paper hold up versus commonly available spatial indexes in Python!

[jay-hennen]

I think this is due to the limited number of neighbors of triangular cells. My gut feeling is that this causes less abrupt shifts between cell sizes during construction of a node.

You know, I think you are pretty much correct on this one. The example that keeps floating up from my memories are fractal-likes such as the Sierpinski triangles. However, then I also remember that such examples are NOT valid 3-neighbor triangular topologies anyway so of course they won't work. But on the other hand, they could be considered valid multi-poylgon topologies! Guess it makes sense that you encountered this bug then!

Just want to say again...hugely impressed you found this bug and the solution.

On the other hand, I can disable JIT entirely and seamlessly debug in Python!

That sounds lovely! Python's debugging and introspection tools are a luxury I'm not sure I could give up anymore. I have to admit that the difficulty of working on the system is part of why it hasn't been developed more.

As for a performance, here's the summary of some comparative testing I did (on a Scipy poster) from many years ago. The VTK implementation is still the gold standard. I believe it is still limited to triangular/trapezoidal cells, but is 3D and is shockingly fast and light. However, it is also stuck in the VTK library and needs VTK faces/vertex constructs.

If I didn't have to shelve this project and move on to other things I expect it would have begun to take lessons from the VTK version in the pursuit of speed. For example the VTK version uses raw C arrays where cell_tree2d uses C++ vectors.

It is worth pointing out that at the time the only other true bounding volume hierarchy in a python library I could find was that TrapezoidMapTriFinder from matplotlib, and as you can see it was disappointing (but actually functional!). The R-Trees didn't have the speed, nor did they provide the actual index of the target cell without post-processing, a major feature deficit.