GoogleMap-delbrot

Google Maps, Python and the Mandelbrot Set

There's some surprising uses for a Mandelbrot set (or similar fractal). With just a small bit of code (and a little bit of computation), you can create an image that has detail at any level, and has varied features over the image so its possible to tell where you are. I've used this at work when debugging a scan system that could scan with fields-of-view from 2000nm to 1nm. It also had the ability to scan sections of the image. By writing a Mandelbrot generator as a test scan unit, I could test all of this functionality (I even had the number of iterations tied to the scan pixel time).

I was recently looking at creating a custom Google map, and I thought that using the Google Maps API to browse a Mandelbrot set might be pretty sweet.

It seems other people have had similar ideas, there's a version here and here, although I couldn't find the source for the first one, and the second doesn't include the continuous-zoom feature that Google maps has. This version only runs locally and needs Python for the webserver, and a browser for viewing the final page.

Creating the Mandelbrot Set

I wanted to use Python to create the Mandelbrot set and serve it locally. I split it into two files, one that creates a Mandelbrot set for a given set of coordinates, and another to respond to HTTP requests and return the requested image.

Creating the Mandelbrot set is relatively easy in Python, although a little different than I'd normally do so. Because whole-image operations are optimized in Python, it's faster to use them to create the set rather than iterate over each pixel individually (this also means that we don't escape from the algorithm once a pixel has diverged, we continue to process it for further iterations). I used a slightly modified version of the code from the scipy website, shown below:

#mandelbrot.py
#creates a mandelbrot for the given coordinates
import numpy as np

def mand(iters, tl, br, xpx, ypx=None):
    ypx = ypx or xpx
    c=np.add(*np.ogrid[tl[0]*1j:br[0]*1j:xpx*1j,
                       tl[1]   :br[1]   :ypx*1j])
    m=c.copy()
    ret=np.ones(c.shape, dtype=np.int32)
    for it in xrange(iters):
        m *= m
        m += c
        ret += np.abs(m) < 2
    #want inf. iterations to be a single color, not based on # of iterations
    ret += -(iters+1)*(np.abs(m) < 2)
    return ret

def mand_and_write(file, iters, tl, br, xpx, ypx=None):
    from matplotlib import pyplot
    m=mand(iters, tl, br, xpx, ypx)
    pyplot.imsave(file, m, cmap='spectral', vmin=0, vmax=255)
    
if __name__=="__main__":
    from matplotlib import pyplot
    m=mand(16, (-2, -2), (2, 2), 256)
    pyplot.imshow(m, cmap='spectral') 
    pyplot.show()

If you run this on it's own, it should show a small 256x256 Mandelbrot set:

The next step is to create the webserver, and route requests to our Mandelbrot generator. I decided to just split the request using '/'s, so a GET request with a path of /-2/-2/2/2/16 returns a 256x256 png from (-2,-2) to (2,2) with 16 iterations. The code is below:

#mandserver.py
#simple amnd server that returns pngs for a mandelbrot
from BaseHTTPServer import BaseHTTPRequestHandler, HTTPServer
import mandelbrot
use_threaded_server = True
if use_threaded_server:
    from SocketServer import ThreadingMixIn

import cStringIO


class MandServer(BaseHTTPRequestHandler):
    cache = {}

    def do_GET(self):
        self.send_response(200)
        self.send_header("content-type", "image/png")
        self.end_headers()
        #we'll skip the leading slash
        if self.path not in self.cache:
            coords = self.path[1:].split('/')
            print self.path
            start = -1, -1
            end = 1, 1
            iters = 16
            if len(coords) >= 2:
                start = float(coords[0]), float(coords[1])
            if len(coords) >= 4:
                end = float(coords[2]), float(coords[3])
            if len(coords) >= 5:
                iters = int(coords[4])
            print "requesting ", iters, start, end
            memFile = cStringIO.StringIO()
            mandelbrot.mand_and_write(memFile, iters, start, end, 256)
            self.cache[self.path] = memFile
        else:
            print "Hitting cache!"
            memFile = self.cache[self.path]
        memFile.seek(0)
        self.wfile.write(memFile.read())


def main():
    httpserverclass = HTTPServer
    if use_threaded_server:
        class ThreadedHTTPServer(ThreadingMixIn, HTTPServer):
            pass
        httpserverclass = ThreadedHTTPServer

    try:
        server = httpserverclass(('', 8080), MandServer)
        print 'started MandServer...'
        server.serve_forever()

    except KeyboardInterrupt:
        print '^C received, shutting down server'
        server.socket.close()

if __name__ == '__main__':
    main()

If you run this, with mandelbrot.py from above in the same folder, it should start serving up the PNGs. You can test it (locally) by going to http://localhost:8080/-2/-2/2/2/16 and confirming you see an image similar to the one above (with a slightly different color scheme). The wbeserver runs threaded, so there's a new thread per request. Images are cached as they're created, so if you run this for a long time in can build up a big cache of viewed images.

Now we just need to get Google Maps to request images from this local address in the right format.

Custom Google Map Source

The first step to creating a custom Google Maps source is getting an API Key. You can do that here. Make sure that the 'Google Maps API v3' service is turned on in the 'Services' tab. Then copy your personal API key from the 'API Access' window.

Actually creating the custom map is relatively easy. We need to create a google.maps.ImageMapType that takes a structure containing the URL to request from and some other options, like the maximum zoom size (we use the highest, 19). I decided to automatically increase the iteration count as the zoom increases. The html is below, make sure to replace %Your-API-Key% with your API Key (surprise):

<!DOCTYPE html>
<html>
  <head>
    <meta name="viewport" content="initial-scale=1.0, user-scalable=no" />
    <style type="text/css">
      html { height: 100% }
      body { height: 100%; margin: 0; padding: 0 }
      #map_canvas { height: 100% }
    </style>
    <script type="text/javascript"
      src="http://maps.googleapis.com/maps/api/js?key=%Your-API-Key%&sensor=false">
    </script>
    <script type="text/javascript">

var mandelbrotTypeOptions = {
  getTileUrl: function(coord, zoom) {
      var bound = Math.pow(2, zoom);
 var left=(coord.y/bound - 0.5)
 var top= (coord.x/bound - 0.5)
      return "http://localhost:8080/" +  left +
      "/" + top +
      "/" + (left + 1.0/bound) +
      "/" + (top  + 1.0/bound) +
      "/" + (32*(zoom+1));
  },
  tileSize: new google.maps.Size(256, 256),
  maxZoom: 19,
  minZoom: 0,
  radius: 2,
  name: "Mandelbrot"
};

var mandelbrotMapType = new google.maps.ImageMapType(mandelbrotTypeOptions);

function initialize() {
  var myLatlng = new google.maps.LatLng(0, 0);
  var mapOptions = {
    center: myLatlng,
    zoom: 1,
    streetViewControl: false,
    mapTypeControlOptions: {
      mapTypeIds: ["mandelbrot"]
    }
  };

  var map = new google.maps.Map(document.getElementById("map_canvas"),
      mapOptions);
  map.mapTypes.set('mandelbrot', mandelbrotMapType);
  map.setMapTypeId('mandelbrot');

      }
    </script>
  </head>
  <body onload="initialize()">
    <div id="map_canvas" style="width:100%; height:100%"></div>
  </body>
</html>

And that's it. Just navigate to the webpage and you should see a fully functional Google Maps Mandelbrot set. As mentioned above, the results are cached, so that things should speed up if you look at areas you've previously been to.

You can find all the files above at github.

Hammer Proficiency

Getting Fluent in New Programming Languages

When I first started to learn Python, I wasn't particularly impressed with it. The problem was that searching for an introduction to Python often led to pages that were an introduction to programming using Python as an introductory language, whereas what I really wanted was an introduction to Python for experienced programmers that told me the things that made it different from, say, C. 

So since I couldn't find the introduction I wanted, I referred to the reference docs. The problem with this approach is that it doesn't really get you fluent in a new language, rather it just allows you to write in a language you already know with a different syntax. I'm going to call this 'hammer proficiency', as in if all you have is a hammer, everything looks like a nail. Basically you've got a new tool, and you've learnt enough about it to hit nails with, but not much else.

Here's an example of being 'hammer proficient' in Python. I'm going to write a function that takes a list of strings and returns another list containing only the strings containing the text 'cat', but doing so from a procedural, C-style viewpoint. I know I'm going to return a list, so I'll look up list creation, and I know I'm going to have to iterate over a list, so I'll look up how to iterate over a range of numbers and to find the length of an array. I'm also going to have to look up how to search for strings in a string and add elements to a list. So I would probably come up with something like:

def filter_cat(strings):
    ret = list()
    for i in range(len(strings)):
        if strings[i].find('cat') != -1:
            ret.append(strings[i])
    return ret

There's nothing fundamentally wrong with this. The problem is that any modern language is likely to be flexible enough that you can write code like this, translating almost line-by-line from a more traditional C-style, procedural approach into the new language. But doing so means you don't get to learn the unique and hopefully advantageous features of the program that you're trying to learn.

If I was to write the same function as above now, after having learnt some more Python, I'd probably do the following, using a lambda expression, a list generator and the 'in' operator:

filter_cat2 = lambda strings: [s for s in strings if 'cat' in s]

Which is now using the features included with Python. Only once you start to know not only the syntax, but also the best approaches to take for a given problem for that language, can you really claim to be fluent.

Getting to Fluency

So, how do you get from hammer proficiency to fluency? Experience (and maybe hammer time?) is obviously one way to get there. I didn't find a wealth of information online answering this problem. Like I said, a lot of the introductory information I found started out from too low a level.

StackOverflow has an (unfortunately closed) list of cool features of Python. There's also a really good list of Python tips and tricks here. Finally, I've been doing a few of the Udacity courses (after taking the Stanford 'Introduction to AI' online class). They vary a bit, but the Peter Norvig classes and lectures are all excellent, and his 'Design of Computer Programs' class is too. And because he's using Python, it's serves as an excellent example of how to use the tool, and not like a hammer.

Hilbert Curve in Python and Image Compression

Hilbert Curve in Python

After playing with creating a Hilbert Curve in CFDG, I started thinking about how to implement it in Python.  I wanted a way to do so which was easy to understand and which only needed a small amount of data to describe the process.

I decided to write one rule for one specific generation of the curve, and then rotate that rule 4 times, and then for each of those rules flip them. I decided to use a two-character string to decide the basic shapes. The shape can be one of u, [, n, ] where they represent the orientation of the basic u-shaped curve:

As well as that, I had to include the direction around the shape I was going, either 'a' for anti-clockwise, or 'c' for clockwise. So to draw the shape above starting at the top left, I would call it 'ua': It's in the 'u' orientation, drawing anti-clockwise.

I then had to make a mapping, that for a certain shape, eg 'ua', tells us how to treat that shape when we split it into 4 smaller quadrants. I made the starting shape the key, and the values tuples of (quadrant, newshape) with the quadrants numbered as shown below:

So the mapping for 'ua' to the next generation is:

{'ua': [(0, ']c'), (2, 'ua'), (3, 'ua'), (1, '[c')]}

Now all that's left is to rotate three times, to give the [, n, ], anticlockwise rules, then mirror each of those until we finally have all 8 rules. Once this is done, I wrote a function that takes a square power-of-2 2d array, splits it into quadrants and applys the next rule to each quadrant. I made the function then yield the value of the array at that point, that way it essentially maps a 2d array of values into a list of values taken in order following a hilbert curve. To check it worked I made it us Tkinter to plot the values it was returning. The full code is here (also available on github):

# mfm 18Mar13 Trying to make a generator that returns values
# from a 2d array in a hilbert curve
import numpy


def make_mapping(m):
    """
    Return a dictionary built up of a single mapping m rotated and flipped.
    The keys are 2 character strings where
    the first is one of [,],u,n and the second is either 'a'nticlockwise or
    'c'lockwise. The values are a list of integer, tuple pair which
    specify how to split up the 4 quadrants and which (shape, direction)
    pair to apply to it.
    The quadrants are
    ---------
    | 0 | 1 |
    ---------
    | 2 | 3 |
    ---------
    """
    quads = [0, 1, 3, 2]  # quadrants in clockwise order
    shapes = ['u', '[', 'n', ']']  # shapes rotated clockwise
    rots = ['a', 'c']  # possible rotations

    def next(item, list, pos=1):
        """
        Given a list and an item, reurns the item in the list
        pos places after the item.
        Wraps around the list at the boundaries
        """
        return list[(list.index(item) + pos) % len(list)]

    other_direction = lambda x: next(x, rots)
    next_quadrant = lambda x: next(x, quads)
    next_shape = lambda x: next(x, shapes)

    def rotate(key):
        rotated_value = [(next_quadrant(quad), next_shape(shape) + dirn)
                         for quad, (shape, dirn) in m[key]]
        shape, dirn = key
        return next_shape(shape) + dirn, rotated_value

    def flip(key):
        flipped_value = list(m[key])
        flipped_value.reverse()
        flipped_value = [(quad, shape+other_direction(dirn))
                         for quad, (shape, dirn) in flipped_value]
        shape, dirn = key
        return shape + other_direction(dirn), flipped_value

    key = m.keys()[0]
    while True:
        key, value = rotate(key)
        if key in m:
            break
        m[key] = value

    for key in m.keys():
        newkey, value = flip(key)
        m[newkey] = value

    return m

hilbert_mapping = make_mapping(
    {'ua': [(0, ']c'), (2, 'ua'), (3, 'ua'), (1, '[c')]})


def apply_mapping(ar, mapping, pos):
    """
    split the 2d array ar into 4 quadrants and call
    apply_mapping recursively for each member of map[pos]
    If ar is 1x1, yield ar
    """
    try:
        y, x = ar.shape[:2]
    except ValueError:
        print ar, ar.shape
    if y <= 1 and x <= 1:
        yield ar[0, 0, ...]
    else:
        # quad indices here are:
        #  0  1
        #  2  3
        # (nb unlike hilbert_v1!)
        first_half = slice(None, x / 2)
        secnd_half = slice(x / 2, None)
        quads = [ar[first_half, first_half, ...],
                 ar[first_half, secnd_half, ...],
                 ar[secnd_half, first_half, ...],
                 ar[secnd_half, secnd_half, ...]]
        for quadrant, newpos in mapping[pos]:
            for x in apply_mapping(quads[quadrant], mapping, newpos):
                yield x


def hilbert(ar):
    for x in apply_mapping(ar, hilbert_mapping, 'ua'):
        yield x


if __name__ == '__main__':
    def showpoints(a, side):
        import Tkinter as tk
        from threading import Thread
        scale = 4
        border = 10
        # we draw points in the order given. We assume the value of the array
        # is it's C style position in a contiguous 2d square array
        w = scale * side + 2 * border
        canvas = tk.Canvas(width=w, height=w)
        canvas.pack(expand=tk.YES, fill=tk.BOTH)
        canvas.create_window(w, w)

        def draw():
            startpos = None
            for p in a:
                nextpos = (border + scale * (p % side),
                           border + scale * (p // side))
                if startpos is None:
                    startpos = nextpos
                else:
                    canvas.create_line(startpos[0], startpos[1],
                                       nextpos[0], nextpos[1])
                startpos = nextpos
        Thread(target=draw).start()
        tk.mainloop()

    s = 128
    x = numpy.arange(s * s).reshape(s, s)

    #xl = list(hilbert(x))
    # print xl
    showpoints(hilbert(x), s)

When run, it produces this output:

PNGs and Gifs

The Hilbert curve is meant to give you positions that, generally, are pretty close to one another. So if we take an image, say a disphenocingulum (which just happened to be 512x512 on wikimedia, and I have no idea what a disphenocingulum is ):

And we display is still as a 512x512 image, but this time the pixels follow the Hilbert curve as you go across each row, starting again at the beginning of the next row, then you get:

As you can see, it's definitely made the image, more, erm, horizontal. I was hoping that this could actually lead to improved compression for lossless formats like PNG and GIF. After experimenting with a few images, the 'Hilberted' images were always slightly bigger than the originals for PNG and slightly smaller than the originals for GIF (using PIL/Pillow for making the images).

Here's the size in kB for each method, for a collection of 512x512 PNG files found on the internet. The alpha channel, if there, was removed and set to black for each image beforehand:

Filename	Regular PNG	Hilberted PNG	Regular GIF	Hilberted GIF
002.png	64	89	78	62
android_sh.png	118	143	160	146
Disphenocingulum.png	91	104	165	154
edbang10.png	203	211	91	84
gear.png	107	122	114	95
imaso-512x512.png	253	272	156	140
iPhone 512x512.png	56	55	81	69
quiet 512x512.png	10	21	70	67
Todo-4-App-Icon-512.png	315	342	252	244

I think the reason for this is that GIF strictly encodes left-to-right and top-to-bottom, and in this case the 'Hilberted' versions do have more similarity left-to-right. However, PNGs have the ability to compare either the previous pixel to the left, or the pixel above (see here for details) or a combination of both. This simple approach to treating the data as 2d instead of 1d like GIF seems to really help with images with solid blocks of color, and although it's a lot simpler than the Hilbert method described above, it seems to produce better compression.

For GIFs, the Hilberted version is always slightly smaller (for the images tested). Note that for 5 out of the 9 images tested, the Hilberted GIF provided the smallest filesize out of all of the possibilities.

I'm including the code below. It needs the code above in a file called 'hilbert.py' to work. It has a variety of command line options for outputting processed and unprocessed images, and selecting the format to use. Hopefully they're documented well enough to make some sense.

The code for these two files can also be found on github.

# 2013-04-03 Trying Code && Beer at Substantial.
# what happens if we take the pixels in an image,
# rearrange using a hilbert curve, then save as png?
# is there any increase (or decrease) in compression
try:
    import Image  # needs PIL for this
except ImportError:
    from PIL import Image

import hilbert
import numpy as np
from StringIO import StringIO


def to_hilbert(im):
    im_dat = np.asarray(im)
    w, h = im.size

    if w & (w - 1) != 0 or h & (h - 1) != 0:
        print "Warning, dimensions", w, h, "not powers of two!"

    # let's create a new image for the transformed data
    imout = Image.new(im.mode, (w, h))
    imout_dat = imout.load()

    imout_dat_pos = 0
    for px in hilbert.hilbert(im_dat):
        imout_dat[imout_dat_pos % w, imout_dat_pos // w] = tuple(px)
        imout_dat_pos += 1
    return imout


def from_hilbert(im):
    w, h = im.size

    if w & (w - 1) != 0 or h & (h - 1) != 0:
        print "Warning, dimensions", w, h, "not powers of two!"

    # let's create a new image for the transformed data
    imout = Image.new(im.mode, (w, h))
    imout_dat = imout.load()

    imout_dat_pos = 0
    im_dat = im.load()
    # we create a matrix of x and y positions which we pass in
    xy = np.dstack(np.mgrid[0:h, 0:w])
    for (px, py) in hilbert.hilbert(xy):
        imout_dat[int(py), int(px)] = tuple(
            im_dat[imout_dat_pos % w, imout_dat_pos // w])
        imout_dat_pos += 1
    return imout


def to_image(im):
    w, h = im.size
    # let's create a new image for the transformed data
    imout = Image.new(im.mode, (w, h))
    imout_dat = imout.load()

    imout_dat_pos = 0
    for px in im.getdata():
        imout_dat[imout_dat_pos % w, imout_dat_pos // w] = tuple(px)
        imout_dat_pos += 1
    return imout


def test(file):
    im = Image.open(file)
    imout = to_image(im)
    assert(all(map(lambda (x, y): x == y, zip(im.getdata(), imout.getdata()))))
    imout = to_hilbert(im)
    imout = from_hilbert(imout)
    assert(all(map(lambda (x, y): x == y, zip(im.getdata(), imout.getdata()))))


if __name__ == '__main__':
    # we can test this on the 512x512 pngs in testpngs folder in bash with:
    # for f in testpngs/*; do python hilbert_png.py -t -f GIF -u "$f"; done;
    # this test, using GIF format, with unprocessing too 
    import argparse
    parser = argparse.ArgumentParser()
    parser.add_argument("file", help='The file to load')
    parser.add_argument(
        "--test", "-t", action='store_true',
        help="Doesn't write files, just calculates the sizes")
    parser.add_argument(
        "--reversehilbert", "-r", action='store_true',
        help="Reverses the effect of the program. Will restore the original image")
    parser.add_argument(
        "--unprocessed", "-u", action='store_true',
        help="Write an unprocessed file as well as a processed one.")
    parser.add_argument(
        "--outfile", "-o",
        help="File to write to.", default="out")
    parser.add_argument(
        "--unprocessedoutfile", "-uo",
        help="File to write unprocessed data to.")
    parser.add_argument(
        "--format", "-f",
        help="The format to use.", default="png")
    args = parser.parse_args()
    infile = args.file  # "out.png"
    im = Image.open(infile)
    # test("out.png")
    func = to_hilbert if not args.reversehilbert else from_hilbert
    if args.test:
        s = StringIO()
    else:
        ofile = args.outfile
        ofile += "." + args.format
        s = open(ofile, 'wb')

    imout = func(im)
    imout.save(s, format=args.format)
    if args.test:
        print "Size of hilberted %s %dKb" % (args.format, s.len / 1024)
    s.close()

    if args.unprocessed:
        if args.test:
            s = StringIO()
        else:
            ofile = args.unprocessedoutfile or 'unprocessed'+args.outfile
            ofile += "." + args.format
            s = open(ofile, 'wb')
        
        imout = to_image(im)
        imout.save(s, format=args.format)
        if args.test:
            print "Size of regular %s %dKb" % (args.format, s.len / 1024)
        s.close()

Hilbert curve in CFDG

Hilbert curve in CFDG

In a previous post, I talked about creating Koch snowflakes using the context-free design grammer (CFDG) program. In this post, I'll talk about how a Hilbert curve can also be created using a trick where the basic shape we use draws part of subsequent generations for us.

A Hilbert curve is a little bit trickier (it was an xkcd comic that first introduced me to this fractal) than a Koch snowflake. The first generation looks like:

And the second generation looks like (I've included the first generation in gray on top of it):

The lines in green connect four copies of the first generation shape together. When we go to the third generation we can see four copies of the second generation connected here with red lines (with again the previous generations in gray):

So to create an nth generation Hilbert curve, we create four copies of the previous generation and connect them together with three additional lines. Unlike the Koch snowflake however, there's no single repeated shape, instead we need repeated previous generations plus connecting lines between them.

The upshot of all this is that there's no easy way to draw the nth generation of the Hilbert curve in CFDG (at least in the old version, it may be easier in the latest version) by repeating a simple shape, we need those extra connectors to join them up, and if we're trying to make it in a recursive fashion, there's no easy way to add them. 

If we look at the connectors we need (the green lines in the second generation, and the red lines in the third generation) we can see that they mirror the first generation shape somewhat - the three lines are in the same orientation, but are scaled and rotated. So if it's not possible to draw any generation with a simple shape, is it possible to draw all generations with a simple (although potentially recursive) shape? Can we come up with a shape which draws the connectors for subsequent generations?

The answer is yes. Take a look at the following CFDG code and it's output:

startshape H
 
shape H
{
 loop 2 [flip 90]
  iline[r 90   x -.5]
 iline [y -.5]
}

shape iline
{
 line[]
 iline [ s .5 y .25]
}

shape line
{
    SQUARE[s 1 0.01]
}

Hopefully this shape, when repeated, will give us all generations of the Hilbert curve (with each generation having some lines drawn by the previous generations). If we now repeat this shape four times scaled down by two, we get the following code:

startshape H
 
shape H
{
 loop 2 [flip 90]
  iline[r 90   x -.5]
 iline [y -.5]
 
 loop 2 [flip 90] {
  H [s .5 r 90 x -.5 y  .5]
  H [s .5      x -.5 y -.5]
 }
}
 

shape iline
{
 line[]
 iline [s .5 y .25 ]
}
shape line
{
    SQUARE[s 1 0.01]
}

And when we run this (and wait for 5 million shapes to be drawn) we get:

Which is exactly what we wanted. So here CFDG has allowed us to use a recursive starting shape, applied itself recursively to draw all generations of the Hilbert curve. There is no simple shape we could have used to create a single generation of the curve, and this approach draws all iterations, and must do as we draw the lines for future generations as we go along.

You can find my small collection of CFDG scripts here.