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

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

-
Perhaps I am getting this backwards, but why not generate random points within a given envelope (say, a circle) and scale the random coordinates with the largest possible distance within that envelope (e.g. the diameter) to your parameter d? Or are p1 and p2 sequential elements of the list? –  Yves Klett Mar 5 '12 at 17:12
You probably meant .. so that NO two points ... has a Euclidean distance lower than say d in the second line? –  kguler Mar 5 '12 at 17:18
"Random" does not mean merely arbitrary; nor, in practice, does it mean (merely) that random numbers were involved. To use "random" points, you need to know their probability distribution. (There exist efficient solutions to this problem having very different distributions.) So: precisely what distribution do you want these points to have? –  whuber Mar 5 '12 at 18:22
I would like to suggest you do have a requirement on the distribution. Otherwise, find two children and ask them to draw pictures of solutions. Check the solutions and digitize them. Ask Mathematica to choose randomly between the two solutions. After the initial computation of the two solutions, this is extremely efficient and is obviously random. I'm not being facetious: this is a nice example of a spatial stochastic process that meets every one of your stated requirements. It might help you see why specifying a distribution is important. –  whuber Mar 5 '12 at 20:11
For practical solutions which may meet your needs, start with the Wikipedia article on low-discrepancy sequences and then look closely at the Halton sequence. –  whuber Mar 5 '12 at 21:02

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 =
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[
EuclideanDistance, {pts, Nearest[pts][#, 2][[-1]] & /@ pts}]] > minD

==> True


Check that I generated the correct number of points..

Length[pts] == npts

==> True

-
You should color the points according to their degree of happiness (or claustrophobia)! –  Yves Klett Mar 5 '12 at 19:38
@Yves That would be yet more fun.. at least when animated we see this to some extent based on their relative laziness :) –  Andy Ross Mar 5 '12 at 19:43
Improvement: instead of calculating a square root for each point, calculate the square of the distance. (Also you could iterate over the points and break at the first one that is found to be too close, instead of calculating all their distances and checking the minimum, but I'm not sure if that'd be faster). –  Gerardo Marset Feb 14 '13 at 23:52
@GerardoMarset that's a good suggestion. Give it a try. If the timings are better feel free to edit. –  Andy Ross Feb 15 '13 at 4:13
@AndyRoss I'm sorry, can't. I don't have Mathematica. –  Gerardo Marset Feb 15 '13 at 5:09

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.

-
This is a very creative answer, +1! There seems to be a relationship between the dimensions of canvas and the maximum number of possible iterations in the Nest function (and therefore number of points). Changing the 150 to 180 without also increasing canvas gives a range of errors including RandomChoice::lrwl: "The items for choice {} should be a list or a rule weights -> choices." –  Verbeia Mar 5 '12 at 21:03
@Verbeia That is correct, when there's no more space for new points, RandomChoice will return an empty list. I just illustrated the context without a robust implementation. Ideally this method could be used to weed out positions where a new point definitely can't be placed, on a square grid. Then choose a cell in the grid, and randomly position a new point inside. Then do another test to see if the point is too close to any points in nearby cells. This will give a method which is both fast and precise, but it's some work to implement it right. –  Archimedes of Syracuse Mar 5 '12 at 21:22
@ArchimedesofSyracuse This is what I am doing with this ind of disjoint squares last one hour. Not an easy implementation for me. On the other hand I also think we can try to compute the probability distribution explicitly by hand. Then use this distribution in MMA to find the points. But may be it is more difficult mentally than efficient code writing.. –  PlatoManiac Mar 5 '12 at 21:31
@PlatoManiac The big problem with calculating the probability distribution and using that is that the points positions are not independent, and can't be independently generated ... so we're back to square one. –  Szabolcs Apr 9 '13 at 14:32

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.

-

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.

-
This scales impressively well +1. –  Andy Ross Mar 6 '12 at 3:31

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$.

-
There, you had to start with the Big O notation. Now everyone will start benchmarking - think of all the entropy ;-) –  Yves Klett Mar 6 '12 at 17:31
No need to benchmark, @Yves: I have already claimed the prize for fastest (c. 500,000,000 point patterns per second). See comments to the question. :-) –  whuber Mar 6 '12 at 18:44

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

-
Is it just my eyes, or is that not uniform? It appears heavy in the middle, and sparse at the top and bottom. –  Mr.Wizard Mar 6 '12 at 10:07
@Mr.Wizard, your eyes are ok. Need to re-check the code and the data. –  kguler Mar 6 '12 at 10:37
@kguler, did you maybe find what was wrong with uniformity? –  VividD Jun 20 '14 at 9:02

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 = sweep[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:

-