Many years ago, while doing my thesis, I played a bit with z-order / Morton coding (see https://en.wikipedia.org/wiki/Z-order_curve). Therefore I share here a little snippet of code from ~1995, that I happened to fall upon while cleaning up an old laptop. I have not noted that the equivalent logic has already been shared, therefore I took an hour to write it up in Python (Jupyter notebook style uploaded here). I would be curious to know if it can be found elsewhere.
First a bit of historical context. My first job (before my PhD) had me writing a Design Rule Checker (DRC) for the world's first analog CAD tool, developed at the CSEM https://www.csem.ch/Home (see ILAC https://ieeexplore.ieee.org/document/18603). Among many things, a DRC does boolean operation on layers of integrated circuit masks, and when these operations produce unexpected results, you know that your CAD tool has made a mistake. Focusing on this blog, a DRC must efficiently "join" geometric datasets, and the question is then how should we structure the data of these geometries for optimal performance. At the CSEM I was using K-d trees (https://en.wikipedia.org/wiki/K-d_tree), yet one of the data structures that came out during this period was the Skip list (https://en.wikipedia.org/wiki/Skip_list), and because the skip list is only a one dimensional structure, I thought it would be fun to use z-ordering, to join 2 or more dimensions structures with skip lists.
The difficulty is that keys in the skip list space only have one dimension, and it is not good enough to compare these keys to do high performance geometry operations. We need a way to take these one dimensional keys and optimally project them into the different dimensions of the geometry space. This is what the function "next_above" does below. You give it a boundary of z-space defined by its maximum z-order value, and as well a direction and a reference value, and it will tell you the next z-order value at the edge of the z-space that is aligned with the reference value along the given direction.
This is a 2d z-order example, but the concept works in any dimension.
We need some basic reference values:
log_bit_count = 5
bit_count = (1<<log_bit_count)
mask = 0xFFFFFFFFFFFFFFFF
mask0 = 0x5555555555555555
mask1 = mask0>>1 # not really but avoids python issues
This function finds the next constrained z-order value past a z-order region. Copyright 2020 James Litsios.
In [3]:
def next_above(x, m, z):
a = (x&~m)|(z&m)
b = x^a# different bits
v = b
for kp in range(0, log_bit_count+1):
k = 1<<kp
v |= (v>>k)
w = ((v+1)>>1)&x
if v==0:
w = 1 # case where already on border
nx = (((x|m)+w)&~m)&~v
y = nx | (z&m)
return y
From the web, some utility code to demonstrate things:
In [4]:
def to_morton2(x, y):
# taken from https://stackoverflow.com/questions/30539347/2d-morton-code-encode-decode-64bits
x = (x | (x << 16)) & 0x0000FFFF0000FFFF
x = (x | (x << 8)) & 0x00FF00FF00FF00FF
x = (x | (x << 4)) & 0x0F0F0F0F0F0F0F0F
x = (x | (x << 2)) & 0x3333333333333333
x = (x | (x << 1)) & 0x5555555555555555
x = (y | (y << 16)) & 0x0000FFFF0000FFFF
x = (y | (y << 8)) & 0x00FF00FF00FF00FF
x = (y | (y << 4)) & 0x0F0F0F0F0F0F0F0F
x = (y | (y << 2)) & 0x3333333333333333
x = (y | (y << 1)) & 0x5555555555555555
return x | (y << 1)
def from_morton1(x):
x = x & 0x5555555555555555
x = (x | (x >> 1))& 0x3333333333333333
x = (x | (x >> 2))& 0x0F0F0F0F0F0F0F0F
x = (x | (x >> 4))& 0x00FF00FF00FF00FF
x = (x | (x >> 8))& 0x0000FFFF0000FFFF
x = (x | (x >> 16)) & 0x00000000FFFFFFFF
return x
def from_morton2(d):
x = from_morton1(d)
y = from_morton1(d >> 1)
return x, y
In [5]:
import numpy as np
Choose a size for the example:
In [6]:
size = 64
Let's color a grid with two regions of lower ranges of z-ordered values:
In [7]:
z0, z1, z2 = (0, size*size//10, 2*size*size//3)
grid1 = np.ones((size+1,size+1,3))
for z in range(z0, z1):
x, y = from_morton2(z)
grid1[x,y,:] = (0,0.5,0)
for z in range(z1, z2):
x, y = from_morton2(z)
grid1[x,y,:] = (0.2,0.7,0.2)
x, y = from_morton2(z1)
grid1[x,y,:] = (1.0,0,0)
x, y = from_morton2(z2)
grid1[x,y,:] = (1.0,0,0)
What we have done is to mark in green the x,y values corresponding to all z-order values up to z2, which correspond to the coordinate (0, 63). For clarity, there are two zones, and the last value (0, 63) is shown in red, as is the coordinate (21, 10), that corresponds to z1.
In [8]:
from_morton2(z1)
Out[8]:
(21, 10)
In [9]:
from_morton2(z2)
Out[9]:
(0, 63)
In [10]:
import matplotlib.pyplot as plt
plt.figure(figsize = (5,5))
plt.imshow(grid1)
plt.show()
Now we choose twenty "lookup keys" to demonstrate what "next_above" does:
In [11]:
keys = np.random.randint(0,size,(20,))
In [12]:
keys
Out[12]:
array([34, 49, 10, 12, 45, 61, 10, 59, 16, 31, 37, 23, 12, 36, 27, 28, 26,
First we find the next z-order value outside the two z-order region that aligns "horizontally" with the given keys:
In [13]:
grid2 = grid1.copy()
for key in keys:
grid2[key, :, :] *= 0.9
grid2[key, 0, :] = (0.7,0.7,0.7)
zabove1 = next_above(z1, mask0, to_morton2(key, 0))
xa1, ya1 = from_morton2(zabove1)
grid2[xa1,ya1,:] = (1,0,1)
zabove2 = next_above(z2, mask0, to_morton2(key, 0))
xa2, ya2 = from_morton2(zabove2)
grid2[xa2,ya2,:] = (1,0,1)
In [14]:
plt.figure(figsize = (5,5))
plt.imshow(grid2)
plt.show()
Then, for demonstration purpose, we find the next z-order value outside the two z-order region that aligns "vertically" with the given keys:
In [15]:
grid3 = grid1.copy()
for key in keys:
grid3[:, key, :] *= 0.9
grid3[0, key, :] = (0.7,0.7,0.7)
zabove1 = next_above(z1, mask1, to_morton2(0, key))
xa1, ya1 = from_morton2(zabove1)
grid3[xa1,ya1,:] = (1,0,1)
zabove2 = next_above(z2, mask1, to_morton2(0, key))
xa2, ya2 = from_morton2(zabove2)
grid3[xa2,ya2,:]= (1,0,1)
In [16]:
plt.figure(figsize = (5,5))
plt.imshow(grid3)
plt.show()
The original ~1995 Prolog version of the algorithm:
xor X Y = (X or Y) and not (X and Y);
shiftDown V = shiftDown2 1 V;
shiftDown2 32 V = V;
shiftDown2 K V = shiftDown2 (2*K) (V or (shr V K));
nextAbove X M Z = Y
where
A = (X and not M) or (Z and M),
B = X xor A,
V = shiftDown B,
Y = nextAbove2 X M Z V (ifelse (V=0) 1 ((shr (V+1) 1) and X));
nextAbove2 X M Z V W = Y
where
NX = (((X or M) + W) and not M) and not V,
Y = NX or (Z and M);
C++ version from ~2003:
int
NextAbove(
int x,
int m,
int z){
int a = (x&~m)|(z&m);
int b = x^a; // different bits
int v= b;
{
for(int k=1; k<32; k<<=1){ //smeared bits down
v |= (v>>k);
}
}
int w = ((v+1)>>1)&x; // keep top bit if it is in x
if(v==0){
w = 1; // case where is alreay on "border"
}
int nx = (((x|m)+w)&~m)&~v;
int y = nx | (z&m);
return y;
}