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

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

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