Fractals Graphics Math Programming Python

Visualizing Data Using the Hilbert Curve

Some time ago, a coworker asked me to help him visualize some data. He had a very long series (many millions) of data points, and he thought that plotting a pixel for each one would visualize it well, so he asked for my help.
I installed Python & PIL on his machine, and not too long after, he had the image plotted. The script looked something like:

data_points = get_data_points()
n =  int((len(data_points)**0.5) + 0.5)
image = Image('1', (n, n))
for idx, pt in enumerate(data_points):
	image.putpixel(pt, (idx/n, idx%n))'bla.png', 'png')

Easy enough to do. Well, easy enough if you have enough memory to handle very large data sets. Luckily enough, we had just enough memory for this script & data series, and we were happy. The image was generated, and everything worked fine.

Still, we wanted to improve on that. One problem with this visualization is that two horizontally adjacent pixels don’t have anything to do with each other. Remembering xkcd’s “Map of the Internet“, I decided to use the Hilbert Curve. I started with wikipedia’s version of the code for the Python turtle and changed it to generate a string of instructions of where to put pixels. On the way I improved the time complexity by changing it to have only two recursion calls instead of four. (It can probably be taken down to one by the way, I leave that as a challenge to the reader :)

Unfortunately, at this point we didn’t have enough memory to hold all of those instructions, so I changed it into a generator. Now it was too slow. I cached the lower levels of the recursion, and now it worked in reasonable time (about 3-5 minutes), with reasonable memory requirements (no OutOfMemory exceptions). Of course, I’m skipping a bit of exceptions and debugging along the way. Still, it was relatively straightforward.

Writing the generator wasn’t enough – there were still pixels to draw! It took a few more minutes to write a simple “turtle”, that walks the generated hilbert curve.
Now, we were ready:

hilbert = Hilbert(int(math.log(len(data_points), 4) + 0.5))
for pt in data_points:
    x,y =
    image.putpixel(pt, (x,y))

A few minutes later, the image was generated.

Graphics Programming Python

Python + Curses + Numpy -> Ascii Art Rotating Cube

Well, the title says it all, or at least most of it.
If you are still not convinced, a picture’s worth 500 dwords:

        / |\
       /  | \
      /   |  \
    //    |   \\
   /      |     \
  /       \      \
 |       / \      |
 |\     /   \    /|
 | \   /     \\ / |
 |  \ /       //  |
 |  /\       /  \ |
 | /  \\    /    \|
 |/     \  /      |
 /       \/       /
  \       |      /
   \      |     /
    \\    |   //
      \   |  /
       \  | /
        \  /

I’ve been toying with the idea of writing a roguelike, and I thought that it would be really fun if it had an ‘Eye-Of-The-Beholder’ kind of view as well, just in ascii-art. There are plenty of things I wanted to do with ascii-art, since it’s such a nice toy. So I started writing some code, and this is a small test on the way.

Nowadays, when I discuss 3d ascii-art, people usually think about a picture transformed with aalib or something similar. Although it is a (very) nice idea, I think it would be more fun to play with ascii-art directly.

Other playing possibilities include a ray-tracer. You could write some fps with ray-tracing ascii-art graphics, because there are so few pixels.
Maybe I’ll write one, when I get to it.

Oh! The code, what about the code?
Here it is.
A few notes:
1. The code ain’t too pretty, as it is a work in progress that doesn’t get much time invested in it. However, I’d still rather publish it.
2. It uses curses. Under windows you’ll have to get some replacement module.
3. It uses numpy. If you don’t have it already, be sure to get it.
4. Run ‘’. Use the direction keys to rotate the cube, ‘q’ to exit.

Fractals Graphics Programming Python

A second look at the dragon fractal

At first when I drew the (twin)dragon fractal, I had a small bug. I used the base 1+i instead of 1-i. This also generated a very similar looking beast. Thinking about that for a while, made me curious. Just like the Mandelbrot set, maybe other interesting fractals could be generated using exactly the same method, with a different complex ‘seed’.

So I patched up my code, and made it output a series of images for complex numbers on a path. I thought the reasonable choice would be a spiral, so using t*(cost+isint) as a base I wrote a spiral that would go around the origin several times and finally land on 1-i.
It might seem obvious to you to try and make this interactive instead – why not move the mouse around and watch the different fractals that emerge? Well, I wanted a video.

I also wanted the fractal to convey a little more information. So each point in the set was given a color according to its generation. I decided after some trial and error that white was better for the older generation, and red for younger generations.
After some PIL and ffmpeg work, it was ready.

While working on the video, I witnessed some very interesting fractals with different seeds.
I was very curious as to why certain shapes emerged. For example, a square pattern was very common:
Squares fractal
This turned out to be not that surprising. Its generator number is of the shape 0+ia. So it made sense. I still didn’t figure out how come hexagon shapes were so common:
Hexagons1 Hexagons2 Hexagons1

Maybe it had something to do with pi/6, I’m not too sure about that. If it did, I would expect to see many more regular polygons.
Here is another curious one:

Another interesting phenomenon was what I started to call ‘the dragon’s heart’. You see, in the final dragon, the starting generations were spread about pretty evenly. However, with other bases, even ones generating something which is pretty similar to the dragon, the starting generations are clustered together – sometimes at the side, sometimes at the middle.

I’ve got a feeling (not proven or anything) that to be space filling, a fractal’s starting generations should be spread about. What do you think?

Fractals Graphics Programming Python

Fractals in 10 minutes No. 2 – Recursive Spring

Imagine a straight wire. Bend this wire until its ends meet. You get a ring. Next stage.
Take another straight wire, bend it as before, but this time, don’t let the ends meet, instead continue bending the wire more and more, until you get a spring. Now, Think of this spring as the first wire, and bend it until the spring’s ends meet. We get a spring-ring. (See pictures below)

At this point it should be obvious – instead of making the ends of the spring meet, we can continue bending the spring into a meta-spring, and again make this meta-spring’s ends meet. Rinse and repeat ad-infinitum.

At first I had a hard time picturing it, and I thought it would be tough to write about it without drawing it. So I set with pen and paper, and tried to draw it. That was even harder than picturing it. So I thought I could write a script to draw it, and that turned out to be pretty easy.

Here are the first four stages. The second and third are shown from the side. Click for bigger versions.

To implement this, I wrote a function that draws any track. A track is defined to be a function takes a float argument between 0 and 1, and returns the appropriate 3d point in between, where 0 is one end of the track and 1 is the other. (Actually, I also decided that tracks should return a “right hand” vector as well). Then I wrote a ring:

def ring(t):
    pi = math.pi
    cos = math.cos
    sin = math.sin
    x = numpy.array((2*cos(2*pi*t),2*sin(2*pi*t),0))
    return (x,-x/3)

Lastly I wrote a wrapper function, that takes any track, and springifies it. It was called springify_track. At this point I could do this:
(The constants are the radius of the spring loops, the sampling delta, and the number loops.)

t1 = ring
t2 = springify_track(t1, 0.25, 0.1, 50)
t3 = springify_track(t2, 0.06, 0.0011, 5000)
t4 = springify_track(t3, 0.011, 0.000012, 500000)
Algorithms Geometry Graphics Math Programming Python

Elegant Line Clipping Algorithm

I actually wrote about that a long time ago, in the wiki that once was Since the wiki is long gone, and I figured this is a very elegant algorithm, I decided I should write about it.

First, some background – while writing Ascii-Plotter, I had to plot lines. To plot lines correctly I had to clip them. Now since this was a fun project, I decided I wanted to figure out how to do that myself, even though the problem is really solved. So I sat down with pen and paper, and did a little thinking, and I came up this nice algorithm.

The idea is simply this: consider a line, that is not parallel to any of the axes. This line intersects all the four lines that make up the limits of our screen. Here is a diagram:

An intersection of a line with a bounding box

Now notice that there are six points of interest: the four intersection points, and the two endpoints of the line. For each point, we find its parameter in the parametric representation of the line, where a and b are the endpoints:

x(t) = a+t(b-a)

To find out the actual endpoints of the clipped line, we sort the six points according to their parameter, and take the two innermost points (with zero based indexing, those would be 2 and 3). What if the line is outside the box you say? We check that the result points’ parameters are between 0 and 1 (which correspond to the original endpoints).

Note that for a line parallel to one of the axes, a similar method applies, but without two of the intersection points. However, it is probably easier to just calculate the clipping directly in that case.

Here is some sample code, without all the surrounding checks (for parallelism with the axes, trivial clipping etc…).
The full version is available in the code of Ascii-Plotter.

points_t = [0.0,1.0]
points_t.append( float(rect_bottom_left[0]-line_pt_1[0])/(line_pt_2[0]-line_pt_1[0]) )
points_t.append( float(rect_top_right[0]-line_pt_1[0])/(line_pt_2[0]-line_pt_1[0]) )
points_t.append( float(rect_bottom_left[1]-line_pt_1[1])/(line_pt_2[1]-line_pt_1[1]) )
points_t.append( float(rect_top_right[1]-line_pt_1[1])/(line_pt_2[1]-line_pt_1[1]) )
if points_t[2] < 0 or points_t[2] >= 1 or points_t[3] < 0 or points_t[2]>= 1:
	return None
result = [(pt_1 + t*(pt_2-pt_1)) for t in (points_t[2],points_t[3]) for (pt_1, pt_2) in zip(line_pt_1, line_pt_2)]
return (result[0],result[1]), (result[2], result[3])