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

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.

• (1) I would use RandomReal instead of Random. (2) Not really clear how "erode" is intended to operate. You might want to elaborate on that. Is the idea to pick a point at random and remove it along with several of its closest neighbors? Feb 6, 2014 at 14:47
• Eroding mean randomly removing points, yes. How can I ask Mma to remove a random point and several of its neihbors, and do it again for many other randomly selected points ?
– Cham
Feb 6, 2014 at 14:54
• According to the eroding process I defined above, the final distribution may contain holes inside. This is an interesting possibility for what I'm trying to achieve. I'm also interested in an eroding process that act "from the exterior" only (no "holes" inside), but this is more ambiguous to me ; I don't know how to define the "surface" of a distribution of points. Select first a point which is farthest from the center (origin) ?
– Cham
Feb 6, 2014 at 15:06

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.

• I didn't had time to try this solution yet, but is it eroding the ball from the outside, or is it also making holes inside ?
– Cham
Feb 6, 2014 at 19:16

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.

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]