Look at the following diagram depicting a parametric line $O +\vec vt$:
The goal is simple. Visit every square (squares numbered 0,1,2,3...) exactly and no more than once in the order they are numbered. And each time you move through the line you want to move the origin to the intersection of the dot with the line (*marked with the dots in my diagram).
However, the above is a simplification. The full problem also involves:
The grid is 3D (so we check intersections with planes, not lines). The side of a cube in this grid can be any real number, not necessarily 1.
The cubes are not "centered" at the origin. So the origin of the grid could be different than the origin of the reference frame, and it not be aligned with the distances of the grid at all (i.e the side length may be 0.4 and the origin could be placed at ($\pi$, $\pi$, $\pi$))
The general idea is simple enough. Check intersections with the 3 planes that you can potentially hit, based on the 3 orthogonal components of the direction vector
Check the intersection of the parametric line with each of the 3 planes and project your current point to the intersection to move it.
This sounds easy, except that generating the planes is turning to be a challenge. Given an arbitrary point $P$ determining the 6 planes on the current grid that contain it is being harder than expected and I wanted to know if people had ideas on how to either do this in a different way, or calculate the planes.
The problem becomes much, much simpler, if you invert the situation: find the functions that yield the distance along the line to the $n$'th intersection with the grid plane, for each grid axis.
If the grid is regular, these are linear functions of form $$\begin{aligned} t_{x}(n) &= x_0 + x_t n \\ t_{y}(n) &= y_0 + y_t n \\ t_{z}(n) &= z_0 + z_t n \end{aligned}$$
In fact, if the grid is axis-aligned, with spacing $\Delta_x$, $\Delta_y$, and $\Delta_z$, and your parametric line is $\vec{p}(t) = \vec{o} + t \vec{v}$ where $\vec{o} = ( o_x , o_y , o_z )$ and $\vec{v} = ( v_x , v_y , v_z )$, then $$x_t = \frac{\Delta_x}{v_x}, \quad y_t = \frac{\Delta_y}{v_y}, \quad z_t = \frac{\Delta_z}{v_z}$$ and if the grid origin is at $\vec{g} = ( g_x , g_y , g_z )$, and we define $f(x) = x - \lfloor x \rfloor$, then $$\begin{aligned} x_0 &= \frac{\Delta_x}{v_x} f\left( v_x \frac{o_x - g_x}{\Delta_x} \right ) \\ y_0 &= \frac{\Delta_y}{v_y} f\left( v_y \frac{o_y - g_y}{\Delta_y} \right ) \\ z_0 &= \frac{\Delta_z}{v_z} f\left( v_z \frac{o_z - g_z}{\Delta_z} \right ) \end{aligned}$$
Even in the general case (arbitrarily-oriented regular grid), you only need to find the two first intersections of each grid plane from $\vec{o}$ along $\vec{v}$. The distances to the first intersections yield $x_0$, $y_0$, and $z_0$; and the distances between the first and second intersections yield $x_t$, $y_t$, and $z_t$, assuming you measure distance in units of $t$ ($\lVert\vec{v}\rVert$).
You iterate starting with $n_x = n_y = n_z = 0$:
If $t_x(n_x) = t_y(n_y) = t_z(n_z)$, intersection occurs at a grid cell corner at $t = t_x(n_x)$.
Increment $n_x$, $n_y$, and $n_z$ for next iteration.
If $t_x(n_x) = t_y(n_y) \lt t_z(n_z)$, intersection occurs at a grid cell edge at $t = t_x(n_x)$.
Increment $n_x$ and $n_y$ for next iteration.
If $t_x(n_x) = t_z(n_z) \lt t_y(n_y)$, intersection occurs at a grid cell edge at $t = t_x(n_x)$.
Increment $n_x$ and $n_z$ for next iteration.
If $t_y(n_y) = t_z(n_z) \lt t_x(n_x)$, intersection occurs at a grid cell edge at $t = t_y(n_y)$.
Increment $n_y$ and $n_z$ for next iteration.
If $t_x(n_x) \lt t_y(n_y)$ and $t_x(n_x) \lt t_z(n_z)$, intersection occurs at a grid cell face at $t = t_x(n_x)$.
Increment $n_x$ for next iteration.
If $t_y(n_y) \lt t_x(n_x)$ and $t_y(n_y) \lt t_z(n_z)$, intersection occurs at a grid cell face at $t = t_y(n_y)$.
Increment $n_y$ for next iteration.
If $t_z(n_z) \lt t_x(n_x)$ and $t_z(n_z) \lt t_y(n_y)$, intersection occurs at a grid cell face at $t = t_z(n_z)$.
Increment $n_z$ for next intersection.
As it happens, a very similar approach can be used to rasterize curves, including cubic Bézier curves, at an arbitrary (but fixed) precision. I find it an interesting future approach to implement true (limited by machine precision only) cubic curves with full velocity control for Cartesian robots (3D printers et cetera).