# How can you iterate though an arbitrary parametric line in a cartesian grid?

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.

• Are the points 0,1,2,3... supposed to be on the line? Looks like they aren't. Mar 6, 2018 at 20:57
• 0,1,2,3... are not points, they number the boxes (squares in the grid) that you visit and the order you visit them Mar 6, 2018 at 21:10
• Given a point $P=(a,b,c)$ the 6 planes are $x=\lfloor a\rfloor$, $x=\lfloor a\rfloor+1$, and so on for $y$ and $z$. Equivalently, you can subtract $(\lfloor a\rfloor,\lfloor b\rfloor,\lfloor c\rfloor)$ from the point and always work in the cube $[0,1]^3$.
– user856
Mar 7, 2018 at 3:23
• That is not correct, the side of a plane is not necessarily an integer, it's any real so doing the floor won't necessarily get you the correct plane coordinate Mar 7, 2018 at 3:37
• Start by doing a coordinate transformation to the standard grid then. Mar 9, 2018 at 6:55

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).