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?

Improve this question
$\endgroup$
9
  • $\begingroup$ try using Nearest $\endgroup$
    – george2079
    Sep 24, 2015 at 16:31
  • 3
    $\begingroup$ pts2 = Union[pts, SameTest -> (Norm[#1 - #2] < 0.05 &)]; $\endgroup$
    – Bob Hanlon
    Sep 24, 2015 at 16:35
  • 1
    $\begingroup$ Possible duplicate: (32409) Related: (2594) $\endgroup$
    – Mr.Wizard
    Sep 24, 2015 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, 2015 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, 2015 at 1:13

5 Answers 5

Reset to default
17
$\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

Improve this answer
$\endgroup$
10
  • 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, 2015 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, 2015 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, 2015 at 22:30
  • $\begingroup$ @march Thank you, that works. $\endgroup$
    – Jack LaVigne
    Sep 25, 2015 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, 2015 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

Improve this answer
$\endgroup$
0
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

Improve this answer
$\endgroup$
1
  • $\begingroup$ Your grid would look a bit uniformer if you applied aspect ratio of 1 $\endgroup$
    – Alexey Bobrick
    Oct 6, 2015 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..

Improve this answer
$\endgroup$
6
  • $\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, 2015 at 4:01
  • $\begingroup$ @ciao Might be of interest for other users though... $\endgroup$
    – anderstood
    Sep 25, 2015 at 5:07
  • $\begingroup$ Doesn't run here :( !Mathematica graphics - Have you copied faithfully? $\endgroup$
    – Dr. belisarius
    Sep 25, 2015 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, 2015 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, 2015 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]

Improve this answer
$\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.