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.

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.

share|improve this question
add comment

1 Answer 1

up vote 11 down vote accepted

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

share|improve this answer
    
Terrific, any idea what Nearest is doing? –  RM1618 Sep 24 '13 at 3:18
    
@RM1618 No, I have no idea about how it is implemented internally. –  belisarius Sep 24 '13 at 3:21
1  
Interesting. The r here is undocumented Nearest[data,x,{n,r}]. –  Murta Sep 24 '13 at 3:22
    
@Murta I didn't see it either when I looked through the documentation. –  RM1618 Sep 24 '13 at 3:28
3  
show 6 more comments

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.