Finding the nearest point

Question

Algorithms

Algorithmic abstract problem:

Given: 1,000,000 black dots and 1,000,000 red dots (x,y,z: float);
— 10,000 < x,y,z < 10,000 (Cartesian coordinate system).

You need to find the sum of the distances from each red dot to the nearest black dot from it.
Example:


I seriously considered 2 options:
KD - trees - but there is a sad counterexample for them, when all black points are located at the same distance from all red ones (red - a segment; black - a ring around any of the points of the red segment);

Voronoi diagram- but it is very sad for any situation, difficult to program (especially for the 3D version, it is problematic to build a tree from a diagram in order to have a logarithmic descent). This algorithm, despite the complexity of finding the nearest point o(log(n)) (after constructing o(???)), will be slow. Determining whether a point belongs to a polyhedron is a complex operation.

Maybe there are alternative algorithms? I know that there are a lot of varieties of this problem, but this particular case and the most efficient algorithm are of interest.

Answer 1

Thank you for the challenge. She interested me a lot.
I sketched one algorithm on my knee, the implementation in C # works pretty well: ~ 4.5 sec on my decrepit machine. If you correctly parallelize (and this should be quite simple for this algorithm), then it should be even faster.
The only difference from your statement of the problem: I set the coordinates of the points as double, not float.
The essence of the algorithm : go to integers in the sense that cover the ODZ (-10,000 < x, y, z < 10,000) with a three-dimensional grid.
For the case of a uniform distribution of points, the easiest way is to use a uniform grid.
Experiments have shown that the optimal grid spacing for the parameters specified in the condition of the problem lies in the range of 100-400, but this is not important.
It is necessary to decompose all points of the same color (let it be red) into cells. Personally, I decided to lay out not the points themselves, but their indices in the original array. This is easy to do, especially in the uniform case. That is, in C# language, the following three-dimensional array is filled:

List<int>[,,] grid

The elements of this array are lists of indices of red dots that fall into the corresponding cell.
Next, we have to do something that I call a "spherical cycle" (I don't know if there is any standard name for such structures). The essence of such a cycle can be explained on the fingers as follows: we go around the cells of the grid, starting from some initial triple of indices i ₀ , j ₀ , k ₀ , and then diverging into concentric spheres. In other words, first the loop hits the initial cell, then runs through all the cells that are 1 (rounded) away from the initial one, then 2 (rounded), and so on.
Further, in the most ordinary cycle, we run through all the black points, calling for each of them this very “spherical cycle”, in which we are looking for the nearest neighbor. True, there is still a trick with a condition for exiting this cycle, but these are trifles.
There is a source code in C#. If there is interest, I can provide it.

Answer 2

The easiest way is by iterating:

var sqr = function(x) { return x * x; };
var summ = 0;
for (var i in red_points)
{
  var minSLen = false;
  for (var j in black_points)
  {
    var sLen = sqr(black_points[j].x - red_points[i].x) + sqr(black_points[j].y - red_points[i].y) + sqr(black_points[j].z - red_points[i].z);
    if (minSLen === false || sLen < minSLen)
      minSLen = sLen;
  }
  if (minSLen >= 0)
    summ += Math.sqrt(minSLen);
}

Its efficiency is rather low, but it works in any conditions :-)