I believe I have a solution for you. Suppose the grid is defined by all lines parallel to the $y$-axis through the points $\{x_k = a + k\Delta x\,:\,k \in \mathbb{Z}\}$ and all lines parallel to the $x$-axis through the points $\{y_j = b + j\Delta y\,:\,j\in\mathbb{Z}\}$. Further suppose the line extends from $(x_s,y_s)$ to $(x_e,y_e)$. The following pseudocode outlines the algorithm. At its conclusion the array $P$ contains a list of index pairs $(k,j)$ that index the bottom left corner of the rectangle with vertices $(x_k,y_j)$, $(x_k,y_j + \Delta y)$, $(x_k + \Delta x, y_j)$, $(x_k+\Delta x,y_j + \Delta y)$ that contains a portion of the line.
rtol = macheps
atol = 1e-14
k = floor((xs - a)/dx)
j = floor((ys - b)/dy)
ke = floor((xe - a)/dx)
je = floor((ye - b)/dy)
x = a + k*dx
y = a + j*dy
tcurr = 0
Add the pair (k,j) to an empty array P
while (k does not equal ke or j does not equal je) do
tmin = 2
if (xe does not equal xs) then
tol = rtol*|x| + atol
t = (x + dx + tol - xs)/(xe - xs)
if (tcurr < t and t <= 1) then
tmin = min(tmin, t)
end
t = (x - xs - tol)/(xe - xs)
if (tcurr < t and t <= 1) then
tmin = min(tmin, t)
end
end
if (ye does not equal ys) then
tol = rtol*|y| + atol
t = (y + dy + tol - ys)/(ye - ys)
if (tcurr < t and t <= 1) then
tmin = min(tmin, t)
end
t = (y - ys - tol)/(ye - ys)
if (tcurr < t and t <= 1) then
tmin = min(tmin, t)
end
end
if (tmin is 2) then // should only be possible at the end of the line
Add the pair (ke,je) to the array P
break from the loop
end
tcurr = tmin
k = floor((xs + tmin*(xe - xs) - a)/dx)
j = floor((ys + tmin*(ye - ys) - b)/dy)
Add the pair (k,j) to the array P
x = a + k*dx
y = b + j*dy
end
This algorithm may have to be tweaked depending on how you want to handle edge cases such as when the line passes through grid points.