 }
 
 // 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

	return value;
}

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

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

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

__forceinline__ __device__ __host__ uint32_t computeMortonCode(float x, float y, float z, float s)
{
	// Discretize the unit cube into a 8 bit integer
	uint32_t dx = (uint32_t)min(max(x * 256.0f, 0.0f), 255.0f);
	uint32_t dy = (uint32_t)min(max(y * 256.0f, 0.0f), 255.0f);
	uint32_t dz = (uint32_t)min(max(z * 256.0f, 0.0f), 255.0f);
	uint32_t dw = (uint32_t)min(max(s * 256.0f, 0.0f), 255.0f);

	uint32_t ddx = expandBitsBy3(dx);
	uint32_t ddy = expandBitsBy3(dy);
	uint32_t ddz = expandBitsBy3(dz);
	uint32_t ddw = expandBitsBy3(dw);

	uint32_t ret = (uint32_t)(ddx << 3 | ddy << 2 | ddz << 1);
	//uint32_t ret = (uint32_t)(ddx << 2 | ddy << 2 | ddz << 1);

	return ret;
}

__forceinline__ __device__ __host__ uint32_t onComputeMortonCode(
	const aten::vec3& vmin,
	const aten::vec3& vmax,
	const aten::aabb& sceneBbox)
{
	// Normalize [0, 1].
	const auto size = sceneBbox.size();
	const auto bboxMin = sceneBbox.minPos();
	
	aten::vec3 center = (vmin + vmax) * 0.5f;
	center = (center - bboxMin) / size;

	auto code = computeMortonCode(center);

	return code;
}

__forceinline__ __device__ __host__ uint32_t onComputeExtendedMortonCode(
	const aten::vec3& vmin,
	const aten::vec3& vmax,
	const aten::aabb& sceneBbox)
{
	// Normalize [0, 1].
	const auto size = sceneBbox.size();
	const auto bboxMin = sceneBbox.minPos();

	aten::vec3 center = (vmin + vmax) * 0.5f;
	center = (center - bboxMin) / size;

	// Compute diagonal.
	//const auto div = sceneBbox.getDiagonalLenght();
	//auto d = aten::length(vmax - vmin);
	//d /= div;

	auto code = computeMortonCode(center.x, center.y, center.z, 0);

	return code;
}

template <typename T>
