# Efficient "radius-search-based averaging" algorithm on a set of 3-D points that lends itself to massive parallelization without locking

I'm looking for a good algorithm for "simple radius-search-based averaging" (or "radius search based sub-sampling") as referred to section 3.4 of the DynamicFusion article.

Relevant excerpt from the paper (to save you some reading time):

The set of all unsupported vertices is then spatially sub-sampled using a simple radius search averaging to reduce the vertices to a set of new node positions ... that are at least $$\epsilon$$ distance apart.

As a complication, I need this to run efficiently on a GPU, hence in a setting with massive parallelism without locking.

Here is the best I could conjure up so far (I apologize in advance if I'm straying from local norms for writing pseudocode using MathJax + Markdown) :

1. Store vertices (set $$V$$) in a KD-Tree structure.

2. For each vertex, find all neighbors within fixed radius $$\epsilon$$ and store their indices in a 2D table, $$N$$, where $$N_{i}$$ is the set of fixed-radius neighbors of vertex $$v_{i}$$

3. Initialize Array $$F := \emptyset$$, for filtered vertex indices, array $$A := \emptyset$$, for averages, two atomic counters $$c_{ready} := 0$$ and $$c_{candidate} := 0$$, and the boolean mask array $$M$$ such that $$|M|$$ = $$|V|$$, with each entry of $$M$$ initialized to $$true$$.

4. For each vertex index $$i$$ (in parallel):

1. Retrieve $$v_{i}$$ from $$V$$, $$v_{i} := V(i)$$

2. $$i_{filtered} = 0$$

3. $$c_{filtered} = 0$$ // regular, local, non-atomic variable

4. do

{

1. Loop $$c_{filtered} := c_{candidate}$$ until $$c_{filtered} = c_{ready}$$

2. while ( $$i_{filtered} < c_{filtered}$$)

exit procedure if $$M(i)=false$$ or $$distance(V(F(i_{filtered}), v_{i}) < \epsilon$$, otherwise $$i_{filtered} := i_{filtered} + 1$$

} while (atomic_compare_exchange($$c_{candidate}$$, $$c_{filtered}$$, $$c_{filtered } + 1$$)

5. $$F(c_{filtered}) := i$$

6. atomic_increment($$c_{ready}$$)

7. $$N_{i} = N(i)$$

8. Vertex $$A_i := V_i$$

9. $$c_{vertices} = 0$$ // local counter for averaging

10. For all $$j$$ in \$N_{i}:

1. $$M(j) := false$$
2. $$A_i = A_i + V_j$$
11. $$A_i = A_i / c_{vertices}$$

12. $$A(i) = A_i$$

In the end, we can just mask $$A$$ using $$M$$ to get the final set of averaged vertices. Aside from utter ugliness, the obvious problem with this algorithm is that it doesn't actually produce points that are guaranteed to be $$\epsilon$$ apart, since the averages of these vertex groups are different from the central point of each group. My hope is there will be negligible and there won't be much need to run this recursively (some preliminary experiments seem to confirm this for my input domain).

Another problem is that I'm not sure the whole schema with the counters is going to work properly (my current testing version is single-threaded, python, dumbed-down version of this).

If you see a problem with this approach stemming from parallelization, please save me from hours and hours of debugging :) Also, if you can recommend any viable alternative (or an approach that is guaranteed to produce averages $$\epsilon$$ apart), please do!

Why the original proposal doesn't work

I implemented my original algorithm, and I'll say right away that although it seems to work fine on the CPU, since the two atomics essentially constitute a critical section, on the GPU it causes deadlock. There do exist workarounds, but they critically impact thread occupancy, which makes the idea not worthwhile in my humble understanding.

A better solution

I found a much more suitable algorithm for this that is much more easy to parallelize on GPU.

The idea is simple. Here it is for a 2D set of points, but it's easily extendable to the third spatial dimension.

1. We're aiming to bin all points in a square grid, with each square bin having size $$r$$. This precludes that a rectangular bounding box of the point cloud is pre-computed, which is easy enough to do with a kernel that keeps track of six atomic variables -- maximum and minimum coordinates of the bounding box -- as it's traversing the points. With the minima and maxima computed, it's simple enough to come up with a bin (integral type) coordinate space. In the most basic/simple version of this algorithm, bins can be represented densely by a linear raster array in memory.
2. Subsequent point-binning can be done in a kernel that runs on every point and determines the point's bin based on its coordinates, the pre-computed point cloud extrema / bin coordinate space, and the value of $$r$$. In the basic version, the bin index can be computed directly from its integral coordinates in the bin space. Instead of storing references to points, each bin can simply contain the (atomic) aggregate and the total (atomic) count of all points that happen to reside in the bin, for easy averaging later.
3. Using the aggregates and the counts, we can then average the points within each bin (in a kernel that runs on all bins with non-zero point counts).
4. Finally, we perform merging of each bin with its immediate neighbor(s) in each direction via weighted averaging (using the counts in each bin). It may be more efficient to run this as a separable kernel that loops the entirety of a spatial dimension (right, down, down-right diagonal), but, for my purposes, I think a more suitable variant will have each kernel work on bin pairs, the pair being neighbors in one dimension, thus running twice for each dimension of interest.

If you just do steps (1-3), you have yourself a regular grid downsampling algorithm.

Caveats

While the algorithm certainly satisfies the requirement that no resulting point is within $$r$$ of any other point, it does not make any guarantee that the set of average points has maximal-possible density. In other words, there certainly might be a more optimal way to pick clusters of points and then average each cluster, still satisfying the $$r$$-size-gap requirement, but coming up with many more points in the end.

Potential Improvements

If the points are clustered in such a way that there are large gaps between the clusters, direct raster approach to storing bins can be sub-optimal, since many of the bins will remain empty and will result in poor thread occupancy. A simple way around this is to use a spatial hashing technique (for an example, refer to this article) or a quadtree/octree structure.