Algorithm to Find the Aggregate Mass of "Granola Bar"-Like Structures?

Question

I'm a planetary science researcher and one project I'm working on is N-body simulations of Saturn's rings. The goal of this particular study is to watch as particles clump together under their own self-gravity and measure the aggregate mass of the clumps versus the mean velocity of all particles in the cell. We're trying to figure out if this can explain some observations made by the Cassini spacecraft during the Saturnian summer solstice when large structures were seen casting shadows on the nearly edge-on rings. Below is a screenshot of what any given timestep looks like. (Each particle is 2 m in diameter and the simulation cell itself is around 700 m across.)

The code I'm using already spits out the mean velocity at every timestep. What I need to do is figure out a way to determine the mass of particles in the clumps and NOT the stray particles between them. I know every particle's position, mass, size, etc., but I don't know easily that, say, particles 30,000-40,000 along with 102,000-105,000 make up one strand that to the human eye is obvious.

So, the algorithm I need to write would need to be a code with as few user-entered parameters as possible (for replicability and objectivity) that would go through all the particle positions, figure out what particles belong to clumps, and then calculate the mass. It would be great if it could do it for "each" clump/strand as opposed to everything over the cell, but I don't think I actually need it to separate them out.

The only thing I was thinking of was doing some sort of N² distance calculation where I'd calculate the distance between every particle and if, say, the closest 100 particles were within a certain distance, then that particle would be considered part of a cluster. But that seems pretty sloppy and I was hoping that you CS folks and programmers might know of a more elegant solution?

---

Edited with My Solution: What I did was to take a sort of nearest-neighbor / cluster approach and do the quick-n-dirty N² implementation first. So, take every particle, calculate distance to all other particles, and the threshold for in a cluster or not was whether there were N particles within d distance (two parameters that have to be set a priori, unfortunately, but as was said by some responses/comments, I wasn't going to get away with not having some of those).

I then sped it up by not sorting distances but simply doing an order N search and increment a counter for the particles within d, and that sped stuff up by a factor of 6. Then I added a "stupid programmer's tree" (because I know next to nothing about tree codes). I divide up the simulation cell into a set number of grids (best results when grid size ≈7 d) where the main grid lines up with the cell, one grid is offset by half in x and y, and the other two are offset by 1/4 in ±x and ±y. The code then divides particles into the grids, then each particle N only has to have distances calculated to the other particles in that cell.

Theoretically, if this were a real tree, I should get order N*log(N) as opposed to N² speeds. I got somewhere between the two, where for a 50,000-particle sub-set I got a 17x increase in speed, and for a 150,000-particle cell, I got a 38x increase in speed. 12 seconds for the first, 53 seconds for the second, 460 seconds for a 500,000-particle cell. Those are comparable speeds to how long the code takes to run the simulation 1 timestep forward, so that's reasonable at this point. Oh -- and it's fully threaded, so it'll take as many processors as I can throw at it.

Answer 1

My first suggestion is to cut up your problem into two problems: first, figure out what you want and then figure out how to efficiently get what you want. You can't efficiently get something you haven't defined yet. I'll put some ideas in this answer that may help you find this definition. I suggest you make an inefficient implementation of the ideas you like first, apply it to a few not-too-large datasets, evaluate the results by hand, adapt your definition and repeat (possibly asking another question here), until you are happy with your definition. After that, I suggest you ask another question on how to efficiently compute the result of your definition (if you still need help).

So, let's see what would correspond with our intuitive idea of a 'strand'. Your strands seem to consist of roughly uniformly distributed points, though you should check this by making a zoomed-in picture (of the original dataset) - the resolution of your picture is too low to say with certainty that the points really are roughly uniformly distributed. I'll assume they are for this answer.

An initial idea may be to look at the nearest neighbor of every point. Let's pick a point X, call its nearest neighbor Y and set D as the distance between X and Y. We then look at the circle C around X with radius D*A, where A is a tuning parameter, say A=3. If X is part of a strand, we expect that for every point Z in C, the distance from Z to its nearest neighbor W is about the same as D. If it is significantly shorter, say more than A (or maybe some other parameter B) then X is apparently near to points that are much nearer to each other than to X, so X is probably not part of a strand.

This criterion is not complete however. It only gives a criterion to detect a 'border' between between areas dense with points and areas less dense with points. We still have to cluster points together into strands.

There is a feature in your picture that shows that this is not simple. In the lower-right corner of your picture, there is a relatively large area with a lot of stray points. These stray points are themselves roughly uniformly distributed, so if we were to remove all points in the strand around it (and all other points), then we would expect any strand-detecting algorithm to mark this set of stray points as a strand! We therefore need to be careful when making our clusters.

An idea may be to do the following. We are going to make a graph on these points, where the vertices are the points and edges signify that two points have similar density. For every point, we check the criterion above. If it checks out, we connect X with an edge to all points in C. If it does not check out, we do not add any edge, and mark X as 'stray'. After doing this for every point, we consider the set of connected components. These should consist of a single (in the case of your picture, but other datasets may have multiple) connected component consisting of all points in strands, plus (potentially a lot) more components consisting of single stray points and these 'stray strands'. However, these stray strands have points in them that have been marked as 'stray', so you can simply ignore any component containing a point that has been marked as 'stray'.

A danger of this idea is that you may have a feature where the density of a strand progressively lowers as you move along the strand, until the density is so low that it is just a set of stray points. As our criterion is 'local', it may not detect this and mark these stray points as part of the strand. I'm not sure if this will be a problem: I'd guess most stray points should be caught by the criterion, as changes in density seem quite abrupt in your picture.

If this problem does occur, you could try an alternative to just taking the connected components. For every point X, we calculate the distance to its nearest neighbor D(X). We start at the point with minimal D(X) and perform a BFS (or DFS, the order doesn't matter). We add any point Y whose D(Y) is not much larger than the D(X) (by a tunable factor) we started with. If we do encounter a point Y that has too large D(Y), we remove the edge (X, Y), mark Y as 'stray' and act as if we never visited Y in our BFS. If tuned right, this should prevent the problem I describe above.

An alternative idea to fix this problem acts a bit more local: you could do a BFS and keep track of the lowest D(X) (I use D(X) as a measure of the density around a point) encountered at most say 10 BFS-steps before, and if we encounter a Y that has D(Y) that is much larger than this D(X), we do the same as the other (potential) solution I offered.

As a disclaimer: all the ideas above I thought up on the spot just now, I don't really know if this particular problem has been studied before, so I may just be sprouting nonsense. Just try out the ideas (whether my ideas or your own) that sound sensible to you and find out if they really work, and only then focus on implementing them efficiently.

Answer 2

Using modular decomposition you can create a tree that will contain all particles as leafs and upper nodes will cluster these. Based on that tree you can define measures that are applied to every node of it from the root to the leafs downwards. You stop this downwards traversal when the measurements reach user defined thresholds. One such measurement may be the density of the convex hull of all particles in a cluster.

Answer 3

I think you're after a machine learning clustering algorithm.

This page from the Python SciKit Learn toolkit has pictures that suggest the DBSCAN algorithm (Wikipedia) could be what you're looking for. It seems ideal as it's input parameter is the neighbourhood size, while most other clustering algorithms want the number of clusters, which you would not know in advance.

"A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise" Ester, M., H. P. Kriegel, J. Sander, and X. Xu, In Proceedings of the 2nd International Conference on Knowledge Discovery and Data Mining, Portland, OR, AAAI Press, pp. 226–231. 1996

Answer 4

I have been thinking about this problem. I am not a physics expert, so bear with me.

It seems that it is not the distance between the particles that counts to determine the clumps. It is whether or not the gravity fields overlap.

Take a particle P, and determine what other particles have overlapping gravity fields.

Then take one of those and do the same thing. Your goal is not to find all the particles in the clump but to find its boundaries.

Repeat this until all the clumps are found.

Now go back and determine the mass of the clumps. You will have eliminated stray particles, and you can use the clump boundaries to find the mass.

I am not sure if this helps, but it is all I could think of.

Answer 5

You could, at the end of each timestep, convert the data into a graph, compute the minimum spanning tree, and then begin to remove edges that exceed a certain threshold. That should give you clumps and an easy way to enumerate through the particles in each clump.