# Filtering a point cloud: How to reject some points in a 3D cloud of particles?

I'm having an issue with a 3D distribution of points. The distribution is really nice looking in Mathematica, but it is used in an OpenGL astronomy program (Celestia) to generate a nebula, with sprites. The problem is caused by a too high density of points at some locations in the distribution. The sprites are then adding, which turns into an ugly rendering. I need to lower the density of points at those locations only, and I really don't know how.

I need to reject some points when the local density gets too high, or to add a short range repulsive force between the particles.

Here's the code I'm currently using with Mathematica 7.0. That code was defined elsewhere on this forum (see Simulation of the Crab nebula)

MinRadius := 0.00; (* min radius of distribution *)
MinSprite := 0.003; (* min radius of sprites *)
MaxSprite := 0.07; (* max radius of sprites *)
Oblateness := 0.8; (* oblateness of the spherical distribution *)
NumberOfVoids := 1000;
NumberOfPoints := 20000;

SpriteSize[r_] :=

voidpts = Select[RandomReal[{-1, 1}, {NumberOfVoids, 3}],
pts = Select[RandomReal[{-1, 1}, {NumberOfPoints, 3}],
nf = Nearest[voidpts];
DistributeDefinitions[nf];

pts = ParallelMap[Nest[0.9975 (# + 0.01 (# - First@nf[#])) &, #, 100] &, pts];

PlotColor = ColorData["SunsetColors"];

Graphics3D[{PointSize[0.005], {PlotColor[Norm[1.3 #]], Sphere[#, 0.005]} & /@ pts},
Boxed -> False, Background -> Black, Lighting -> "Neutral", SphericalRegion -> True]


The number of points generated at the intersection of several filaments is too large : there's too much points accumulated there, and this is a problem for rendering in Celestia (too much sprites at the same location is very ugly).

So is there a way to reduce or dilute these points ? I need a rejection mecanism when the local density gets higher than a certain value (depending on the local sprite size).

The shortest distance between points shouldn't be smaller than the local sprite size (which is a function of distance to the origin). How can I add this as a constraint in the iteration process ?

Take note that I'm working with Mathematica 7.0 only (I don't have access to newer versions), so please consider suggestions compatible with this version.

• A couple of friendly suggestions: Firstly, ask just one question - do you want answers on speeding up the code, or on reducing the point density? Secondly, try not to include too much non-essential code. For example the colouring of the sprites is not really relevant to the problem, nor is the creation of the final output data for Celestia. Whilst these things are vital for your application, in the context of what you're asking here they are just clutter. Mar 11 '13 at 11:41
• Maybe some of the suggestions here could help? mathematica.stackexchange.com/questions/2594/… Mar 11 '13 at 13:13
• @Simon : Agreed. I edited a bit my first message and simplified the code. The important thing to me is to reduce the density of points in the critical locations (intersections of the "tree-like" filaments).
– Cham
Mar 11 '13 at 13:49
• @Simon : How can I add a directive about shortest distance between points = local sprite size ? Maybe there should be a constraint in the iteration process ? The sprite size is a function of the distance from the origin. Is there a way to use it as a constraint ?
– Cham
Mar 11 '13 at 14:14

Here's a pretty simple approach for filtering the point cloud. How it works:

• For each point p, find its nearest neighbour q.
• If they are closer than the minimum required separation, and p is further from the origin than q, remove p from the list. (The "further from the origin" part is to make sure that for any given pair of nearest neighbours, we only remove one from the list.)
• Repeat the procedure until the list of points no longer changes.

The minimum separation between points is a function of distance from the origin, defined by minSep, which can be adjusted to your needs.

minSep[r_] := 0.01 r

keep[{p_, q_}] := Norm[p] < Norm[q] || Norm[p - q] > minSep[Norm[p]]

filterOnce[pts_] := With[{nf = Nearest[pts]}, Select[pts, keep[nf[#, 2]] &]]

pts = FixedPoint[filterOnce, pts];

• This code appears to work very well ! I'll have to do more tests, to be sure that it's a complete solution to my problem. Apparently, it is ! 8-)
– Cham
Mar 11 '13 at 21:04