Dear Andrew, I realised I dropped the ball on responding to an earlier query on this. Your implementation looks neat and is along the lines of what I was thinking (that reference you linked is one of the all knn refs I had looked at). However I have to say that I am rather disappointed about the final speed-up. I had honestly expected something in the range of an order of magnitude given consistent locality for many points. Is it possible in the end that the density of points in neuron data is just a bit low? Do you know where the algorithm is spending time right now? Is it very different from the naive implementation of scanning all points with regular knn?

Hi Greg, I think you may have meant to respond atclbarnes/nblast-rs#20, so I'll respond there.

Reverting this to draft status as there are performance pitfalls for 3d data that will need to be looked into.

Thank you very much for this contribution!
This looks really promising and the feature would fit in well with the rstar library.

I'm still trying to understand how all of this works. Thanks a lot for the detailed description and comments, those make this much easier.
One thing I'm wondering about is the failing unit test. The distance between the expected / the actual point seems to be a bit too large to be explainedjust by rounding errors. OTOH,#40 fixes the issue, which may be either a coincidence or the actual fix. I'll try to investigate a bit in that direction.

There are a few points to consider wrt#40 and precision issues. Most important is that the size of the discrepancy caused by rounding errors, etc., can be a much smaller magnitude than the discrepancy of the distance of the returned neighbor. A discrepancy near the epsilon level in the wrong direction will causedistance_2_if_less_or_equal to returnNone. This can cause the correct neighbor to be disregarded, and other neighbors to be returned instead, which may be arbitrarily far from the correct neighbor. The same is true for pruning in this PR.

I'm running into this right now after making several fixes and refactors to this PR, in that the all nearest neighbors iterator and the single nearest neighbor function return different results rarely. This isn't because one is wrong, but because both are: sometimes one yields the correct neighbor, sometimes the other. This only shows up consistently when exercised with large (> 100K node) 3D trees, but of course may be happening undetected with one or both methods more frequently.

As a real example, here's a node being ignored by the existing single nearest neighbor lookup viadistance_2_if_less_or_equal due to a small difference (0.00000000000000109):

[rstar/src/object.rs:226] distance_2 = 0.008569106766881232[rstar/src/object.rs:226] max_distance_2 = 0.00856910676688123

which results in this being ignored and thus a much larger difference (4%!) in the returned neighbor vs the correct one given byall_nearest_neighbors:

[rstar/src/algorithm/nearest_neighbor.rs:538] single_neighbor.unwrap().distance_2(neighbors.query) = 0.008943686353908992[rstar/src/algorithm/nearest_neighbor.rs:538] neighbors.target.distance_2(neighbors.query) = 0.008569106766881232thread 'algorithm::nearest_neighbor::test::test_all_nearest_neighbors' panicked at 'assertion failed: `(left == right)`  left: `Some([0.016142428043963264, 0.3308523413774229, 0.08940400251818637])`, right: `Some([0.12212996229489037, 0.18978760533844952, 0.07978876040511107])`', rstar/src/algorithm/nearest_neighbor.rs:540:13

Whilemin_max_dist_2 could be further changed to something even more identical to the results oflength_2 (with#40 it differs fromlength_2 only in the order the final additions occur, but this is already enough for small discrepancy due to non-associativity), the underlying issue remains that distances generated different ways (such as throughEnvelope::distance_2 vsPointExt::distance_2) are compared. While returning differentneighbors that differ only by an epsilon is not a problem, comparisons elsewhere will be.

One solution could be adding anepsilon associated const or method toRTreeNum, whichdistance_2_if_less_or_equal and other hard thresholding comparisons can use. For int types this would beZero::zero and for floats it could be their associated::EPSILON. How does that sound?

For now a small change to#40 makes it identical to itslength_2 results before#35 while retaining most of the performance there, so I will push that amended change to#40 later this weekend. That prevents most of the discrepancy issues between different nearest neighbor methods.