Counting spheres within a certain radius of a sphere in a dense packing

I have a large set of X,Y,Z coordinates which define the centre points of my spheres in a certain packing. I want to be able to count the number of spheres within a particular distance of each sphere, i.e. how many sphere's are within 1R or 2R of each sphere.

At the moment I'm just checking which other coordinates are nearby using an if loop to check if delta x, y and z are less than 1R. If all 3 deltas are within than value then I count the point.

However as I do this with a large data set I want to try and see if there is a better way of doing this.

Any Suggestions?

JP

[Walter Roberson]

There is no such thing as an if loop.

[Image Analyst]

No sphere will be within 1R or 2R of another unless they can penetrate each other, or unless you're talking about outer surface to center instead of center to center, or unless that R is not the R of the spheres but something greater than twice the sphere radius.

By the way, an interesting book on spatial statistics is here: http://www.csiro.au/resources/~/media/CSIROau/Divisions/CSIRO%20Mathematics%20Informatics%20%20Statistics/Rspatialcourse_CMIS_PDF%20Standard.pdf

Answer 1

[Sean de Wolski]

I'm going to think out loud on this one:

if you make three lists of you coordinates, each sorted by a different primary and secondary dimension (x,y,z), (y z x), (z x y) and had another list of the corresponding rows in each dimension. You could begin traversing at the row closest to the coordinate center you want, and work out in each dimension. Once a dimension has failed (e.g. y is too big, you quit traversing there, until all dimensions fail.

Edit Never mind, don't bother doing that. /edit

How many points do you have? I'd imagine a well vectorized brute-force approach could be pretty quick.

Point_we_care_about = [42 71 -30];
Max_radius = 200;
xyz_data = (rand(1000,3)-0.5)*1000;
rows_inside = sqrt(sum(bsxfun(@minus,Point_we_care_about,xyz_data).^2,2))<Max_radius;
%rows_inside is a logical vector of rows inside

find(rows_inside) %row Indices

Making xyz_data have 10000 rows it takes .0014 seconds on my system; 1/100 of a second with 100000 rows.

[John]

Thanks for that answer! Significantly faster than my if/else loop! It takes just over 1.5 secs to find all the neighbours of each individual sphere, which is a huge improvement over my 25sec loop.

Next step is to set it to run for every timestep interval rather than just one. That may take a while!!

Only problem with it is that it lists the object row in the answer as well. But that can easily be removed.

[Sean de Wolski]

Yes, I would do that at the end though if you can. Resizing arrays takes time.

[John]

Do you think it would be faster to include it with the max radius (i.e 0 > x < Max_radius)?

[Sean de Wolski]

Yes. You could always store the radii and the indices of the rows corresponding to small enough radii, then extract those radii and remove the one equal to 0 quickly and painlessly later.

Not knowing your code structure it's hard to say.

Answer 2

[Walter Roberson]

A close relative of this question came up several times last fall (almost as if it had been assigned as a class project at the time!). The question then had to do with building simulators for soil deposition, through either water or wind action. As the deposition was random, the placement of new particles depended upon the ability to detect that there were no existing particles at the proposed new site.

Unfortunately in those series of posts, no-one had a proposal that was faster than basic trial-and-error. Except that I proposed that octrees could make that particular task more efficient.

I have kept the problem simmering in the back of my mind since, hoping that some day I might come up with a nifty mathematical transformation that would make it much faster, but unfortunately I have not come up with one.

In your current question, are you doing this probing for a small number of points, or for all (or the majority of) points? It might not be out of the question that a "sliding window" style approach might help.

Alternately, perhaps a vectorized calculation of the distance from every sphere to every other might work better than one might expect, as it could be largely BLAS'd and you could skip the sqrt() step. repmat() the squares of the individual sphere radii (if the 1R and 2R vary per sphere), sum(A < B), subtract 1 because the self-distance is 0... might turn out to be efficient enough, perhaps much more efficient than trying to be "smart" about it.

[John]

What I'm trying to do initially is to create a list of all the neighbouring spheres (those in contact with sphere x). This will allow me to look at the distribution of contacts within a packing structure (and how it evolves over time).

I also want to be able to use it to calculate the number of spheres within a certain radius (say 10R) and then use this for calculating values of void ratio/porosity, which I can average then over a grid.

Answer 3

[Image Analyst]

The stats program R has a plug-in called spatstat. spatstat has a function called pairdist that calculates the distances. There are a bunch of other types of metrics it also uses to characterize the spatial statistics. Maybe if you get that package it might have the source code. See chapter 19 of this book: http://www.csiro.au/resources/~/media/CSIROau/Divisions/CSIRO%20Mathematics%20Informatics%20%20Statistics/Rspatialcourse_CMIS_PDF%20Standard.pdf Even if you can't get source code, this book is an excellent resource for all kinds of ways to analyze spatial distributions for "regularity," spatial uniformity, clumping, etc.

[John]

Thanks a lot for this. R looks quite a useful package to have so I'm going to explore it later in the week (or over Christmas). I like that there is a pair correlation function available, hopefully it will be significantly faster than my MATLAB routine, which struggles above a few thousand points.