## 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.
---------
| 0 | 1 |
---------
| 2 | 3 |
---------
"""
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_shape = lambda x: next(x, shapes)

def rotate(key):
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()
for quad, (shape, dirn) in flipped_value]
shape, dirn = key
return shape + other_direction(dirn), flipped_value

key = m.keys()
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:
#  0  1
#  2  3
# (nb unlike hilbert_v1!)
first_half = slice(None, x / 2)
secnd_half = slice(x / 2, None)
ar[first_half, secnd_half, ...],
ar[secnd_half, first_half, ...],
ar[secnd_half, secnd_half, ...]]
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
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, startpos,
nextpos, nextpos)
startpos = nextpos
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_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_pos = 0
# 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_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()
"--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()
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()

```

#### 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!