Efficient way to generate random points with a predefined lower bound on their pairwise Euclidean distance

Question

Using Mathematica what is an efficient way to generate a list of $n$ random two dimensional points $\{x_i,y_i\}$ where $i=1,...,n$ so that no two points $p_1$ and $p_2$ in the list has an Euclidean distance lower than say $d$ meaning $\|p_1-p_2\|\leq d$. I came up with the following solution. Though it works I wanted to know if some better method exists. Code

NodeGenetrator[LowerBound_, UpperBound_, DistanceBound_,SampleLength_] := 
Block[{list},
      list = RandomReal[{LowerBound, UpperBound}, {1, 2}];
      For[i = 0, Length@list <= SampleLength - 1, i++,
          list = Module[{NewVal, dist},
                        NewVal = RandomReal[{LowerBound, UpperBound}, 2];
                        (*NewVal=RandomVariate[NormalDistribution[
                        Mean@{LowerBound,UpperBound},DistanceBound],2];*)
                        dist = Map[EuclideanDistance[NewVal, #] &, list];
                        If[Min[dist] >= DistanceBound, 
                           AppendTo[list, NewVal],
                            list
                          ]
                       ];
         ];
      list
     ];
(* Define function parameters *)
LowerBound = 0;
UpperBound = 100;
DistanceBound = 5;
SampleLength = 60;
sample = NodeGenetrator[LowerBound, UpperBound, DistanceBound,SampleLength];
(ListPlot[#, Frame -> True, Axes -> False,PlotStyle -> PointSize[Large],
AspectRatio -> 1] &)@sample

Output

1. As one can see this code uses For loop and also keeps on iterating until number of points requested is not met. This somehow makes the execution time for this function unpredictable!
2. As expected the parameter DistanceBound has a important effect on the function behavior. If we try DistanceBound=12.5 function evaluation becomes very time consuming.
3. Here we check only one mutual distance condition but is it possible to use a more general test function that checks for more than one mutual characteristics involving those requested number of points that are to be generated.

BR

Answer 1

This is not an efficient answer but it is fun to play with so I thought I'd post it. For efficiency the use of Nearest might provide a good starting point.

g[n_, {low_, high_}, minDist_, step_: 1] := 
 Block[{data = RandomReal[{low, high}, {n, 2}], temp, happy, sdata, 
   hdata},

  While[True,
   temp = ((Nearest[data][#, 2][[-1]] & /@ data));
   happy = 
    Boole@Thread[MapThread[EuclideanDistance, {temp, data}] > minDist];
   hdata = Pick[data, happy, 1];
   sdata = Pick[data, happy, 0];
   If[sdata === {}, Break[]];
   sdata += RandomReal[{-step, step}, {Length[sdata], 2}];
   data = Join[sdata, hdata];
   ];
  data
  ]

The idea is to do an initial sampling and then allow the points that are too close to "walk" somewhere else. The function takes a desired number of points n, a low and high value for the data range, the minimum acceptable distance between points minDist and a step argument which allows points to "walk" up to a certain distance in the x and y directions.

Its especially fun to watch dynamically.

g[150, {0, 30}, 1.5, 1]

Edit: Per suggestion of Yves Klett the points are colored by relative happiness (green being more happy, red being less happy).

Edit 2:

Now for a more serious attempt at something efficient..

findPoints =
  Compile[{{n, _Integer}, {low, _Real}, {high, _Real}, {minD, _Real}},
   Block[{data = RandomReal[{low, high}, {1, 2}], k = 1, rv, temp},
    While[k < n,
     rv = RandomReal[{low, high}, 2];
     temp = Transpose[Transpose[data] - rv];
     If[Min[Sqrt[(#.#)] & /@ temp] > minD,
      data = Join[data, {rv}];
      k++;
      ];
     ];
    data
    ]
   ];

And to test it for benchmarking...

npts = 350;
minD = 3;
low = 0;
high = 100;

AbsoluteTiming[pts = findPoints[npts, low, high, minD];]

==> {0.0312004, Null}

Check that the min distance is less than the threshold.

Min[
  MapThread[
   EuclideanDistance, {pts, Nearest[pts][#, 2][[-1]] & /@ pts}]] > minD

==> True

Check that I generated the correct number of points..

Length[pts] == npts

==> True

Answer 2

If you don't need high precision, you can do something along these lines:

canvas = Image@ConstantArray[0, {100, 100}];

distance = 6;

{img, {pts}} = 
  Reap[Nest[
    ImageCompose[#, 
      SetAlphaChannel[#, #] &@Image@DiskMatrix[distance], 
      Sow@RandomChoice@
        Position[Transpose@ImageData[#, DataReversed -> True], 0.]] &,
     canvas, 150]];

This is not fast in general, but it does not slow down when the density is very high. Note how the region is nearly completely filled up.

{Framed[img], Graphics@Point[pts]}

The idea can be used to speed up other methods when the point density becomes high.

Answer 3

Well, a very simple but very limited way could be to generate a random set of points, calculate all distances between them and scale the minimum distance to mindist:

mindist = 1;
npts = 200;
pts = RandomReal[{0, 100}, {npts, 2}];
scaledpts = mindist/Min[Norm /@ Subtract @@@ Subsets[pts, {2}]]*pts;

Graphics[{Green, Point[pts], Red, Point[scaledpts]}, Frame -> True]

Edit: On the plus side this should work for arbitrary distributions

mindist = .01;
npts = 1000;
pts = RandomVariate[NormalDistribution[1, 3], {1000, 2}];
scaledpts = mindist/Min[Norm /@ Subtract @@@ Subsets[pts, {2}]]*pts;

Graphics[{Green, Point[pts], Red, Point[scaledpts]}, Frame -> True]

Of course, this becomes ugly fast (and hopefully no point is duplicate!). But for small sample numbers, why not? Also, a nice way to see how much memory you´ve got (my machine with 16GB maxed out between 9500 and 10000 points).

See the humorous memory graph with successive runs (1000 pts increase). Rather interesting sawtooth there which releases the memory nicely each time.

Answer 4

Not very pretty, but you could do something like this

generate[pts_, bnds_, mindist_] := Module[{x, plist, ylist, intervals},
  x = RandomReal[bnds[[1]]];
  plist = Pick[pts, UnitStep[Abs[x - pts[[All, 1]]] - mindist], 0];
  If[Length[plist] == 0, 
   Return[Join[pts, {{x, RandomReal[bnds[[-1]]]}}]]];
  ylist = 
   List @@ IntervalIntersection[Interval[bnds[[1]]], 
     IntervalUnion @@ (Interval[#[[2]] + 
           Sqrt[mindist^2 - (#[[1]] - x)^2] {-1, 1}] & /@ plist)];
  intervals = If[Length[ylist] >= 2,
    Transpose[{ylist[[;; -2, 2]], ylist[[2 ;;, 1]]}], {}];
  If[ylist[[1, 1]] > bnds[[1, 1]], 
   PrependTo[intervals, {bnds[[1, 1]], ylist[[1, 1]]}]];
  If[ylist[[-1, 2]] < bnds[[1, -1]], 
   AppendTo[intervals, {ylist[[-1, 2]], bnds[[1, -1]]}]];
  If[intervals =!= {},
   Join[pts, {{x, 
      RandomReal@
       RandomChoice[(#2 - #1) & @@@ intervals -> intervals]}}],
   pts]
  ]

points[bnds_, mindist_, n_] :=
  NestWhile[generate[#, bnds, mindist] &, {}, Length[#] < n &, 1, 2 n]

Here, bnds is the domain of the random points in the form {{xmin, xmax}, {ymin, ymax}}, mindist is the minimum distance between two points, and n is the number of generated points.

Example

bnds = {{0, 100}, {0, 100}};
mindist = 5;
npts = 150;
pts = points[bnds, mindist, npts];

ListPlot[pts, Frame->True]

Explanation

The idea is to iteratively generate new points by choosing an arbitrary value $x_0$ in the interval $[x_\text{min}, x_{\text{max}}]$. For each already generated point p that is closer than mindist to the line $x=x_0$ we find the interval $[y_1, y_2]$ such that $||p-(x_0,y)||\leq \text{mindist}$ for $y_1\leq y\leq y_2$ and choose an arbitrary value for $y$ from the complement of the union of all these intervals (provided this set is non-empty).

If the plane starts filling up, the chance of choosing a value for x for which there is no available y value becomes larger, so we still need to generate extra points. To reduce the number of extra generated points, this method could be combined with Archimedes' solution, for example by periodically checking for which intervals for x there are still y values available, and choosing only from those.

Answer 5

Here is an outline of a fairly efficient way to do this.

Say you want to generate $m$ points with a minimum separation of $d$, in a box of size $s \times s$. Let $d_2=\max(d,\sqrt{m}/s)$. Partition your box into subsquares of size $d_2 \times d_2$; these will be used as bins. We will generate, say, $3m$ random points in the full box. Now place each in its appropriate bin. This can be done in $\mathcal{O}(m)$.

Now iterate over the generated points. For each, locate its bin, and gather all points in neighboring subsquare bins. For each such neighbor, see if it is within distance $d$ of the candidate point and has already been chosen. If not, the candidate point is chosen as a member of our random set, and a down value is set to indicate this (so checking each neighbor is $\mathcal{O}(1)$).

We stop once we have $m$ points collected, or else run out of points. If the latter, generate more, add them to the existing binds, and proceed as before.

The memory hoof-print is $\mathcal{O}(d_2^2)$. The speed will typically be $\mathcal{O}(m)$ unless the sizes of $d$ and $s$ make it either difficult or impossible to have $m$ points inside with minimal separation of distance $d$.

Answer 6

A more recent question got pointed here as it being a duplicate of this one, so I am adding one of my answers from there to here.

None of the answers here so far make efficient use of Nearest[], which is quite powerful for a large number of points as it makes use of a partitioning algorithm. Nearest[] can generate a NearestFunction[] for a set of points, which does that partitioning once. Then the generated function can be used repeatedly over many points.

Here cull[] takes a list of points and removes enough of those points to meet the condition that no two points are closer than d from each other. Then many candidate points can be added and a single NearestFunction[] used to cull the points. Note that this will work for any number of dimensions, as well as measures other than Euclidean.

cull[pts_, d_] :=
 Module[{p, f, c},
  p = pts;
  While[
   f = Nearest[p -> Automatic];
   c = f[p, {2, d}];
   Last[Dimensions[c]] != 1,
   p = Pick[p, First[#] <= Last[#] & /@ c];
  ];
  p
 ]

The use of the efficient algorithm becomes important when generating a large number of points, and especially when trying to fit as many randomly selected points as possible into a region while meeting the condition.

spread[] uses cull[] as described, adding n dim-dimensional points to the list each time, and culling, repeating until n points remain. n and d can be chosen to be too large relative to high-low, making it impossible to fit that many points in the region. For that reason, spread[] also has a maximum number of iterations for adding the last point. Without a maximum, the algorithm can provide no assurance that it will terminate. Each iteration tries n samples. When n*maxiter points are tried with zero points successfully added, spread[] returns the points it has, which meet the bounds and distance condition, but the count is short of n.

spread[n_, dim_, d_, low_, high_, maxiter_] :=
 Module[{p, m, k},
  p = {};
  m = 0;
  k = maxiter;
  While[
   p = Join[p, RandomReal[{low, high}, {n, dim}]];
   p = cull[p, d];
   k = If[Length[p] > m, maxiter, k - 1];
   m = Length[p];
   m < n && k > 0
  ];
  If[Length[p] > n, Take[p, n], p]
 ]

This example use of spread[] casts the result in the form of circles, which makes it easy to visualize the result and see that the points have the required separation, since the circles don't overlap.

circles[n_, r_, maxiter_] := 
 Circle[#, r] & /@ spread[n, 2, 2 r, -1, 1, maxiter]

Then:

circles[100, 0.1, 1000] // Graphics

gives 80 circles (100 is too many random 0.1 radius circles in a 2x2 box):

Here is an example with spheres:

spheres[n_, r_, maxiter_] := 
 Sphere[#, r] & /@ spread[n, 3, 2 r, -1, 1, maxiter]

spheres[150, 0.2, 1000] // Graphics3D

Returning 127 radius 0.2 spheres, since 150 is too many in a 2x2x2 box:

Answer 7

A partial answer based on low-discrepancy Halton sequences (see the Quasi- Random Number Generators section of the tutorial CUDALink/tutorial/Applications in the docs):

Using

 VanDerCorput[base_][len_] := Table[
 With[{digits = Reverse@IntegerDigits[n, base]},
  Sum[2^(-ii)*digits[[ii]], {ii, Length[digits]}]], {n, len}]  

as the generator, one can generate 2D Halton sequences using two different bases and combining the resulting numbers in pairs:

 sample = Transpose[{VanDerCorput[2][1000], .5 VanDerCorput[3][1000]}]

where I rescaled the second component to constrain the points to unit square. Halton sequences fill the unit square uniformly as seen in the plot below:

Halton sequences are good candidates for generating the raw samples on which further selections using many cool ideas put forth in other answers.

For instance, using Yves Klett's scaling idea:

 mindist = .02;
 scaledsample =  mindist/Min[Norm /@ Subtract @@@ Subsets[sample, {2}]]*sample;

we get