 }
 
 // NOTE
+// https://fgiesen.wordpress.com/2009/12/13/decoding-morton-codes/
 // http://www.forceflow.be/2013/10/07/morton-encodingdecoding-through-bit-interleaving-implementations/
+// http://devblogs.nvidia.com/parallelforall/thinking-parallel-part-iii-tree-construction-gpu/
 
-__forceinline__ __device__ __host__ unsigned int expandBits(unsigned int value)
+__forceinline__ __device__ __host__ uint32_t expandBits(uint32_t value)
 {
 #if 0
 	// NOTE

#else
	return discretized.x << 2 | discretized.y << 1 | discretized.z;
#endif
}

template <typename T>

	return value;
}

#if 0
__forceinline__ __device__ __host__ uint64_t expandBitsBy3(uint32_t value)
{
	uint64_t x = value;
	x = (x | x << 36) & 0x000f000000000ffful;
	x = (x | x << 24) & 0x000f000f000000fful;
	x = (x | x << 12) & 0x000f000f000f000ful;
	x = (x | x << 6) & 0x0303030303030303ul;
	x = (x | x << 3) & 0x1111111111111111ul;
	return x;
}
#else
__forceinline__ __device__ __host__ uint32_t expandBitsBy3(uint32_t value)
{
	uint32_t x = value;
	x = (x | x << 12) & 0x000f000fu;
	x = (x | x << 6) & 0x03030303u;
	x = (x | x << 3) & 0x11111111u;
	return x;
}
#endif	

__forceinline__ __device__ __host__ uint32_t computeMortonCode(aten::vec3 point)
{
	// Discretize the unit cube into a 10 bit integer
	uint3 discretized;
	discretized.x = (uint32_t)min(max(point.x * 1024.0f, 0.0f), 1023.0f);
	discretized.y = (uint32_t)min(max(point.y * 1024.0f, 0.0f), 1023.0f);
	discretized.z = (uint32_t)min(max(point.z * 1024.0f, 0.0f), 1023.0f);

	discretized.x = expandBits(discretized.x);
	discretized.y = expandBits(discretized.y);
