Generating Fractals with LINQ

Jump to the code

After seeing Luke Hoban's excellent LINQ ray tracer, I wondered if the same techniques could be used in other algorithmic image generators - specifically in generating the Mandlelbrot Set fractal.

I've never written a fractal generator before, so I used the Wikipedia article on the Mandelbrot Set, and some Java sources I found on the web as a base. Using a complex number class (strangely .NET does not provide its own) the algorithm is fairly simple.

Escape Time Algorithm

The simplest method of generating the image uses something called the escape time algorithm. Simply put, it iterates over every pixel in the image and calculates a number which is used to select a colour to draw that pixel - the number comes from repeating a calculation in the complex number plane and determining if during each repetition it satisfies an 'escape' condition - hence the name. Points within the mandelbrot set do not 'escape' (or escape to infinity), and so after a maximum number of iterations is hit the point is coloured black.

The maximum number of iterations is defined by the programmer, and determines the detail level of the fractal. The more iterations allowed, the better quality the output.

It's probably easier to describe in code.

public IEnumerable<PixelData> GenerateSet(int width, int height, int maxIterations)
{
  List<PixelData> result = new List<PixelData>();

  double xmin = -2.5;
  double xmax = 1;
  double ymin = -1.25;
  double ymax = 1.25;
  double dx = (xmax - xmin) / width;
  double dy = (ymax - ymin) / height;
  double x = xmin + dx / 2;
  
  for (int i = 0; i < width; i++)
  {
    double y = ymax - dy / 2;

    for (int j = 0; j < height; j++)
    {
      y = ymax - dy / 2;
      int iterations = countIterations(x, y, maxIterations);
      
      result.Add(new PixelData()
      {  
        X = i,
        Y = j,
        Iterations = iterations < maxIterations ? (int?)iterations : null
      });
      
      y -= dy;
    }
    
    x += dx;
  }
  
  return result;
}


int countIterations(double x, double y, int maxIterations)
{
  int result = 0;
  Complex c = new Complex(x, y);
  Complex z = new Complex(x, y);
  
  while (result <= maxIterations && z.Magnitude() < 2)
  {
    z = z * z + c;
    result++;      
  }

  return result;
}

This simple algorithm can produce some pretty fractals, but it is an iterative method, and this method doesn't translate easily to LINQ.

Eliminating the Outer Loops

The outer loops of the code above are the ones that iterate over the X and Y co-ordinates of the image. Using the Enumerable.Range() iterator we can replace those loops with LINQ statements that have the same effect:

from i in Enumerable.Range(0, width)
   let dx = (xmax - xmin) / width
   let dy = (ymax - ymin) / height
   let x = (xmin + dx / 2) + (i * dx)
   select from j in Enumerable.Range(0, height)
      let y = (ymax - dy / 2) - (j * dy)
      let iterations = countIterations(x, y)
      select new
      {
        X = i,
        Y = j,
        Iterations = iterations < maxIterations ? (int?)iterations : null
      };

This has the same effect as the two nested for-loops in the first example that run over the width and height of the image. But, we're still relying on the countIterations() method, that does the actual calculation of the escape count of the pixel. Clearly, if we are to call this a LINQ fractal generator, we must find some way to do this within in the context of a LINQ query.

This is where it start to get a bit tricky.

Eliminating the inner loop

The countIterations() method essentially counts the number of times it needs to calculate the value of z = z² + c before the magnitude, or absolute value, of z 'escapes' and becomes greater than 2. Since LINQ (or lambda expressions) do not afford us while-loops we must find another way to calculate this.

What we can do is replace the iterative calculation of z with a recursive one. Instead of escaping the algorithm, we can calculate the magnitiude of z for every value between 0 and the maximum number of iterations - then from that set select the maximum value that satisifes our escape condition, that the magnitude of z (| z |) < 2.

We can define a lambda expresssion that calculates the magnitude of z after any iteration by recusively calling itself, like this:

Func<int, Complex, Complex> f = 
    (int n, Complex z) => n > 0 ? f(n - 1, z) * f(n - 1, z) + c : z;

Phew. That's a bit nasty. It also won't compile, as the function f has not yet been fully defined. To get round this, we need to resort to using the Y-combinator. I used the same method as LukeH, and it results in the following code to generate the magnitude of z in the range 0...maxIterations:

from n in Enumerable.Range(0, maxIterations)
  let c = new Complex(x, y)
  let func = (Func<Func<int, Complex, Complex>, Func<int, Complex, Complex>>)
    (f => (n, z) => n > 0 ? f(n - 1, z) * f(n - 1, z) + c : z)
  let escape = Y(func)
  select escape(n, new Complex(x, y))

So, adding a where clause, an orderby clause and using the FirstOrDefault() method, we can get the expression to return the escape iteration count, or 0 if it escapes to infinity, as below.

let c = new Complex(x, y)
let func = (Func<Func<int, Complex, Complex>, Func<int, Complex, Complex>>)
  (f => (n, z) => n > 0 ? f(n - 1, z) * f(n - 1, z) + c : z)
let escape = Y(func)
let iterations = (from n in Enumerable.Range(0, maxIterations)
  where escape(n, new Complex(x, y)).Magnitude() < 2
  orderby n descending
  select n).FirstOrDefault()

So, now all we need to do is fit that into the outer loop query we wrote before.

Final result: the code

The final result looks something like this:

from i in Enumerable.Range(0, width)
  let dx = (xmax - xmin) / width
  let dy = (ymax - ymin) / height
  let x = (xmin + dx / 2) + (i * dx)
  select from j in Enumerable.Range(0, height)
    let y = (ymax - dy / 2) - (j * dy)
    let c = new Complex(x, y)
    let func = (Func<Func<int, Complex, Complex>, Func<int, Complex, Complex>>)
      (f => (n, z) => n > 0 ? f(n - 1, z) * f(n - 1, z) + c : z)
    let escape = Y(func)
    let iterations = (from n in Enumerable.Range(0, maxIterations)
      where escape(n, new Complex(x, y)).Magnitude() < 2
      orderby n descending
      select n).FirstOrDefault()
    select new PixelData()
    {
      X = i,
      Y = j,
      Iterations = iterations < maxIterations ? (int?)iterations : null
    };

Please note the above code is horribly inefficient. Using the iterative version I have been merrily running it over 256 iterations in almost real-time. but the LINQ version takes half a minute or so to execute over 8 - and I haven't managed to have the patience to run it over 16 yet. Since the number of iterations affects the quality of the picture, 8 iterations yields a rather sorry fractal, but given the time this method produces fractals every bit as pretty as the iterative version.