How to randomly erode a 3D distribution of points with Mma 7.0?

Question

Suppose we have a 3D distribution of points, like the spherical ball example below :

Ball[num_]:=Table[
{
    #1 Sqrt[1-#2^2]Cos[#3],
    #1 Sqrt[1-#2^2]Sin[#3],
    #1 #2
}

&[
    Random[NormalDistribution[1, 0.5]],
    Random[Real,{-1,1}],
    Random[Real,{0,2Pi}]
],{num}]

Graphics3D[{AbsolutePointSize[2],Point[Ball[10000]]},Boxed->True,BoxRatios->{1,1,1},ImageSize->800,SphericalRegion->True]

This code produces the ball shown here :

Obviously, this ball has a spherical symmetry and could be considered as a "smooth" distribution (under the limit of an infinite set of points).

Now, I would like to randomly "erode" the distribution, to get something that look a bit like a fractal shape (think about how an eroded mountain could be obtain from a smooth hill). How can I do that ? What Mathematica procedures could modify the Ball distribution that act like a random eroding process ?

To me, "eroding" mean removing a randomly selected point, and several of its neighbors, then repeat the process several times.

Please, the suggestions should be compatible with Mathematica 7.0.

Answer 1

FindClusters might give you a starting point. For simplicity, let's start with your definition for Ball:

Ball[num_] := 
 Table[{#1 Sqrt[1 - #2^2] Cos[#3], #1 Sqrt[
       1 - #2^2] Sin[#3], #1 #2} &[
   RandomReal[NormalDistribution[1, 0.5]], Random[Real, {-1, 1}], 
   Random[Real, {0, 2 Pi}]], {num}]

By using the function FindClusters with the Method -> "Agglomerate" you can split the points into clusters (I chose to split into 500 clusters, but you can tune this to your preferences). After a sort which puts the largest cluster first (which will likely be the cluster of particles near the center) you can "erode" by included the first n components you desire.

pts = Ball[10000];    
clusters = 
  SortBy[FindClusters[pts, 5000, Method -> "Agglomerate"], 
   1/Length[#] &];

For example I can "erode" by about 50% by including the first 500 components

N[Total[Map[Length,clusters[[1;;500]],2]]/10000]
(* 0.4674 *)

GraphicsRow[{ListPointPlot3D[pts, 
   PlotStyle -> Directive[Black, AbsolutePointSize[2]], 
   BoxRatios -> 1, PlotRange -> {{-2, 2}, {-2, 2}, {-2, 2}}, 
   Boxed -> False, Axes -> False, PlotLabel -> "Original"], 
  ListPointPlot3D[clusters[[1 ;; 500]], 
   PlotStyle -> Directive[Black, AbsolutePointSize[2]], 
   BoxRatios -> 1, PlotRange -> {{-2, 2}, {-2, 2}, {-2, 2}}, 
   Boxed -> False, Axes -> False, PlotLabel -> "Eroded"]}, 
 ImageSize -> 800]

There is likely quite a bit to improve, but again this might provide a good initial guess for your particular needs.

Answer 2

I'll work with a built-in for the random point generation. If this is not in version 7 you can still use RandomReal.

ball[n_] := RandomVariate[NormalDistribution[], {n, 3}]

bl = ball[10^4];

Graphics3D[{AbsolutePointSize[2], Point[bb]}, Boxed -> True, 
 BoxRatios -> {1, 1, 1}, SphericalRegion -> True]

Now we can erode clumps of points at random as follows.

erode[pts_, holes_, hsize_] := 
 Module[{nf = Nearest[pts], pt, nbrs, rmove},
  rmove = Reap[Do[
      pt = pts[[RandomInteger[{1, Length[pts]}]]];
      nbrs = nf[pt, RandomInteger[{1, hsize}]];
      Map[Sow, nbrs];
      , {holes}]][[2, 1]];
  DeleteCases[pts, Alternatives @@ rmove]
  ]

erodedball = erode[bl, 100, 100];

Eroding from the outer regions first is slightly trickier. You might first invert all points e.g. dividing by distance squared from origin, so that the furthest are now near the center. Randomly choose from those, and use erode in neighborhoods of the original values. For funkier effects, perhaps erode from the inverted set, so "neighbors" are not really neighbors in the Euclidean sense.

Answer 3

To visualize the deleting process.

(* generate random data*)
data = RandomVariate[NormalDistribution[1, 3], {5000, 3}];
(* generate random delete centers*)
pt = RandomChoice[data, 30];
(* generate random outer delete centers say at least away from mean by 10 *)
ptouter =RandomChoice[Complement[data, Nearest[data, Mean@data, {Infinity, 10}]], 30];
f = Nearest[data, #, 100] &;
removes = FoldList[Join, {}, f /@ pt];
removesouter = FoldList[Join, {}, f /@ ptouter];
(* list of the rest points after each delete *)
afterRemove = Complement[data, #] & /@ removes;
afterRemoveouter = Complement[data, #] & /@ removesouter;

If you use ListPointPlot3D and ListSurfacePlot3d we get this cool animation. We highlight the points to be deleted during the current step with the red sphere denoting the active delete center.

Separately looking at the two cases now in the following first image we can animate the process for random delete in the cluster. Next one is the the animation when we delete points from the outer region of the cluster.

Code:

pics = Table[
   Graphics3D[{{{PointSize[Medium], Directive[Opacity[.9], White], 
       Point /@ afterRemove[[1]], PointSize[Medium], 
       Directive[Opacity[.5], Black], 
       Point /@ afterRemove[[i]]}}, {Directive[Red], 
      Sphere[#, .51] &@pt[[i]]}, {Directive[Opacity[.8], Yellow], 
      Thin, Tube[Take[pt, {1, i}]]}}, SphericalRegion -> True, 
    ImageSize -> 800], {i, 1, -1 + Length@afterRemove, 1}];
Export["path"<>"try.gif", pics, 
Background -> None, "Interlaced" -> True, "DisplayDurations" -> .75]