Quickly selecting for points in a set that are within a critical distance of points in another set

Question

Imagine two lists of two-dimensional coordinates:

listA = RandomReal[{0,100},{202,2}];
listB = RandomReal[{0,100},{97,2}];

I'm attempting to quickly generate a new series of lists, outputListA and outputListB consisting of the indices for the set of points in listA and listB, respectively, that are within some Euclidean distance $D$ of a point in a list for which they are not a member (i.e. points in listA that are at most a distance distCut from at least one point in listB and vice versa).

This isn't the right way to do things (it takes $\approx 88$ milliseconds for sizes of listA and listB shown), but it hopefully illustrates what I'm trying to do:

listA = RandomReal[{0, 100}, {202, 2}];
listB = RandomReal[{0, 100}, {97, 2}];

outputList = {};
distCut = 1;

For[x = 1, x <= Length[listA], x++,
 For[y = 1, y <= Length[listB], y++,
   If[EuclideanDistance[listA[[x]], listB[[y]]] <= distCut,
     outputList = Append[outputList, {x, y}];
     ];
   ];
 ];

outputListA = Intersection[outputList[[All, 1]], outputList[[All, 1]]];
outputListB = Intersection[outputList[[All, 2]], outputList[[All, 2]]];

Length[outputListA]
Length[outputListB]

A smarter way to proceed might be to round values in listA and listB to a multiple of distCut, and then check for values in the rounded lists that are equal. However, I can't think of a good way to do this that avoids unnecessary attrition / misses points.

A thought - Can we repeatedly apply Nearest for each point $p_i$ in listA until we find a point that is more than a distance distCut from $p_i$? Is Nearest doing anything more sophisticated than sequentially scanning through all of the points in a comparison list and checking Euclidean distances?

Update - I've specified now that outputListA and outputListB should consist of the indices of the points in listA and listB satisfying the distance cutoff criterion.

Dr. belisarius · Accepted Answer · 2013-09-24T05:52:46.407

11

Nearest[]supports the following syntax:

Nearest[data,x,{n,r}] give up to the n nearest elements to x within a radius r

So:

listA = RandomReal[{0, 100}, {202, 2}];
listB = RandomReal[{0, 100}, {97, 2}];
distCut = 1;
f = Nearest[listB];
Timing[s = Union[f[#, {Infinity, distCut}] & /@ listA]]

{0., {{}, {{1.80442, 77.5901}}, {{5.51155, 53.6578}}, {{9.51508, 16.3514}}, {{31.804, 24.5855}}, {{35.4298, 36.839}}, {{65.3657, 85.783}}, {{76.6778, 92.0471}}, {{84.3255, 77.5874}}}}

Manipulate[
 (s = Flatten[Union[f[#, {Infinity, distCut}] & /@ listA], 1]; 
  Show[ListPlot[{listA, listB}, PlotStyle -> {Red, Blue}], 
       ListPlot[s, PlotStyle -> {PointSize[Medium], RGBColor[0, distCut/10, 0]}]]),
    {distCut, 1, 10}]

enter image description here

edited Sep 24 '13 at 05:52

answered Sep 24 '13 at 03:14

Dr. belisarius

115,881
13
203
453

Terrific, any idea what Nearest is doing? – RM1618 Sep 24 '13 at 03:18
@RM1618 No, I have no idea about how it is implemented internally. – Dr. belisarius Sep 24 '13 at 03:21
1

Interesting. The r here is undocumented Nearest[data,x,{n,r}]. – Murta Sep 24 '13 at 03:22
@Murta I didn't see it either when I looked through the documentation. – RM1618 Sep 24 '13 at 03:28
3

@Murta http://reference.wolfram.com/mathematica/tutorial/UsingNearest.html – Dr. belisarius Sep 24 '13 at 03:28
@belisarius Ah, very good. It's just not in the documentation for my Mathematica v9 install. – RM1618 Sep 24 '13 at 03:29
@belisarius Exactly. Another interesting point is that Nearest creates an index when used as a function. It's not documented too. Time do not growth lineary when points number increase. I didn't know about that. – Murta Sep 24 '13 at 03:33
@Murta Ah, that explains the non-linear scaling with the point set size I was observing. – RM1618 Sep 24 '13 at 03:44
Thanks for reminding me about this. I believe I first saw it used by Daniel Lichtblau here: http://stackoverflow.com/a/7207446/618728 – Mr.Wizard Sep 24 '13 at 07:30
+1. This is certainly going to be useful to me. – RunnyKine Sep 24 '13 at 09:35
@Mr.Wizard I've learned it from that answer, but I forgot that :). What a shame. – Dr. belisarius Sep 24 '13 at 19:13

Quickly selecting for points in a set that are within a critical distance of points in another set

1 Answers1

Linked