interleave the zvalue bits all at once; big speed improvement in building the bounding volume hierarchy

1 files changed, 10 insertions, 6 deletions

diff --git a/geometry.py b/geometry.py
index a65dbc4..48d19fc 100644
--- a/geometry.py
+++ b/geometry.py
@@ -28,8 +28,14 @@ class Node(object):
 
         self.children = children
 
-def interleave(*args):
-    return int("".join(chain(*zip(*[bin(x)[2:].zfill(x.nbytes*8) for x in args]))), 2)
+def interleave(arr):
+    if len(arr.shape) != 2 or arr.shape[1] != 3:
+        raise Exception('shape mismatch')
+
+    z = np.zeros(arr.shape[0], dtype=np.uint64)
+    for i in range(arr[0][0].nbytes*8):
+        z |= (arr[:,2] & 1 << i) << (2*i) | (arr[:,1] & 1 << i) << (2*i+1) | (arr[:,0] & 1 << i) << (2*i+2)
+    return z
 
 def morton_order(mesh, bits=3):
     lower_bound = np.array([np.min(mesh[:,:,0]), np.min(mesh[:,:,1]), np.min(mesh[:,:,2])])
@@ -40,11 +46,9 @@ def morton_order(mesh, bits=3):
     def quantize(x):
         return np.uint64((x-lower_bound)*max_value/(upper_bound-lower_bound))
 
-    zvalues = np.empty(mesh.shape[0], dtype=np.uint64)
-    for i, triangle in enumerate(mesh):
-        zvalues[i] = interleave(*quantize(np.mean(triangle, axis=0)))
+    mean_positions = quantize(np.mean(mesh, axis=1))
 
-    return zvalues
+    return interleave_arr(mean_positions)
 
 class Solid(object):
     def __init__(self, mesh, inside, outside):