I have a large array of 2d points, e.g.:


nPoints = 100000;

xmin = 0; ymin = 15;
xmax = 45; ymax = 32;

points = Transpose[{RandomReal[{xmin, xmax}, nPoints],
                    RandomReal[{ymin, ymax}, nPoints]}

Now I want to count the number of neighbouring points within a certain radius throughout a regular grid in certain steps (dx and dy).

My simple counting code ist the following:

number = Array[0 &, {Round@((xmax - xmin)/dx), Round@((ymax - ymin)/dy)}];

dx = 0.125; dy = 0.125;
radius = 2.;

     number[[i, j]] = 
        Length@Nearest[points, {xmin + (i - 1)*dx + dx/2, ymin + (j - 1)*dy + dy/2}, 
           {All, radius}];
     , {j, 1, Round@((ymax - ymin)/dy), 1}
   , {i, 1, Round@((xmax - xmin)/dx), 1}
   ]; // AbsoluteTiming

{22.7748, Null}

The corresponding ListDensityPlot of number can be plotted with:

ListDensityPlot[Transpose@number, ColorFunction -> "Rainbow", AspectRatio -> Automatic, 
   PlotRangePadding -> None, Frame -> True, InterpolationOrder -> 0]

enter image description here

The color represents regions of a certain number of points within radius.

Since I have to process the Do loops a few thousand times this solution is simply far too slow.

What is a fastest possible solution to calculate the Array number?


It's already a very good idea to use Nearest. You can improve that a little by defining the NearestFunction once. Moreover, a NearestFunction knows its input dimension so that it can automatically thread over its input. No loop constructs needed. Using also points-> Automatic releaves Nearest from the duty of reading all the coordinates from points. On my machine, this leads to a 15-fold speedup:

nPoints = 100000;
xmin = 0.; ymin = 15.; xmax = 45.; ymax = 32.;
points = Transpose[{RandomReal[{xmin, xmax}, nPoints], 
    RandomReal[{ymin, ymax}, nPoints]}];
dx = .125;
dy = .125;
radius = 2.;

grid = Developer`ToPackedArray@Table[
    {x, y},
    {x, xmin + dx/2, xmax, dx},
    {y, ymin + dy/2, ymax, dy}
number2 = Map[Length, Nearest[points-> Automatic, grid, {All, radius}], {2}]; // AbsoluteTiming

number = ConstantArray[0, Dimensions[grid][[1 ;; 2]]];
Do[Do[number[[i, j]] = 
         points, {xmin + (i - 1)*dx + dx/2, 
          ymin + (j - 1)*dy + dy/2}, {All, radius}];, {j, 1, 
      Round@((ymax - ymin)/dy), 1}];, {i, 1, Round@((xmax - xmin)/dx),
     1}]; // AbsoluteTiming
number == number2

{1.62238, Null}

{24.2462, Null}



Setting nPoints = 1000000 leads to a pretty bad blow up of memory usage (not to mention runtime). Clearly, it is suboptimal to store the list of all points within given radius around each grid point when we only require the length of these lists...

The following is slightly slower but more robust under cranking up nPoints. The idea is to find for each point in points the grid points with distance lower than radius and to add 1 to the corresponding positions in gridcounter. Moreover, its still persisting high memory requirements can easily be adressed by cutting points into smaller chunks and handing over these chunks to the NearestFunction one after another.

number3 = Module[{g, gridcounters, data},
    g = Nearest[Flatten[grid, 1] -> Automatic];
    gridcounters = ConstantArray[0, Times @@ Dimensions[grid]];
    data = g[points, {All, radius}];
    Do[gridcounters[[idxlist]] += 1;, {idxlist, data}];
    g =.;
    data =.;
    ArrayReshape[gridcounters, Dimensions[grid][[1 ;; 2]]]
    ]; // AbsoluteTiming
    number3 == number2

{2.27033, Null}


Another remark

If you are also satisfied with square-shaped balls (balls with respect to the $\infty$-norm) then you could also apply BinCounts in conjunction with a suitable scaling.

  • $\begingroup$ Thank you ... sorry I forgot to set dx, dy, radius - I updated my question with the values that you have chosen. btw you can replace Max[Abs[number - number2]] by number == number2 which gives True. $\endgroup$ – mrz Nov 24 '17 at 14:25
  • $\begingroup$ Yes, with integer matrices, == is a good way. But with real matrices, this leads often to false negatives. So, it's just a habit of mine. $\endgroup$ – Henrik Schumacher Nov 24 '17 at 14:30
  • $\begingroup$ Both of your solutions, especially where you caclulate number3, are very sophisticated. $\endgroup$ – mrz Nov 27 '17 at 14:48
  • 1
    $\begingroup$ Thanks! (I hope that is meant positively...) $\endgroup$ – Henrik Schumacher Nov 27 '17 at 14:55
  • $\begingroup$ Absolutely ... I never thought about the memory usage problem and also your idea to search for each point the corresponding grid point is interesting. And finally the performance of your solution is much better. Once again ... thanks a lot for the code. $\endgroup$ – mrz Nov 27 '17 at 16:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.