Triangle Plasma

A couple of months ago, out of boredom, I implemented the old plasma fractals algorithm that is emblematic of the late 1980s/early 1990s demoscene, “rectangle subdivision fractals” or whatever you want to call it. It’s an old hack. I wrote some code that would generate animated GIFs that were cross-fades of plasma fractals but the GIF files came out too large and I lost interest before figuring out how to crank up the compression. It got me thinking, however, about an analog of the diamond-square algorithm that does not use rectangles.

Diamond-square is a modification of the old algorithm I implemented. It is supposed to not exhibit artifacts on grid lines the way that the simpler version does. (You can see these artifacts in the upper left in the this image from the Wikipedia article — which looks to me like it was generated with the simpler version, actually.) I’ve always thought the diamond-square algorithm is kind of inelegant because the version in which you only displace rectangle centers and just interpolate rectangle sides does not in my opinion produce output that is that much worse and is trivial and fun to implement.

In any case, if you really want zero of these kinds of artifacts then don’t use a grid at all. Instead of using a recursive hierarchy of subdivided tilings of a rectangle, my idea is to use a recursive hierarchy of Delaunay triangulations. I implemented this algorithm last week, initially exactly as below.

Generate a small set of random vertices V. Assign a random value to each vertex in V.
Calculate the Delaunay triangulation T of V.
For each triangle t in T:
- If the area of t is greater than some threshold, randomly generate a point p in t.
- Assign a value to p that is the average value of t‘s vertices weighted by their distance from p displaced by a random amount that is proportional to the relative area of t.
- Add p to V.
If any new points were added to V goto 2. otherwise terminate.

The above works but is horribly inefficient. The slow step is 2., generating full Delaunay triangulations at each iteration. The only thing the above algorithm has going for it is ease of implementation. It allows you to treat a typical Delaunay triangulation library call as a black box.

To achieve greater efficiency you need to delve into the inner workings of the triangulation process. Specifically if you have an incremental implementation of Delaunay triangulation you can modify the above to not need to re-triangulate everthing at each iteration.

I chose a modified version of the Bowyer-Watson algorithm for this purpose. Briefly Bowyer-Watson goes like this: given a set of vertices, iteratively add each vertex v to an in-progress triangulation by finding all triangles with circumcircles that contain v, deleting those triangles, inserting v, and wrapping v with triangles connecting v to each edge of the star-shaped polygon that deleting the triangles created. Done naively this is O(n^2). Done while maintaining the triangle adjacency graph and using some other data structure to efficiently find the triangle containing a given vertex, it is O(n log n). There is actually some confusion online about the time complexity of the Bowyer-Watson algorithm, (see my answer to this StackOverflow question), but it is straight-forward in my case because at every step I already know which triangle I randomly generated a point in — so I do not need some extra data structure to find triangles efficiently.

The modified algorithm is as below:

Generate a small set of random vertices V. Assign a random value to each vertex in V.
Calculate the Delaunay triangulation T of V.
while the area of t, the largest triangle in T, is greater than k do
1. generate a random point in t. Assign a random value to it as above.
2. perform a Bowyer-Watson iteration step on p in t, i.e. searching the triangle adjacency graph starting at t for all triangles with circumcircle’s containing p, deleting them. etc. This set of “bad triangles” is guaranteed to be contiguous.

The above is very fast. My implementation is here. I wrote it in C# to .NET Framework so that I would get GDI graphics routines for free. This implementation is indebted to Rafael Kuebler’s Delaunay/Voronoi code which I modified for my purposes and Renaud Bédard’s implementation of uniform poisson-disk sampling (I found that if the initial points are uniformly distributed there can be little cusps between seed points that are too close together — poisson-disk sampling yields better results) Also I use Mono.Options for command line argument parsing.

It’s a command line program with options as follows. When I mention the area of a triangle below, area is always represented as a percentage of the area of the mean of the area of the initial seed triangles:

w : Width of output (pixels)
h : Height of output (pixels)
i : Initial vertex distance. This a parameter going to the poisson disk sampler for the creation of seed points. It’s a percentage of min(w,h) defining a distance between sampled points. Higher percentages yield fewer seed points. Lower percentage yield more, which in turn yields more perlin-noise-like output
m : Minimum area cutoff threshold — area of the smallest of triangle that the program will continue subdividing i.e. the largest triangle allowed in the output .
c: contrast factor. Apply sigmoid contrast to the output values. I use a similar function to the one used by ImageMagick. High numbers mean more contrast. 10 is a lot of contrast.
v : Perturbation method, an enumerated type parameter declaring the function to use to displace triangle mid-point values. Can be “Normal”, “Uniform”, “ClampedNormal”, or “ClampedUniform”. The default is Normal, meaning given the weighted mean value of some triangle t‘s vertices v generate a new value by randomly sampling a normal distribution with mean v and with standard deviation that is the area the t raised to the power of the perturbation parameter. Uniform is similar. Clamped perturbation methods cap values at 1.0 and disallow negative values. Unclamped methods do not but I normalize all values before generating the output. Clamping tends to increase the perceived contrast in the output. I prefer using unclamped perturbation and then playing with the contrast directly, however.
p : Perturbation parameter, exponent used by the perturbation method function above.
s : Scale factor applied to output. Useful for generating output that has triangles that are effectively subpixel in size. Generate large and then scale down.
b : Color-blend, if specified generates colored output. See examples below. Otherwise the output is gray scale.

Below is some output along with the command line arguments that generated it. The program will also generate SVG based on the file extension. The first two examples are generated at high resolution and then scaled down so that the actual triangles are subpixel in size. The latter ones are color blends with the triangles showing.

-w 2400 -h 2400 -s 0.25 -i 0.6 -m 0.000005 -c 10