Counting number of points in fluctuating lines that fall within a circle

Question

I have a list of points pts that fluctuate around some mean line. pts is pre-sorted by the projection of the points onto this mean line (for instance, in the special case where the points fluctuate about a horizontal line, then the points would be sorted by their x-values). I want to find the number of elements in pts that land within a circle centered at fixed point centre and for various radii.

My problem is that my current code is too slow, as I work with systems with a huge number of points, and I feel that there must be a way to further exploit the geometry of the problem to get significant speed-ups.

Here's some sample data for pts, which I'll have fluctuate about the mean line defined by meanY and unit vector vec:

(*The numbers I use in this block are not important*)

vec = {0.951057, 0.309017};
orthoVec = {-vec[[2]], vec[[1]]}; (*the normal vector*)
meanY = -10;

px = FoldList[Plus, RandomReal[{0.5, 1}, 100]];
px = px - Median[px]; (*just for some symmetry in the data*)
py = RandomReal[meanY + {-0.5, 0.5}, 100];
pts = Transpose@{px, py};
pts = pts . {vec, orthoVec}; (*This line is the change of basis so that the data is along our mean line*)

and the circle parameters:

centre=RandomReal[{-1,1},2];
radii=Table[r, {r, 1., 40., 0.5}];

This gives a configuration that looks like this:

Here's my current attempt: I find the norm of all the points from the circle center, sort them, and them finding the index of the largest norm less than the circle radius. To find the indices, I use the function cf from Henrik Schumacher's answer here. The code for my attempt is:

norms = Sqrt[(pts[[All,1]] - centre[[1]])^2 + (pts[[All,2]] - centre[[2]])^2]; (*take the norm of the points*)
norms=Sort[norms]; (*sort the points*)
ctr = cf[norms, radii]; (*count the number of points inside the circle for each radius in radii*)

But given the structure of my point set (i.e., along a line), this feels like a waste, primarily because sorting the norms is very costly when I already have a point set that is sorted along a line. Any ideas on how I can exploit this fact to speed up my code?

Here is my own start of an idea on how this algorithm might be improved: because the points in pts are already sorted by their projections onto the mean line, we get the following plot for the norms by index:

This leads me to think that there should be a way to first rough sort the points according to just some function of their index, and then apply Sort to the resulting list, which I think should result in fewer computations, but I'm not sure how to go about coding this up.

Edit

Thanks for the comments everyone. Here's some more info about what I'm trying to do that will hopefully clarify my question.

The points fluctuating about a line is just one piece of my data. What I'm actually trying to do is to compute the number of points inside the circle from a 2D point system as a function of circle radius, where my point system can be broken into these fluctuating lines:

ptsArray = Table[
   px = FoldList[Plus, RandomReal[{0.5, 1}, 100]];
   px = px - Median[px];
   py = RandomReal[meanY + {-0.5, 0.5}, 100];
   pts = Transpose@{px, py};
   pts.{vec, orthoVec}, {meanY, -20, 20, 2.5}];

A more complete picture of my point system would look like this:

I COULD have just used all of the data combined, computed the norms, and found the number of points in the circle that way, but it slows down the code significantly when we go to very large systems (on the order of millions of points in all) and I've found it much more efficient to break it into these 1D sequences. That is why the code needs to be so fast: I'm not just computing this for 1 line of points, I'm looping through thousands of such lines and then adding up the number of points in the circle from each.

Answer 1

I might be very well not understanding the issue because with even 100,000 points obtaining the number of points within all possible radii doesn't seem to take much time:

n = 100000;
px = FoldList[Plus, RandomReal[{0.5, 1}, n]];
px = px - Median[px]; (*just for some symmetry in the data*)
py = RandomReal[meanY + {-0.5, 0.5}, n];
pts = Transpose@{px, py};
pts = pts . {vec, orthoVec}; 

centre = {0, 0};
AbsoluteTiming[
  distance = Sqrt[(Total[(# - centre)^2] & /@ pts)] // Sort;
  index = Range[Length[distance]];]
(* {0.0292905, Null} *)

ListPlot[Transpose[{distance, index}], Frame -> True,
 FrameLabel -> {"Radius", "Number of points within radius"}]

One can then turn that into a function:

points = Interpolation[Transpose[{distance, index}]]
Floor[points[5000]]
(* 13337 *)

Answer 2

Perhaps the cleanest way to find this distribution is with Nearest. I.e.

pts=RandomReal[{-10,10},{10000000,2}];
center={1,1};
rad=1;
Timing@Length@Nearest[pts,center,{\[Infinity],rad}]
(* .08 seconds, 78624 *)
Timing@Length@Select[pts,Norm[#-center]^2<rad&]
(* upwards of 40 seconds, 78624 *)

with the added bonus that if you omit the preceeding Length you have a list of points ordered by closeness to center. Not sure yet how to calculate the variance of the distribution Nearest[pts,center,10000], but I'm sure it can be done. Regardless, we can plot it (this time with a uniform grid, idk how that should affect the variance):

ListPlot[
  Norm[#-{1,1}]&/@ 
    Nearest[
      Join@@CoordinateBoundsArray[{{-10,10},{-10,10}},0.0063210],
    {1,1},{\[Infinity],1}]
]

Answer 3

EDIT: I've updated the answer to use BinCounts, which suits this problem in multiple ways.

Most of the old description below applies to this solution, but here the counts for intervals between (squared) radii is performed with BinCounts, and they're summed up to correspond with the UnitStep solution with Accumulate:

Needs["Developer`"];
With[{length = 1000000},
 With[
  {pts = ToPackedArray[RandomReal[{-30, 30}, {length, 2}]],
   centre = RandomReal[{-1, 1}, 2]},
  Accumulate[
   BinCounts[(pts + ConstantArray[-centre, length])^2 . {1, 1},
    {Prepend[Range[1., 40., 0.5], 0.]^2}]]]]

The actual computation takes about 0.10 seconds for 79 radii and million points on my laptop, or 8.2 seconds for 100 million points(!). (Using packed arrays gives just a couple-percent benefit on calculation side but it's so easy to add that I did it.)

EDIT: This can be sped even further!

If we use regular binning, the binning of the large squared distance list is even faster. Thus we first perform bin counts for bins between zero and square of the largest radius, with a bin size of square of the step, and then re-bin this data using a second HistogramList run (which is a cousin of BinCounts, but also handles WeightedData).

fixedHistogramList corresponds to HistogramList, but fixes a slow implementation for specific inputs - maybe in the future it is unnecessary.

Needs["Developer`"];

(* Fix an odd performance issue which is present at least up to v13.1.0. *)
ClearAll[fixedHistogramList];
fixedHistogramList[list_List, {xmin_?NumericQ, xmax_, dx_}] :=
  {Range[xmin, xmax, dx], BinCounts[list, {xmin, xmax, dx}]};

With[{binrange = {1., 40., 0.5}},
 With[{length = 1000000},
  With[{
    pts = ToPackedArray[RandomReal[{-30, 30}, {length, 2}]], 
    centre = RandomReal[{-1, 1}, 2]},
   fixedHistogramList[(pts + ConstantArray[-centre, length])^2 . {1, 1},
     {0, binrange[[2]], binrange[[3]]}^2] //
    Accumulate[Round[Last[
        HistogramList[WeightedData[Most@First@#, Last@#],
         {Prepend[Range @@ binrange, 0.]}^2]]]] &]]]

The actual computation (that is, excluding generation of the pts list) takes 3.6 seconds on my laptop for 100 million points and 79 radii (end result bins).

This is very slightly faster than putting Sqrt around the vectored squared distance calculation and dropping ^2s from range definitions. The square-rootless version has the benefit of avoiding computing million square roots (although they can be quite fast in a vectored manner), while one with a square root has a benefit of having less "unnecessary" bins on regular binning step.

---

Old answer:

I simplify the problem and don't care about the lines - just those random points and circles of different radii sharing the same centre.

This computes counts for a single set of points with different radii:

With[{length = 1000000},
 With[
  {pts = RandomReal[{-30, 30}, {length, 2}],
   centre = RandomReal[{-1, 1}, 2]},
  With[{sqdistances = ((pts + ConstantArray[-centre, length])^2) . {1, 1}},
   Table[Total[UnitStep[r^2 - sqdistances]], {r, 1., 40., 0.5}]]]]

Here centre is subtracted from all point coordinates, individual resulting $x$ and $y$ values squared with a vectored way, and sums computed using .. This results squared distances from centre. After this UnitStep is used to convert each distance inside r to $1$, outliers $0$ and a total of these is computed.

Computing these counts for 79 radii with million points takes about 0.6 seconds on my laptop.

When number of radii being tested is large, this computation is pretty much dominated by the counting part, in which the UnitStep approach seems to be unfortunately dramatically more efficient that Select or Count approaches for this sort of a problem. So, the following, easier-to-understand approach is only marginally slower in the end:

With[{length = 1000000},
 With[
  {pts = RandomReal[{-30, 30}, {length, 2}],
   centre = RandomReal[{-1, 1}, 2]}, 
  With[{distances = Norm[# - centre] & /@ pts}, 
   Table[Total[UnitStep[r - distances]], {r, 1., 40., 0.5}]]]]

Answer 4

To keep things simple I give an example for just one radius. You may adapt it for several radii yourself.

First we create the points and the circle according to your example:

SeedRandom[1];

px = FoldList[Plus, RandomReal[{0.5, 1}, 100]];
px = px - Median[px]; (*just for some symmetry in the data*)
py = RandomReal[meanY + {-0.5, 0.5}, 100];
pts = Transpose@{px, py};
pts = pts . {vec, orthoVec};

circle = Circle[{0, 0}, 15];

Graphics[{Point[pts], circle}, Axes -> True]

Then we fit a line through the points:

line = LinearModelFit[pts, {x}, x]["BestFit"]

Then subtract from all points the y intercept of the line: pts1. As the points now go through the origin we can rotate them so that they lay along the x axis: pts1:

pts[[All, 2]] -= line[[1]];
pts1 = RotationMatrix[-line[[2, 1]]] . Transpose@pts // Transpose;

We also shift and rotate the circle: circle1:

circle1 = circle;
circle1[[1]] = 
 RotationMatrix[-line[[2, 1]]] . (circle1[[1]] - {0, line[[1]]});

This now looks like:

Graphics[{Point[pts1], circle1}, Axes -> True]

Next we determine the points where the transformed circle cuts the x axis. These points determine approx. the points to count. As the points a fuzzy, we add some margin for safety:

limit = Sqrt[circle1[[2]]^2 - circle1[[1, 2]]^2] + 1;

Possible points to count lay inside the limits:

possible = Select[pts1, (-limit <= #[[1]] <= limit) &];

Now because the points do not lay on a straight line we need to check and count more carefully:

Select[possible, Norm[circle1[[1]] - #] <= circle1[[2]] &] // Length

(* 27 *)

And for convenience, here is all together:

SeedRandom[1];

vec = {0.951057, 0.309017};
orthoVec = {-vec[[2]], vec[[1]]}; (*the normal vector*)
meanY = -10;

px = FoldList[Plus, RandomReal[{0.5, 1}, 100]];
px = px - Median[px]; (*just for some symmetry in the data*)
py = RandomReal[meanY + {-0.5, 0.5}, 100];
pts = Transpose@{px, py};
pts = pts . {vec, orthoVec};

circle = Circle[{0, 0}, 15];

line = LinearModelFit[pts, {x}, x]["BestFit"];
pts[[All, 2]] -= line[[1]];
pts1 = RotationMatrix[-line[[2, 1]]] . Transpose@pts // Transpose;
circle1 = circle;
circle1[[1]] = 
  RotationMatrix[-line[[2, 1]]] . (circle1[[1]] - {0, line[[1]]});
Graphics[{Point[pts1], circle1}, Axes -> True]
limit = Sqrt[circle1[[2]]^2 - circle1[[1, 2]]^2] + 1;
possible = Select[pts1, (-limit <= #[[1]] <= limit) &];
Select[possible, Norm[circle1[[1]] - #] <= circle1[[2]] &] // Length