import numpy as np
from itertools import chain

def interleave(*args):
    return int("".join(chain(*zip(*[bin(x)[2:].zfill(x.nbytes*8) for x in args]))), 2)

def morton_order(mesh, dtype=np.uint8):
    vertices = mesh.reshape(mesh.shape[0]*3, 3)
    
    upper_bound = np.max(vertices, axis=0)
    lower_bound = np.min(vertices, axis=0)

    quantize = lambda x: dtype((x-lower_bound)*np.iinfo(dtype).max/(upper_bound-lower_bound))

    zvalue = []
    for triangle in mesh:
        center = np.mean(np.vstack((triangle[0], triangle[1], triangle[2])), axis=0)
        zvalue.append(interleave(*quantize(center)))

    ordered_mesh = np.empty(mesh.shape)

    for i, idx in enumerate(zip(*sorted(zip(zvalue, range(len(zvalue)))))[-1]):
        ordered_mesh[i] = mesh[idx]

    return ordered_mesh