Analysing space filling curves

gistfile1.py

from math import sin, cos, pi, sqrt, log
import random
import Image
import ImageDraw

algs  = (('gosper', 'x', 6, {
                'x': 'x+yf++yf-fx--fxfx-yf+',
                'y': '-fx+yfyf++yf+fx--fx-y'}),
         ('square', 'x', 4, {
                 'x': 'xf-f+f-xf+f+xf-f+f-x'}),
         ('hilbert', 'l', 4, {
                'l': '+rf-lfl-fr+',
                'r': '-lf+rfr+fl-'}),
         ('hilbert2', 'x', 4, {
                'x': 'xfyfx+f+yfxfy-f-xfyfx',
                'y': 'yfxfy-f-xfyfx+f+yfxfy'}),
         ('peano', 'f', 4, {
                'f': 'f+f-f-f-f+f+f+f-f'}))

def reduce_sentence(sentence, rules, n):
        prev = sentence
        for i in xrange(n):
                cur = []
                for c in prev:
                        if c in rules:
                                cur.extend(rules[c])
                        else:
                                cur.append(c)
                prev = cur
        return cur

def to_points(sentence, angl_frac):
        angl_frac = float(angl_frac)
        point = (0,0)
        angl = 0
        ret = [(0,0)]
        for x in sentence:
                if x == '+':
                        angl += 1
                elif x == '-':
                        angl -= 1
                elif x == 'f':
                        angl = angl % angl_frac
                        point = (point[0] + cos(angl / angl_frac * 2 * pi),
                                 point[1] + sin(angl / angl_frac * 2 * pi))
                        ret.append(point)
        return ret

def normalize_points(points, width, height):
        min_x = min([point[0] for point in points])
        max_x = max([point[0] for point in points])
        min_y = min([point[1] for point in points])
        max_y = max([point[1] for point in points])
        scale_x = width / (max_x - min_x) 
        scale_y = height / (max_y - min_y) 
        return [((point[0] - min_x) * scale_x,
                 (point[1] - min_y) * scale_y) for point in points]

def render(points, width, height):
        points = normalize_points(points, width, height)
        im = Image.new('1', (width, height))
        draw = ImageDraw.Draw(im)
        prev = points[0]
        for point in points[1:]:
                draw.line(prev +  point, fill=128)
                prev = point
        del draw
        im.show()

def d2(v1, v2):
        assert len(v1) == len(v2)
        ret = 0
        for i in xrange(len(v1)):
                ret += (v1[i] - v2[i]) ** 2
        return sqrt(ret)

def render_locality_graph(points):
        im = Image.new('1', (len(points), len(points)))
        pxbuf = im.load()
        data = []
        m = 0
        for i in xrange(len(points)):
                data.append([])
                for j in xrange(len(points)):
                        if points[i] == points[j]:
                                r = 0
                        else:
                                r = abs(i - j) / d2(points[i], points[j])
                        data[-1].append(r+0.1)
                m = max(m, max(data[-1]))
        for i in xrange(len(points)):
                for j in xrange(len(points)):
                        pxbuf[i,j] = int(data[i][j] / m * 255)
        im.show()

def calculate_average_locality(points):
        s = 0
        for i in xrange(len(points)):
                for j in xrange(len(points)):
                        if points[i] == points[j]:
                                r = 0
                        else:
                                r = ((points[i][0] - points[j][0])**2 +
                                   (points[i][1] - points[j][1])**2) / abs(i-j)
                        s += r
        return s / len(points)**2

def unpart1by1(n):
        n&= 0x55555555
        n = (n ^ (n >> 1)) & 0x33333333
        n = (n ^ (n >> 2)) & 0x0f0f0f0f
        n = (n ^ (n >> 4)) & 0x00ff00ff
        n = (n ^ (n >> 8)) & 0x0000ffff
        return n

def deinterleave2(n):
        return unpart1by1(n), unpart1by1(n >> 1)


if __name__ == '__main__':
        #for alg in algs:
        #        for i in xrange(1,8):
        #                sentence = reduce_sentence(alg[1], alg[3], i)
        #                points = to_points(sentence, alg[2])
        #                print alg[0], i, len(points), calculate_average_locality(points)
        #                if len(points) > 300:
        #                        break
        dumb_points = []
        for i in xrange(400):
                dumb_points.append(map(float, deinterleave2(i)))
        print calculate_average_locality(dumb_points)