Uniform Grid Ray Casting

VoxelGridIterator.h

/**
 * An iterator for a cubic grid of voxels.
 *
 * The iterator traverses all voxels along a specified direction starting from
 * a given position.
 */
mxBIBREF(
"Ray casting technique described in paper:"
"A Fast Voxel Traversal Algorithm for Ray Tracing - John Amanatides, Andrew Woo [1987]"
"http://www.cse.yorku.ca/~amana/research/grid.pdf"
);
class VoxelGridIterator {
public:
	VoxelGridIterator(
		/// the start position of the ray
		const V3f& _origin,
		/// the direction of the ray
		const V3f& _direction,
		/// the side dimensions of a single cell
		const V3f& _cellSize
	);

	// Post-increment
	VoxelGridIterator operator++();

	// Dereferencing
	Int3 operator*();
	const Int3 & operator*() const;

private:
	/// The integral steps in each direction when we iterate (1, 0 or -1)
	/*const*/ Int3 m_iStep;

	/// The t value which moves us from one voxel to the next
	/*const*/ V3f m_fDeltaT;

	/// The t value which will get us to the next voxel
	V3f m_fNextT;

	/// The integer indices for the current voxel
	Int3 m_iVoxel;
};

/// @todo: vectorize and use branchless DDA?: https://www.shadertoy.com/user/fb39ca4/sort=popular&from=0&num=8
inline
VoxelGridIterator::VoxelGridIterator(
									 /// the start position of the ray
									 const V3f& _origin,
									 /// the direction of the ray
									 const V3f& _direction,
									 /// the side dimensions of a single cell
									 const V3f& _cellSize
)
{
#define ZERO_EPSILON 0.0001f

	mxASSERT(V3_IsNormalized(_direction));

	// 1. Initialization phase.

	const V3f invCellSize = V3_Reciprocal(_cellSize);

	const V3f absDir = V3_Abs(_direction);
	// for determining the fraction of the cell size the ray travels in each step:
	const V3f invAbsDir = {
		(absDir.x > ZERO_EPSILON) ? (1.0f / absDir.x) : 0.0f,
		(absDir.y > ZERO_EPSILON) ? (1.0f / absDir.y) : 0.0f,
		(absDir.z > ZERO_EPSILON) ? (1.0f / absDir.z) : 0.0f
	};

	// 1.1. Identify the voxel containing the ray origin.
	// these are integer coordinates:
	const Int3 startVoxel = {
		Float_Floor( _origin.x * invCellSize.x ),
		Float_Floor( _origin.y * invCellSize.y ),
		Float_Floor( _origin.z * invCellSize.z )
	};

	m_iVoxel = startVoxel;

	// 1.2. Identify which coordinates are incremented/decremented as the ray crosses voxel boundaries.

	/// The integral steps in each primary direction when we iterate (1, 0 or -1)
	m_iStep.x = Sign(_direction.x);	// signum
	m_iStep.y = Sign(_direction.y);
	m_iStep.z = Sign(_direction.z);

	// 1.3. Determine the values of t at which the ray crosses the first voxel boundary along each dimension.

	const float minX = startVoxel.x * _cellSize.x;
	const float maxX = minX + _cellSize.x;
	m_fNextT.x = ((m_iStep.x >= 0) ? (maxX - _origin.x) : (_origin.x - minX)) * invAbsDir.x;

	const float minY = startVoxel.y * _cellSize.y;
	const float maxY = minY + _cellSize.y;
	m_fNextT.y = ((m_iStep.y >= 0) ? (maxY - _origin.y) : (_origin.y - minY)) * invAbsDir.y;

	const float minZ = startVoxel.z * _cellSize.z;
	const float maxZ = minZ + _cellSize.z;
	m_fNextT.z = ((m_iStep.z >= 0) ? (maxZ - _origin.z) : (_origin.z - minZ)) * invAbsDir.z;

	// 1.4. Compute how far we have to travel in the given direction (in units of t)
	// to reach the next voxel boundary (with each movement equal to the size of a cell) along each dimension.

	/// The t value which moves us from one voxel to the next
	m_fDeltaT = V3_Multiply( _cellSize, invAbsDir );

	// 2. Traversal phase is minimal.
	// ...

#undef ZERO_EPSILON
}

inline
VoxelGridIterator VoxelGridIterator::operator++()
{
	// Find the minimum of tMaxX, tMaxY and tMaxZ during each iteration.
	// This minimum will indicate how much we can travel along the ray and still remain in the current voxel.
	const int iMin = Min3Index( m_fNextT.x, m_fNextT.y, m_fNextT.z );

	m_iVoxel[ iMin ] += m_iStep[ iMin ];
	m_fNextT[ iMin ] += m_fDeltaT[ iMin ];

	return *this;
}

inline
Int3 VoxelGridIterator::operator*() {
	return m_iVoxel;
}

inline
const Int3 & VoxelGridIterator::operator*() const {
	return m_iVoxel;
}