## Introduction

The Hilbert curve is calculated iteratively. Check out this: Hilbert_curve: Representation_as_Lindenmayer_system

## The code

# comment this to enable interactive mode
import matplotlib
matplotlib.use('AGG')

import numpy as np
import matplotlib.pyplot as plt
import os

def apply_rules(s):
""" Hilbert Curve as a Lindenmayer system (L-system)
https://en.wikipedia.org/wiki/Hilbert_curve#Representation_as_Lindenmayer_system"""
s=s.replace("a","-Bf+AfA+fB-")  # capital letters "A" and "B" so that the second operation
return s.lower() # make everyone lowercase

axiom = "a"
n=3 # number of iterations
# displacements, ordered in a counter-clockwise direction
dxdy = np.array([[ 1, 0],    # right
[ 0, 1],    # up
[-1, 0],    # left
[ 0,-1] ])  # down
# displacement is of size 1, so the higher n is, the greater the domain
length = 2**n-1; margin = 0.05*length
domain = [0-margin,length+margin,0-margin,length+margin] # a 5% margin around the curve
s = axiom # string to iterate upon

for i in np.arange(n):
s = apply_rules(s)

make_movie=True
plt.ion() # interactive mode disabled if "matplotlib.use('AGG')"
fig = plt.figure(figsize=(6,6))
ax.axis('off')                        # no frame
ax.axis(domain)                       # domain size
ax.set_aspect('equal')                # square look
ax.set_xticks([]); ax.set_yticks([])  # no ticks
ax.set_title(r"$n = {:d}$".format(n))
plt.show()

# "a" and "b" can be erased now
s=s.replace("a","")
s=s.replace("b","")

frame_names = []  # these two are only relevant if make_movie==True
frame_counter=0

p = np.array([[0.0,0.0]]) # this is the starting point (0,0)
p_plot, = plt.plot(p[:,0],p[:,1],color="black")

# iterate on the string s
for i,c in enumerate(s):
# uncomment to see how fast things are going
# print("{:d}/{:d}".format(i,len(s)))

# rotations "+" and "-" change the displacement array dxdy
# "+" means clockwise rotation
if c == '+': dxdy = np.roll(dxdy,+1,axis=0)
# "-" means counter-clockwise rotation
if c == '-': dxdy = np.roll(dxdy,-1,axis=0)
# forward "f"
if c == 'f':
# add one more point to array p
p = np.vstack([p, [p[-1,0]+dxdy[0,0],p[-1,1]+dxdy[0,1]] ])
# update p_plot data, this is MUCH faster that plotting
# several line segments separately
p_plot.set_data(p[:,0],p[:,1])
fig.canvas.draw()
if make_movie:
fname = "_tmp{:05d}.png".format(frame_counter)
frame_names.append(fname)
fig.savefig(fname,bbox_inches='tight',resolution=300)
frame_counter += 1

if make_movie:
frames = "_tmp%5d.png"
# movie_command = "mencoder mf://*.png -mf fps=24:type=png --ovc lavc -lavcopts vcodec=mpeg4:mbd=2:trell -oac copy -o hil{:d}.avi".format(n)

# we might have other .png figures in the directory
# in this case, use the code below
f = open("png_list.txt", "w")
for i in frame_names:
f.write(i+"\n")
f.close()
movie_command = "mencoder mf://@png_list.txt -mf fps=24:type=png -ovc lavc -lavcopts vcodec=mpeg4:mbd=2:trell -oac copy -o hil{:d}.avi".format(n)

err=os.system(movie_command)
if err!=0:
raise RuntimeError("Couldn't run mencoder.  Data in tmp*.png files")
for fname in frame_names:
os.remove(fname)

# we now have one video ready.
# if you want to join several videos, use this:
# sudo apt-get install gpac
# MP4Box -cat part1.avi -cat part2.avi -new joinedfile.avi