Nearest Points in a Point Cloud

I have a large amount of points, typically well over 10^6, in a 3D point cloud - usually quite scattered. I want to carry out various queries on this cloud, the two main ones are: find the nearest points in the cloud to a bunch of other query points (typically 100-1,000 other scattered 3D points) or find all points within a radius of some point.

A dummy example of doing this naively would be:

%set of 10^6 points in my 'point cloud'
N = 1e6;
x = rand(N ,1); y = rand(N ,1); z = rand(N ,1);
%xq,yq,zq, set of 10^3 'query' points
Nq = 1e3;
xq = rand(Nq,1); yq = rand(Nq,1); zq = rand(Nq,1);
d2 = (x-xq').^2 + (y-yq').^2 + (z-zq').^2;
[mind2,inear] = min(d2);
%set of points in the point cloud nearest to the query points
xc = x(inear); yc = y(inear); zc = z(inear); 
dist = sqrt(mind2); %distance between nearest cloud point and query point

%find points in the cloud within some radius of a single point
xp = rand(1); yp = rand(1); zp = rand(1);
radius = 0.5;
d2 = (x-xp).^2 + (y-yp).^2 + (z-zp).^2;
near = d2<=radius^2;
xnear = x(near); ynear = y(near); znear = z(near);

Obviously these methods work but really start to fail for large point clouds (ideally I'd like to be heading to 10^7 or larger) and large sets of query points - I have to create very large matrices for d2 or put it in a slow 'for' loop.

I've seen lots of examples of how to create Octrees to box up points but very few examples of how to then use the Octree for these applications but it would be great to see any examples of these applications if they exist - or any other methods that will help.

[Matt J]

I have to create very large matrices for d2 or put it in a slow 'for' loop.

Doesn't seem like it. You would only need a 10 iteration loop to go from 10^6 to 10^7.

[Nathan]

But 10^6 points and 10^3 queries already takes quite a while and lots of memory and things like Octree are supposed to scale very well with large numbers of points - also yes 'slow' is subjective but if I can get this to run in less than 1s then it'd be very useful.

Answer 1

[Matt J]

See pdist2, knnsearch, and rangesearch.

[Nathan]

Thanks, knnsearch is a partial solution as it does the first without creating a N x Nq matrix - to ammend my above example:

[inear, dist] = knnsearch([x,y,z],[xq,yq,zq]);

I think there are still more scalable methods for large data sets. Creating an Octree should make searching more like log(N) - I just can't find a nice worked example

[John D'Errico]

But MATT did suggest rangesearch, which does exactly what you asked, find all points within a specified distance. knnsearch will also be useful, since you were not clear about exactly which problems you need to solve.