13
$\begingroup$

Consider a list of points:

pts = Partition[RandomReal[1, 10000], 2];
ListPlot[pts]

enter image description here

I'd like to delete points so that the minimum distance between two points is 0.05. The following code does the job:

pts2 = {pts[[1]]};
Table[If[Min[Map[Norm[pts[[i]] - #] &, pts2]] > 0.05, 
AppendTo[pts2, pts[[i]]]], {i, 2, Length[pts], 
1}]; // AbsoluteTiming (* -> 1.35 *)
ListPlot[pts2]

enter image description here

But it becomes slow for large lists, probably because of AppendTo which does not know what type is going to come next.

How could this be done more efficiently? Note: there is no uniqueness of the resulting list, but that's not a problem.

Just for better referencing, let me give another formulation of the question: How to delete points in a neighbourhood of other points of a list?

$\endgroup$
  • $\begingroup$ try using Nearest $\endgroup$ – george2079 Sep 24 '15 at 16:31
  • 3
    $\begingroup$ pts2 = Union[pts, SameTest -> (Norm[#1 - #2] < 0.05 &)]; $\endgroup$ – Bob Hanlon Sep 24 '15 at 16:35
  • 1
    $\begingroup$ Possible duplicate: (32409) Related: (2594) $\endgroup$ – Mr.Wizard Sep 24 '15 at 17:10
  • 2
    $\begingroup$ Just for you to keep in mind: Union with SameTest option set explicitly, has quadratic complexity in the number of points, because it performs pairwise conparisons. $\endgroup$ – Leonid Shifrin Sep 24 '15 at 22:30
  • $\begingroup$ Is there a reason you're generating then deleting, vs just properly generating in the first place? There are algorithms for blue noise that will easily generate millions of points / sec with the conditions of OP... $\endgroup$ – ciao Sep 25 '15 at 1:13
15
$\begingroup$

The following is a much faster, but not optimal, recursive solution:

pts = RandomReal[1, {10000, 2}];
f = Nearest[pts];

k[{}, r_] := r
k[ptsaux_, r_: {}] := Module[{x = RandomChoice[ptsaux]}, 
                      k[Complement[ptsaux, f[x, {Infinity, .05}]],  Append[r, x]]]

ListPlot@k[pts]

Mathematica graphics


Some timings show this is two orders of magnitude faster than the OP's method:

ops[pts_] := Module[{pts2},
  pts2 = {pts[[1]]};
  Table[If[Min[Map[Norm[pts[[i]] - #] &, pts2]] > 0.05, 
    AppendTo[pts2, pts[[i]]]], {i, 2, Length[pts], 1}];
  pts2]

bobs[pts_] := Union[pts, SameTest -> (Norm[#1 - #2] < 0.05 &)]

belis[pts_] := Module[{f, k},
  f = Nearest[pts];
  k[{}, r_] := r;
  k[ptsaux_, r_: {}] := Module[{x = RandomChoice[ptsaux]}, 
                        k[Complement[ptsaux, f[x, {Infinity, .05}]], Append[r, x]]];
  k[pts]]


lens = {1000, 3000, 5000, 10000};
pts = RandomReal[1, {#, 2}] & /@ lens;
ls = First /@ {Timing[ops@#;], Timing[bobs@#;], Timing[belis@#;]} & /@  pts;
ListLogLinePlot[  MapThread[List, {ConstantArray[lens, 3], Transpose@ls}, 2], 
               PlotLegends -> {"OP", "BOB", "BELI"}, Joined ->True]

Mathematica graphics

$\endgroup$
  • 3
    $\begingroup$ @belisarius And therein lies one of my biggest gripes... the "I really need this, but I'm accepting something else..." switcheroo that happens here way too often. So now, a user searches for "fast...", gets this question as a hit, and is led to the slowest solution. DOH! +1 on you, shame on that kind of nonsense. $\endgroup$ – ciao Sep 24 '15 at 19:47
  • $\begingroup$ @belisarius Your code'efficiency impressed me a lot.But if the element is a Disk with diffrence radius instead of Point.Can you do same thing efficiencly?I'm confusion with this problem in long time,Could you help me to try it? $\endgroup$ – yode Sep 24 '15 at 20:58
  • $\begingroup$ @belisarius - Amazing. How do you think of these things? One question I have is I thought that Append was to be avoided., How is it that you are able to get away with this? Second question, I tried Trace but it pretty much upchucked. How can I figure out how many steps are used in the recursion? $\endgroup$ – Jack LaVigne Sep 24 '15 at 22:30
  • $\begingroup$ @march Thank you, that works. $\endgroup$ – Jack LaVigne Sep 25 '15 at 0:26
  • 1
    $\begingroup$ @JackLaVigne Knowing that Nearest[ ] is the most efficient Mathematica function for ranking distances helps to think these things :). WRT your concerns about Append[ ], you are right but .... I'm appending only a few points. Most of the points get discarded. If you want to improve this little program try to optimize the DeleteCases[ ]instead of the Append[ ] part :) $\endgroup$ – Dr. belisarius Sep 25 '15 at 2:21
15
$\begingroup$
pts = Partition[RandomReal[1, 10000], 2];

ListPlot[pts]

enter image description here

Use SameTest option with Union

pts2 = Union[pts, SameTest -> (Norm[#1 - #2] < 0.05 &)];

Length[pts2]

326

ListPlot[pts2]

enter image description here

$\endgroup$
5
$\begingroup$

The following "solution" has the benefits of:

  • making a very a uniform grid.

  • being fast.

It has the (perhaps mortal) drawbacks of:

  • not being automated.

  • being pretty liberal about kicking out points.

Nonetheless, I wanted to play a little. Here's my take: generate a square grid of points and use Nearest to pick out the points nearest to the gridpoints:

pts = Partition[RandomReal[1, 10000], 2];
nearestOnGrid[points_, d_] := Nearest[points, Outer[List, Range[0, 1, d], Range[0, 1, d]]~Flatten~1]~Flatten~1
testDistances[grid_, leastD_] := Min[EuclideanDistance @@@ grid~Subsets~{2}] < leastD

Then, if we do

grid = nearestOnGrid[pts, 0.074]; // AbsoluteTiming
testDistances[grid, 0.05] // AbsoluteTiming
(* {0.000957, Null} *)
(* {0.016401, True} *)

Note that the choice of 0.074 was not automated. I used testDistances to find a value for the grid-spacing that made it True. However, since this takes 0.016 seconds, trying to automate the procedure with some sort of bracketing method will definitely make this slower than the rest of the options above.

Nonetheless, the results are:

GraphicsRow[{ListPlot[pts], ListPlot[grid]}, ImageSize -> 600]

enter image description here

$\endgroup$
  • $\begingroup$ Your grid would look a bit uniformer if you applied aspect ratio of 1 $\endgroup$ – Alexey Bobrick Oct 6 '15 at 15:40
4
$\begingroup$

one more.. this I think fares well for very large n.

result = NestWhile[ 
            Nest[ Complement[#, Rest@Nearest[ # , RandomChoice[#] ,
              { Infinity, .05}]] & , #, Ceiling[(Length@#)/100] ] &, pts,
                  Min[EuclideanDistance @@@ Nearest[#, #, 2]] < .05 & ];

Kind of ugly to double Nest but the convergence test is the expensive part..

$\endgroup$
  • $\begingroup$ No need to do such an expensive test - I cobbled up a similar idea (did not bother posting since OP accepted a slow answer) - just nest ~Sqrt[numpts], then do a pass over remaining points. 3-4X faster than belisarius's neat answer. +1, o/c... $\endgroup$ – ciao Sep 25 '15 at 4:01
  • $\begingroup$ @ciao Might be of interest for other users though... $\endgroup$ – anderstood Sep 25 '15 at 5:07
  • $\begingroup$ Doesn't run here :( !Mathematica graphics - Have you copied faithfully? $\endgroup$ – Dr. belisarius Sep 25 '15 at 5:08
  • $\begingroup$ @belisarius: New form for Nearest... does not work on earlier V of MMA (not sure when added, but it's the Nearest[{...},{...}] form, equivalent (but I'd imagine faster than) Nearest[{...},#]&/@{...} $\endgroup$ – ciao Sep 25 '15 at 5:35
  • $\begingroup$ @ciao Thanks! They will keep overloading those functions syntax until one doesn't know what the heck a program is doing anymore :) $\endgroup$ – Dr. belisarius Sep 25 '15 at 5:45
1
$\begingroup$

Another style of coding by NestWhile instead of recursion in current accepted answer.

pts = RandomReal[1, {10000, 2}];
f = Nearest[pts];

pts = Last[NestWhile[Apply[{Complement[#, f[x = RandomChoice[#], {All, .05}]], 
       Append[#2, x]} &], {pts, {}}, Length@First[#] != 0 &]];
ListPlot[pts]

$\endgroup$

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.