Tuesday, April 9, 2013

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 = [(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:
        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
        y, x = ar.shape[:2]
    except ValueError:
        print ar, ar.shape
    if y <= 1 and x <= 1:
        yield ar[0, 0, ...]
        # 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
                    canvas.create_line(startpos[0], startpos[1],
                                       nextpos[0], nextpos[1])
                startpos = nextpos

    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:

Regular PNG
Hilberted PNG
Regular GIF
Hilberted GIF
iPhone 512x512.png
quiet 512x512.png

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
    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')
        "--test", "-t", action='store_true',
        help="Doesn't write files, just calculates the sizes")
        "--reversehilbert", "-r", action='store_true',
        help="Reverses the effect of the program. Will restore the original image")
        "--unprocessed", "-u", action='store_true',
        help="Write an unprocessed file as well as a processed one.")
        "--outfile", "-o",
        help="File to write to.", default="out")
        "--unprocessedoutfile", "-uo",
        help="File to write unprocessed data to.")
        "--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()
        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)

    if args.unprocessed:
        if args.test:
            s = StringIO()
            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)

1 comment:

  1. Very cool! Just got here from this youtube vid, https://www.youtube.com/watch?v=DuiryHHTrjU and this is a perfect compliment to it. Thanks!