I'm trying to implement this algorithm in C#, but I have very inconsistent results at the moment, and I can't figure out why.
Here is my current implementation :
public static IEnumerable<Point> FastVoxelTraversal(Point a, Point b)
{
var ret = new List<Point>();
Vector direction = ((Vector)b - (Vector)a);
int stepX = direction.x >= 0 ? 1 : -1;
int stepY = direction.y >= 0 ? 1 : -1;
int stepZ = direction.z >= 0 ? 1 : -1;
var iterator = a;
float nextVoxelBoundaryX = (iterator.x + stepX);
float nextVoxelBoundaryY = (iterator.y + stepY);
float nextVoxelBoundaryZ = (iterator.z + stepZ);
float tMaxX = (direction.x != 0) ? (nextVoxelBoundaryX - a.x) / direction.x : float.MaxValue;
float tMaxY = (direction.y != 0) ? (nextVoxelBoundaryY - a.y) / direction.y : float.MaxValue;
float tMaxZ = (direction.z != 0) ? (nextVoxelBoundaryZ - a.z) / direction.z : float.MaxValue;
float tDeltaX = (direction.x != 0) ? 1f / direction.x * stepX : float.MaxValue;
float tDeltaY = (direction.y != 0) ? 1f / direction.y * stepY : float.MaxValue;
float tDeltaZ = (direction.z != 0) ? 1f / direction.z * stepZ : float.MaxValue;
while (iterator != b)
{
if (tMaxX < tMaxY)
{
if (tMaxX < tMaxZ)
{
iterator.x += stepX;
tMaxX += tDeltaX;
}
else
{
iterator.z += stepZ;
tMaxZ += tDeltaZ;
}
}
else
{
if (tMaxY < tMaxZ)
{
iterator.y += stepY;
tMaxY += tDeltaY;
}
else
{
iterator.z += stepZ;
tMaxZ += tDeltaZ;
}
}
ret.Add(iterator);
}
return ret;
}
The cell size is 1.
This is the type of result I get. As you can see, sometimes the results are as expected, and sometimes not.
This does not seem to be related to the coordinates being positive or negative (I have strange results in striclty positive spaces, and good results in negative spaces, and vice versa).

I struggled to find more documentation on this algorithm, so I'm thinking that my variables could not be initialized correctly? What should I correct?
Thank you very much.
EDIT :
Thank you for your answers. I painted the results I'm expecting :

This algorithm should perform a "supercover", meaning that it should include all cells it traverses. This is different from Bresenham, which is supposed to give an approximative rasterization.
nextVoxelBoundaryX = (iterator.x + stepX);looks strange. Initial point may be inside cell, example of initialization