algorithm - How can you iterate linearly through a 3D grid?

Question

Assume we have a 3D grid that spans some 3D space. This grid is made out of cubes, the cubes need not have integer length, they can have any possible floating point length.

Our goal is, given a point and a direction, to check linearly each cube in our path once and exactly once.

So if this was just a regular 3D array and the direction is say in the X direction, starting at position (1,2,0) the algorithm would be:

for(i in number of cubes)
{
     grid[1+i][2][0]
}

But of course the origin and the direction are arbitrary and floating point numbers, so it's not as easy as iterating through only one dimension of a 3D array. And the fact the side lengths of the cubes are also arbitrary floats makes it slightly harder as well.

Answer 1

Assume that your cube side lengths are s = (sx, sy, sz), your ray direction is d = (dx, dy, dz), and your starting point is p = (px, py, pz). Then, the ray that you want to traverse is r(t) = p + t * d, where t is an arbitrary positive number.

Let's focus on a single dimension. If you are currently at the lower boundary of a cube, then the step length dt that you need to make on your ray in order to get to the upper boundary of the cube is: dt = s / d. And we can calculate this step length for each of the three dimensions, i.e. dt is also a 3D vector.

Now, the idea is as follows: Find the cell where the ray's starting point lies in and find the parameter values t where the first intersection with the grid occurs per dimension. Then, you can incrementally find the parameter values where you switch from one cube to the next for each dimension. Sort the changes by the respective t value and just iterate.

Some more details:

cell = floor(p - gridLowerBound) / s   <-- the / is component-wise division

I will only cover the case where the direction is positive. There are some minor changes if you go in the negative direction but I am sure that you can do these.

Find the first intersections per dimension (nextIntersection is a 3D vector):

nextIntersection = ((cell + (1, 1, 1)) * s - p) / d

And calculate the step length:

dt = s / d

Now, just iterate:

if(nextIntersection.x < nextIntersection.y && nextIntersection.x < nextIntersection.z)
    cell.x++
    nextIntersection.x += dt.x
else if(nextIntersection.y < nextIntersection.z)
    cell.y++
    nextIntersection.y += dt.y
else
    cell.z++
    nextIntersection.z += dt.z
end if
if cell is outside of grid
    terminate

I have omitted the case where two or three cells are changed at the same time. The above code will only change one at a time. If you need this, feel free to adapt the code accordingly.