Selecting for 2D points that are within a threshold distance of an upper- and lower-bound number of points

Question

I have a very large set of 2D points:

numberOf2DPoints = 10^6;
pointList = RandomReal[{0, 1000}, {numberOf2DPoints, 2}];

I'd like to find a way to quickly generate a distribution I can study for the number of points within a distance $r$ from each point, and then I'd like to select points that have at least a lowerbound $k_a$ and an upperbound $k_b$ number of points within a distance $r$ of themselves. Is there a way to use a function like Nearest to accomplish this?

Clarification --- The lowerbound $k_a$ and upperbound $k_b$ refers strictly to the count for the number of points in a circular disk of radius $r$ centered on a particular point (hopefully this makes sense). So I'd want basically a simple histogram for what this distribution of point counts looks like, and to select points that have satisfy the upper- and lowerbound point count criterion.

Answer 1

It's not easy to find in the documentation on Nearest and NearestFunction but they can return all points within a certain radius.

From tutorial/UsingNearest

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

So you can get all points that lie between a distance of 2 and 3 like so:

numberOf2DPoints = 10^6;
pointList = RandomReal[{0, 1000}, {numberOf2DPoints, 2}];
nf = Nearest[pointList];

Complement[
 nf[pointList[[31]], {Infinity, 3}],
 nf[pointList[[31]], {Infinity, 2}]]

Perhaps there is yet another way to call a NearestFunction that removes the need for Complement

Answer 2

This is the distribution you're after. Not as fast as one might want, but:

numberOf2DPoints = 10^5;
pointList = RandomReal[{0, 1000}, {numberOf2DPoints, 2}];
f = Nearest[pointList];
leuc = EuclideanDistance[#, f[#, 2][[2]]] & /@ pointList;
h[leuc_, min_, max_] := Length@Select[leuc, min <= # <= max &]
Plot3D[h[leuc, min, max], {min, 0, 7}, {max, 0, 7}, PlotRange -> All]