find all points within a range to any point of an other set

Question

I have two sets of points A and B.

I want to find all points in B that are within a certain range r to A, where a point b in B is said to be within range r to A if there is at least one point a in A whose (Euclidean) distance to b is equal or smaller to r.

Each of the both sets of points is a coherent set of points. They are generated from the voxel locations of two non overlapping objects.

In 1D this problem fairly easy: all points of B within [min(A)-r max(A)+r]

But I am in 3D.

What is the best way to do this?

I currently repetitively search for every point in A all points in B that within range using some knn algorithm (ie. matlab's rangesearch) and then unite all those sets. But I got a feeling that there should be a better way to do this. I'd prefer a high level/vectorized solution in matlab, but pseudo code is fine too :)

I also thought of writing all the points to images and using image dilation on object A with a radius of r. But that sounds like quite an overhead.

Answer 1

You can use a k-d tree to store all points of A.

Iterate points b of B, and for each point - find the nearest point in A (let it be a) in the k-d tree. The point b should be included in the result if and only if the distance d(a,b) is smaller then r.

Complexity will be O(|B| * log(|A|) + |A|*log(|A|))

Answer 2

I archived further speedup by enhancing @amit's solution by first filtering out points of B that are definitely too far away from all points in A, because they are too far away even in a single dimension (kinda following the 1D solution mentioned in the question).

Doing so limits the complexity to O(|B|+min(|B|,(2r/res)^3) * log(|A|) + |A|*log(|A|)) where res is the minimum distance between two points and thus reduces run time in the test case to 5s (from 10s, and even more in other cases).

example code in matlab:

r=5;
A=randn(10,3);
B=randn(200,3)+5;

roughframe=[min(A,[],1)-r;max(A,[],1)+r];

sortedout=any(bsxfun(@lt,B,roughframe(1,:)),2)|any(bsxfun(@gt,B,roughframe(2,:)),2);
B=B(~sortedout,:);
[~,dist]=knnsearch(A,B);
B=B(dist<=r,:);

Answer 3

bsxfun() is your friend here. So, say you have 10 points in set A and 3 points in set B. You want to have them arrange so that the singleton dimension is at the row / columns. I will randomly generate them for demonstration

A = rand(10, 1, 3);                    % 10 points in x, y, z, singleton in rows
B = rand(1, 3, 3);                     %  3 points in x, y, z, singleton in cols

Then, distances among all the points can be calculated in two steps

dd = bsxfun(@(x,y) (x - y).^2, A, B);  % differences of x, y, z in squares
d = sqrt(sum(dd, 3));                  % this completes sqrt(dx^2 + dy^2 + dz^2)

Now, you have an array of the distance among points in A and B. So, for exampl, the distance between point 3 in A and point 2 in B should be in d(3, 2). Hope this helps.