'(Mandelbrot Set Attractors)

As you have no doubt already seen, the Mandelbrot Set is one of the most popular mathematical functions to calculate with a computer. If you have any doubts, just take a look at this Google search of Mandelbort Set images. Although there is plenty of room for creating stunning images, most people tend to calculate the set the same basic way.

On this page, I want to demonstrate an alternative view of this remarkable set. I iterate the same basic function in more or less the same way everyone else does. The difference is, I am looking for the behavior of Z. For a long time, I've wondered what it looks like when Z never escapes to infinity. Here I start to find out.

I started with the standard map. The real axis goes from -2.0 to 2.0 going left to right. The imaginary axis goes from -2.0 to 2.0 going top to bottom. Yes, this is upside down. This doesn't matter because the M-set is symmetrical about the real axis. The array that represents the pixels is initialized to zero. Z is initialized to 0.0 and C gets the value of the pixel coordinates. Top left is (-2.0, -2.0). The next value of Z is Z squared plus C. Instead of coloring each pixel according to when Z exceeds a certain value, I increment the value represented by each pixel that Z lands in. The more times Z is in a pixel, the brighter it gets.

After determining that this procedure does indeed generate something of interest, I produced the following image.

Z has a tendency to hit the same pixels over and over. This makes the popular pixels quite bright compared to the others. I am using eight bit grey scale to draw the image, so I want to map the brightest pixel to 255 without letting the others be too dark. To do this, I used a logarithmic mapping. This worked fairly well as you can see, but the image is "grainy". To get around this, I added sub-pixel anti-aliasing to my algorithm. This is a lot more work as I am now computing nine sub-pixels for each drawn pixel. The result though is much smoother.

As you can see, the anti-aliasing worked out fairly well. I don't know where those fine lines are comming from though. The next thing I decided to try was just tracking the Z values that never escape. After all, I'm interested mainly in the orbits that are trapped in the set. To do this, I added some simple but expensive code to keep a list of undo information for the Z values. If Z escaped, I would go back and undo its increments. This leaves me with just the Z values that don't escape. It makes quite a difference to the picture. There is more detail, but the artistic merits are debatable.

The code that produces these images is very simple. I've made the raw source available if you wish to play with it. The listing as it now stands is below.

;;;; mset-attractors.lisp -- finding the points of attraction in the mandelbrot set
;;;; 
;;;; Copyright (c) 2005
;;;;
;;;; By David Steuber <david@david-steuber.com>

(defun make-axis-array (n)
  (let ((aa (make-array n :element-type 'double-float))
        (m (1- n)))
    (dotimes (i n)
      (setf (aref aa i)
            (+ -2.0d0 (* 4.0d0 (/ i m)))))
    aa))

(defun make-attractor-map (resolution max-iterations)
  (declare (type fixnum resolution max-iterations)
           (optimize speed))
  (let* ((r3 (* 3 resolution))
         (a3 (make-axis-array r3))
         (aa (make-axis-array resolution))
         (am (make-array (list resolution resolution) :element-type 'fixnum :initial-element 0))
         (interval (- (aref aa 1) (aref aa 0)))
         (min-val (aref aa 0)))
    (declare (type double-float interval min-val)
             (type fixnum r3))
    (labels ((index-from-float (float)
               (declare (type double-float float))
               (when (<= -2.0d0 float 2.0d0)
                 (round (- float min-val) interval)))
             (escape-fn (cr ci)
               (declare (type double-float cr ci))
               (loop
                  for iterations of-type fixnum below max-iterations
                  for z-imag of-type double-float         = 0.0D0 then (+ (* 2.0d0 z-real z-imag) ci)       
                  for z-real of-type double-float         = 0.0D0 then (+ (- z-real-squared z-imag-squared) cr)
                  for z-real-squared of-type double-float = 0.0D0 then (* z-real z-real)
                  for z-imag-squared of-type double-float = 0.0D0 then (* z-imag z-imag)
                  with points = ()
                  while (< (+ z-real-squared z-imag-squared) 4.0D0) do
                    (when (< 0 iterations)
                      (let ((x (index-from-float z-real))
                            (y (index-from-float z-imag)))
                        (and x y (incf (aref am x y))
                             (push (cons x y) points))))
                  finally (when (< (length points) (1- max-iterations))
                            (mapc (lambda (p)
                                    (decf (aref am (car p) (cdr p))))
                                  points)))))
      (dotimes (y r3)
        (dotimes (x r3)
          (escape-fn (aref a3 x) (aref a3 y))))
      am)))

(defun write-pgm-file-from-map (map filespec)
  ;; http://astronomy.swin.edu.au/~pbourke/dataformats/ppm/
  (let* ((width (array-dimension map 1))
         (height (array-dimension map 0))
         (max-val (loop for row below height
                     maximize (loop for col below width
                                 maximize (aref map row col))))
         (b (expt (coerce max-val 'double-float) #.(/ 1.0d0 255.0d0))))
    (with-open-file (pgm filespec :direction :output :if-exists :supersede)
      (format pgm "P2 ~A ~A ~A~%" width height 255)
      (dotimes (y height)
        (dotimes (x width)
          (format pgm "~A~%" (let ((n (coerce (aref map x y) 'double-float)))
                               (if (zerop n)
                                   0
                                   (round (log n b))))))))))

Melinda Green has done something similar with The Buddhabrot Technique. The main difference is she is dealing with the points that escape while I was looking at the ones that don't.