The Development Jungle

Hacking a path through the complexities of software

Rendering the Mandelbrot Set in Clojure

with 2 comments

Clojure Mandelbrot SetClojure is a functional language that brings dynamic types and lisp syntax to the JVM. It has a number of features that make it interesting which will be the subject of another post but for now, let’s dive right into some code.

The Mandelbrot Set

The Mandelbrot set is an example of a fractal which is basically something that can be described using a simple formula and yet exhibit fascinating levels of detail and beauty. I thought it would be a good application to get me started learning Clojure.

The Mandelbrot formula is:

Z = Z2 + C

Where Z is a complex number and C is a constant. To find out whether or not C is in the set, we apply the formula iteratively until the magnitude of Z becomes larger than 2, or we reach a certain number of maximum iterations.

Most 2D fractals are described mathematically using complex numbers which are just like regular (real) numbers but with an added 2nd dimension. They can be represented by two regular numbers so another way to describe the formula above might be:

[x, y] = [x, y]2 + [a, b]

Where [a, b] is equivalent to C and represents the coordinates on the complex plane we’re interested in. By applying the formula we can find out whether or not C is in the Mandelbrot set. If it is, we draw it with a black pixel, otherwise with some other colour depending the number of iterations taken by the formula.

Defining a Functional Mandelbrot Set

Let’s try to write this in Clojure code. First of all, how do you square a complex number? First let’s represent Z as a vector of regular numbers [x, y]. To square them we just apply the formula:

[x, y]2 = [x2 – y2, 2xy]

In Clojure, it looks like this:

(defn mandel [[x y]]
  [(- (* x x) (* y y)) (* 2 x y)])

But this represents just half of the Mandelbrot formula, we need to add that constant C and also somehow call our mandel function iteratively. Normally, in an imperative language, this would be easy. We just assign C to our starting value, then calculate Z in a for loop. In Clojure, it’s a bit different:

(defn mandelformula [xcoord ycoord]
  (iterate
    #(vec (map + (mandel %) [xcoord ycoord]))
    [xcoord ycoord]))

This function uses the iterate function to create successive calls to the mandel function passing our initial starting value [xcoord, ycoord] which represents C in the original formula. The mandel function is actually wrapped in an anonymous function which adds the return value to the initial value hence giving us the Z2 + C. The reason we have to use map is that we cannot directly add the two vectors together. The individual x and y values must be added separately and then constructed back into a vector with the vec function.

That’s quite a dense piece of code but I think represents the mathematical version of the formula fairly well. It’s certainly not the fastest implementation but the goal of functional programming is to favour correctness and readability rather than performance. The great thing about this implementation is that that iterate function returns a lazy infinite sequence. This is very much like the mathematical definition of the formula given at the start; it says nothing about how many times we’re actually going to keep applying the formula.

The next part of the definition of the set says that we need to apply the formula to a given point until either the magnitude of Z becomes greater than 2, or we reach a given number of iterations.

(defn mag [[x y]]
 (+ (* x x) (* y y)))

(take 50 (take-while #(<= (mag %) 4) (mandelformula xcoord ycoord)))

Firstly, we define the function mag which returns the magnitude (actually, the magnitude squared) of the given vector pair. Next we use this function with the take-while function so we’ll only evaluate the terms of the mandelbrot formula while the value of each term is less than or equal to 4 (we’ve squared 2 to avoid having to perform an expensive square root in the mag function). take-while returns another lazy sequence (which might also be infinite) so we use take to restrict the number of terms evaluated to 50.

Now we’re pretty much done with all the maths stuff. All we need to do is count the number of terms calculated and from that number, calculate a colour.

(defn coord-colour
  [[xcoord ycoord] max-iters]
  (let [num-iters (count (take max-iters (take-while #(<= (mag %) 4) (mandelformula xcoord ycoord))))]
    (if (= max-iters num-iters)
    *set-colour*
    (iter-colour num-iters max-iters))))

This function takes a coordinate, a defined number of maximum iterations and returns the colour that the coordinate should be rendered with by applying the Mandelbrot formula. You can see that we’re assigning num-iters to the count of the terms returned by our sequence calculated earlier. We then compare num-iters with max-iters, if they are equal, then the given coordinate is likely to be in the Mandelbrot set so we return *set-colour*. If num-iters is less than max-iters then the coordinate cannot be within the set so we pass the values to another function iter-colour which applies a gradient to the colours.

Calculating a Colour Gradient

First of all we define three vectors representing the colours we want to draw with as RGB values.

;; Colours used to draw the set
(def *grad-colour-a* [255, 255, 0]) ;yellow
(def *grad-colour-b* [0, 0, 255]) ;blue
(def *set-colour* [0, 0, 0]) ;black

We want to calculate a gradient colour such that if num-iters is 0, we return yellow, if it is equal to max-iters we return blue and if it is somewhere in between, we return a mix between the two colours.

(defn grad-colour
  "Returns the colour that is the given fraction of the way between
  the first and second colours given. Returns as a vector of three
  integers between 0 and 255."
  [colA colB frac]
  (vec (map #(+ (* frac (- %2 %1)) %1) colA colB)))

(defn iter-colour
  "Returns the colour needed to paint a point with the given number
  of iterations"
  [num-iters max-iters]
  (grad-colour *grad-colour-a* *grad-colour-b*
    (/ (double num-iters) max-iters)))

The iter-colour function calculates the number of iterations as a fraction of the maximum number of iterations to give a value between 0 and 1 and passes this along with the two colours to the grad-colour function. This function simply calculates the difference between each of the R, G and B values of the two colours and weights it by multiplying by frac. Again, similar to how we calculate the sum of a vector in the mandelformula function, we use map to apply this calculation to each of the three values in the two vectors.

Surfing the Complex Plane

Now we have all we need to calculate the colour of a given coordinate on the complex plane, we’re close to being able to render a fractal in Clojure! But how do we go from the pixels you find on your screen to the mathematical coordinates of the complex plane? Clearly, we need another function.

(defn get-coord
  "Returns the coordinates of the given pixel in the complex plane"
  [x y xstart ystart xsize ysize width height]
  [(+ xstart (* (/ x width) xsize))
  (+ ystart (* (/ y height) ysize))])

This function takes the x and y pixel coordinates and returns the corresponding complex coordinates. The other parameters are the complex coordinate of the top-left of our screen (xstart and ystart), the complex dimensions of the area we’re drawing (xsize and ysize) and the dimensions of the area we are drawing in pixels (width and height).

Now all we need is a way to get all the pixels in our drawing area as vectors to pass on to the this function.

(defn get-pixels [width height]
  (for [y (range height) x (range width)]
    [x y]))

This function returns a lazy sequence of pixel coordinates for the given width and height. We just need to map the get-coord function onto this sequence and we’ll have a lazy sequence of all the complex coordinates in our drawing area. Then all that’s left is to get the colour of each of those complex coordinates (according the mandelbrot formula) and assign the colour to the corresponding pixel.

(defn render [xstart ystart xsize ysize width height max-iters wr]
  (dorun
    (pmap (fn [pixel]
            (let [[x y] pixel]
              (.setPixel wr x y
                (int-array (coord-colour
                  (get-coord (double x) (double y) xstart ystart xsize ysize width height)
                  max-iters)))))
          (get-pixels width height))))

Here we are in the heart of the fractal renderer. Note that we use pmap to process the sequence of pixels which works just like map but runs across multiple threads which allows us to make the most of multi-core system and render each pixel in parallel. Also, we must use dorun because get-pixels returns a lazy sequence, so pmap is forced to actually process each pixel (rather than return a lazy sequence as it would do normally). Finally, the colour returned by each call to coord-colour is converted into a Java int[] array so that it can be passed to the .setPixel function of wr which is an instance of a WritableRaster.

Note also how we coerce the pixel values of x and y into doubles. This is because by default, Clojure will use arbitrary precision using Java’s BigDecimal class which would incur a huge performance hit when it came to applying the formula. Coercing to double ensures we use either double primitives or boxed Double objects.

Setting up the Swing Stuff

All the remaining code is to initialise the Java Swing components needed to display the image generated.

(defn get-img [width height]
  (BufferedImage. width height (BufferedImage/TYPE_INT_RGB)))

(defn get-panel [width height img]
  (proxy [JPanel] [] (paint [g] (.drawImage g img 0 0 (Color/red) nil))))

(defn construct-frame [width height panel]
  "Creates and displays a JFrame of the given dimensions with
  the panel added to it"
  (let [frame (JFrame.)]
    (.setPreferredSize panel (Dimension. width height))
    (doto frame
      (.add panel)
      .pack
      (.setLocationRelativeTo nil)
      .show)))

(defn mandelbrot [xstart ystart xsize ysize width height max-iters]
  "Returns a function to render the mandelbrot set with the
  given parameters on a frame"
  (let [img (get-img width height)
        panel (get-panel width height img)
        wr (.getRaster img)]
    (construct-frame width height panel)
    (fn []
      (do
        (render xstart ystart xsize ysize width height max-iters wr)
        (.repaint panel)))))

We have defined a number of functions for creating the Swing components and also a function for putting them all together called mandelbrot. Notice, rather than call render immediately, the mandelbrot creates and returns an anonymous function instead. This function is actually a closure which means we can call it repeatedly without having to keep passing in the parameters for the window size and number of iterations. When running a program in a REPL, being able to call individual parts of your code can be very handy.

Finally, we assign the function to a var, and then call it!

(def my-mandelbrot (mandelbrot -2 -1.25 3 2.5 1200 1000 50))

(future (my-mandelbrot))

We wrap the call to my-mandelbrot in a future which returns immediately and then executes the function passed to it in a separate thread. The reason for this is that rendering the fractal can take a significant amount of time. Using a future allows you to continue to use the REPL while the rendering happens in the background which can be helpful for debugging or general development. After about 30 seconds or so, you’ll see something like this:

Mandelbrot Set rendered in Clojure

Conclusion

We’ve seen how to create a working fractal rendering program by building on mathematical formulas and making use of Clojure features such as lazy, infinite sequences, parallel processing and closures to do it. These features allow us to keep the code close to the problem domain and hopefully, improve its readability and correctness. In the next post we’ll see how to profile the renderer and drastically improve its performance.

The full source code can be found here on github.

Stay tuned for more Clojure goodness in the future!

Advertisements

Written by alex

November 14, 2009 at 20:08

Posted in programming-languages

Tagged with ,

2 Responses

Subscribe to comments with RSS.

  1. Thanks for this post, it’s awesome. Me and a friend of mine are trying to implement a graphical renderer that you could use with Clojure to draw and display 3D Object. We would like to use it as the bottom layer for a fractal 3D renderer and manipulator, and I think your post will be useful to us.
    Bye!
    Alfredo

    Alfredo Di Napoli

    October 20, 2010 at 17:31

    • That’s great that the post has inspired you to explore the world of fractals and Clojure. 3D fractals are another level complexity and fascination. You might find this page interesting if you haven’t already seen it. There’s tons of useful information and amazing renderings: http://www.skytopia.com/project/fractal/mandelbulb.html

      I would suggest as a start you could try extending the Clojure code I wrote so that it can render any part of the fractal and allow you to zoom into any depth. You will probably find there is still plenty of room for optimizations as well as lots of other interesting challenges in order to do this. If you do start work on any changes it would be great if you did them as a fork of the code I have on github. 🙂

      alex

      October 20, 2010 at 21:14


Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

w

Connecting to %s

%d bloggers like this: