morton.py

%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib.path import Path
import matplotlib.patches as patches
width = 9
height = 3

def next_greater_power_of_2(x):  
    return 2**(x-1).bit_length()

r_width = next_greater_power_of_2(width)
r_height = next_greater_power_of_2(height)
print width, r_width, r_width.bit_length()
print height, r_height, r_height.bit_length()


print "{0}x{1} rounded to {2}x{3}".format(width, height, r_width, r_height)

def morton_encode2d(x, y, width, height):
    r_width = next_greater_power_of_2(width)
    r_height = next_greater_power_of_2(height)
    min_axis = min(r_width.bit_length(), r_height.bit_length())
    
    d = 0
    pos = 0
    for i in xrange(min_axis-1):

        a = x & 0x1
        b = y & 0x1
        x = x >> 1
        y = y >> 1
        temp = ((b<<1)|a)
        d |= temp << pos
        pos = pos + 2

    d |= (x|y) << pos

    return d

x = 6
y = 2
width = 9
height = 3
d = morton_encode2d(x, y, width, height)

def morton_decode2d(d, width, height):
    r_width = next_greater_power_of_2(width)
    r_height = next_greater_power_of_2(height)
    min_axis = min(r_width.bit_length(), r_height.bit_length())

    x = 0
    y = 0
    pos = 0

    for i in xrange(min_axis-1):
        x |= (d&0x1) << pos
        d = d>>1
        y |= (d&0x1) << pos
        d = d>>1
        pos = pos +1

    if width > height:
        x |= d << pos
    if height > width:
        y |= d << pos
    
    return(x,y)

(a,b) = morton_decode2d(d, width, height)


def draw_morton_order(width, height):
    r_width = next_greater_power_of_2(width)
    r_height = next_greater_power_of_2(height)
    verts = []
    codes = []
    
    d = 0

    for x in xrange(r_width):
        for y in xrange(r_height):

            (a, b) = morton_decode2d(d, width, height)
            d = d+1
            
            if (a < width) and (b < height):
                # add to path
                verts.append((a,b))
#                 print "({0},{1}) -> ({2},{3})".format(x,y,a,b)
                codes.append(Path.LINETO)
#             else:
#                 print "({0},{1}) -> skip".format(x,y)
                
    codes[0]  = Path.MOVETO
    print len(verts)
    path = Path(verts, codes)

    fig = plt.figure()
    ax = fig.add_subplot(111)
    patch = patches.PathPatch(path, facecolor='orange', lw=2)
#     ax.add_patch(patch)
    ax.set_xlim(-1,width)
    ax.set_ylim(-1,height)
    ax.grid()
    xx, yy = zip(*path.vertices)
    line, = ax.plot(xx, yy, 'go-')
    plt.gca().invert_yaxis()
    plt.show()

width = 6
height = 3
draw_morton_order(width, height)



for y in xrange(height):
    for x in xrange(width):
        print x,y, morton_encode2d(x, y, width, height)