Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

enter image description here

  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

share|improve this question
    
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
3  
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
1  
"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
2  
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
4  
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

6 Answers 6

up vote 21 down vote accepted

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

enter image description here

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
share|improve this answer
2  
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

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]

Mathematica graphics

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]

Mathematica graphics

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.

Mathematica graphics

share|improve this answer

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]}

Mathematica graphics

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

share|improve this answer
1  
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

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]

Mathematica graphics

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.

share|improve this answer
    
This scales impressively well +1. –  Andy Ross Mar 6 '12 at 3:31

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

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

halton scaled for min distance

share|improve this answer
2  
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 at 9:02

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

share|improve this answer
    
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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