I'm writing a lesson for students on gravity erosion of terrain. The students code a short program in the Wolfram language that turns an image into a 3D model of the erosion. In earlier lessons they coded 2D and very low-resolution 3D models. Next I'd like them to make a higher resolution 3D model.

Here is the code to make the 252x252 grayscale test image:

let1 = Rasterize[Style[" \[ClubSuit] ", FontFamily -> "Arial Black", 270, Black], ImageSize -> 360];
let2 = Rasterize[Style[" \[SpadeSuit] ", FontFamily -> "Arial Black", 270, Gray], Background -> None, ImageSize -> 360];
img = ImagePad[ImageCrop[ColorNegate[ColorConvert[ImageCompose[let1, let2, {0, 0}, {-90, 75}], "Grayscale"]], {180, 180}], 36]

And here's the image:

Here is the code that erodes the image:

dim = ImageDimensions[img];
ht = 48;
dirt = .33;
scree = .66;
gens = 300;
orig = ht*ImageData[ImageReflect[img]];
erode[ar_] := (
   dat = ReleaseHold[ar];
   in = Total[Flatten[Map[(If[#[[1]] > dat[[2, 2, 1]] + scree && #[[2]] > 0, dirt/#[[2]], 0] &), dat, {2}]]];
   out = If[dat[[2, 2, 2]] > 0, dirt, 0];
   Return[dat[[2, 2, 1]] + in - out]);
new = Nest[(
     falls = ArrayFilter[(Count[Flatten[#],n_ /; n < #[[2, 2]] - scree] &), #, 1];
     ArrayFilter[erode, MapThread[Hold[{#1, #2}] &, {#, falls}, 2], 1]) &,
   orig, gens];

Basically, it is looking at every pixel in the height field and removing dirt that falls away while adding dirt that falls in from neighboring pixels. The model it produces has all the terrain features I want the students to see (scree slopes, mantled peaks, V-shaped valleys, etc.). It also preserves the quantity of dirt, so the model has the same amount of dirt (or height sum) at any given iteration. Here's a screen shot:

The problem is that the erode function takes way too long. On my fairly powerful Mac it takes about 12 minutes to produce the data for the model you see above. That's not feasible for a lesson. I've tried reducing the size of the image, but the model becomes too blocky. Fewer iterations (gens) means less erosion. A larger value for dirt causes uneven erosion and spikes on the slopes. I think it can be streamlined either in Nest or in erode, but after two days of trying, I'm unable to optimize it. Can anyone please help me out on this one?

Thanks in advance,

Mark Greenberg