find grid points inside the parallelogram defined by an origin and two vectors

Question

I hope someone knows an efficient computational approach to the following 2D problem:

Given two vectors $\mathbf{A}$ and $\mathbf{B}$, find all grid points that lie within the parallelogram spanned by these vectors.

This feels like it should be a "known" problem, but I lack the vocabulary to search for the right terms. I know how to determine that a given point is inside; I would like not to have to test lots of points... are there any tricks I should know?

As I look at the picture, I am thinking "starting at the origin, you need to move to the left to find points; at a certain X coordinate you can go up a step without crossing A". But as A and B can be swapped, and pointing in any direction, the approach I need has to be a little more robust.

Ultimately I need to know not only the coordinates, but actually the values of $a$ and $b$ for each valid (green, in the diagram) grid point $\mathbf{P}$ such that

$$\mathbf{P} = a\mathbf{A} + b\mathbf{B}$$

for all integer-valued P(i,j) where $a,b \in [0,1\rangle$. If that is actually easier... that's the problem I ultimately need to solve (so I need coordinates (i,j) and their transform (a,b); obviously when I have one, I can find the other).

Answer 1

I suggest that you break the query into two parts, points_inside(parallelogram) := points_inside(triangle1) union points_inside(triangle2), where the triangles are formed by T(origin,origin+A,origin+A+B) andT (origin,origin+A+B,origin+B).

Regarding the find-all-points-in-a-triangle query, this is basically a scanline-conversion problem. You can find solutions in computer-graphics like sources (for drawing triangles onto grids of pixels, that sort of thing).

One of the (many) hits for "scanline conversion of a triangle":

http://vis.uky.edu/~ryang/Teaching/CS335-spr05/Lectures/g_05_fill.pdf

You can probably find better ones. You could probably work with the polygon directly, it'll just be a little more complicated and might be a bit harder to find good tutorials.

You could also try drawing the edges into a small image, then do some floodfill-like search until you hit them. Bresenhams is the prototypical algorithm for scanline conversion of a line.

Answer 2

Depending on the performances that you want, there are an easy and a more complicated algorithm.

Algo 1: interior test in bbox (simple)

For each point p with integer coordinates in the bounding box of the parallelogram (q1,q2,q3,q4)
  If sign(det(p-q1, p-q2)) = sign(det(p-q2, p-q3)) = sign(det(p-q3, p-q4)) = sign(det(p-q4, p-q1))
     mark p

where det( (x1,y1), (x2, y2)) = x1*y2 - x2*y1 denotes the determinant and sign() its sign (positive, zero or negative).

Algo 2: scanline rasterization (more delicate)

Compute the interval [Ymin Ymax] that bounds the Y's of the parallelogram
Determine for each value of Y in [Ymin Ymax] the leftmost XL[Y] and rightmost XR[Y] intersection between the horizontal line and the border of the parallelogram
    For Y in Ymin ... Ymax  
       For X in XL[Y] ... XR[Y]
           Mark X,Y

Step 2 (computing the scalines) can be done using Bresenham algorithm. The algorithm is described in full details in [1]. The second algorithm was used in early 3D game engines (before graphic boards were available) and is very efficient.

More notes (if performance is really an issue)

Starting from Algo 1 and refining it, it is possible to obtain an even faster algorithm: the determinant computations can be split into what depends on X and what depends on Y, and one can compute at the beginning of the loop how the determinant vary. Then in the loop there are only additions and sign comparisons to do. The method is fully detailed in [2] (excellent article by Nick Capens, who implemented a full software emulation of DirectX 10 compatible graphic boards ages ago).

[1] http://www.idav.ucdavis.edu/education/GraphicsNotes/Rasterizing-Polygons/Rasterizing-Polygons.html

[2] http://forum.devmaster.net/t/advanced-rasterization/6145

Answer 3

I'm also unaware of the standard methods of doing this, but I think there's a straightforward way to do it by exploiting linear algebra, assuming both $\mathbf{A}$ and $\mathbf{B}$ can be taken to have their tail at the origin.

If $\mathbf{A} = (A_x, A_y)$, then define $\mathbf{A}_\perp := (A_y, -A_x)$, which is a vector perpendicular to $\mathbf{A}$. The projection of a position vector $\mathbf{P}$ on to $\mathbf{A}_\perp$ gives the distance between $\mathbf{P}$ and the line on which $\mathbf{A}$ lies, up to a sign which indicates whether $\mathbf{P}$ is above or below that line. $$d_{\mathbf{A}}(\mathbf{P}) := \frac{\mathbf{P}\cdot\mathbf{A}_\perp}{\lVert\mathbf{A}_\perp\rVert}$$ The points in your parallelogram are precisely those for which $d_{\mathbf{A}}(\mathbf{P})$ lies between zero and the equivalent value for $\mathbf{B}$, taken as a position vector, and the same is true with the roles of $\mathbf{A}$ and $\mathbf{B}$ reversed: $$\bigl\{\mathbf{P}\ :\ 0 \leq d_{\mathbf{A}}(\mathbf{P}) \leq d_{\mathbf{A}}(\mathbf{B})\ \land\ 0 \leq d_{\mathbf{B}}(\mathbf{P}) \leq d_{\mathbf{B}}(\mathbf{A})\bigr\}$$ (you may have to reverse the signs of one or both inequalities depending on whether $d_{\mathbf{A}}(\mathbf{B})$ and $d_{\mathbf{B}}(\mathbf{A})$ are positive or negative).

If your number of points is not too large, you can devise some enclosing shape - for example a circle with radius $\lVert\mathbf{A}\rVert + \lVert\mathbf{B}\rVert$ (although it shouldn't be hard to do better) - and test all integer grid points within that shape to see which ones satisfy the inequalities. With access to a vectorized math library it should be only a brief piece of code. The computation you wind up doing boils down to computing two dot products and comparing each to the endpoints of a range, for example

(0 <= c1 * px + c2 * py <= c3) and (0 <= c4 * px + c5 * py <= c6)

where the cn's are constants.

Answer 4

I think you're just looking for a straight-forward way of mapping from your original space, call it $(x,y)$ space to your $(a,b)$ space, where $(x,y)=(i\Delta x, j\Delta y)$. Once you're in $(a,b)$ space, you just need to check that $0 \le a \le 1$ and $0 \le b \le 1$.

The reverse is easier to see. Given vectors $\mathbf{A}=(A_x,A_y)$ and $\mathbf{B}=(B_x,B_y)$, some origin $\mathbf{O}=(O_x,O_y)$ where $\mathbf{A}$ and $\mathbf{B}$ start (just to be complete), and scalars $a$ and $b$, a point $\mathbf{P}=(x,y)$ can be written

$\mathbf{P} = a \mathbf{A} + b \mathbf{B} + \mathbf{O}$

or in matrix notation

$\left( \begin{matrix} x \\ y \end{matrix} \right) = \left( \begin{matrix} A_x & B_x\\ A_y & B_y\end{matrix} \right) \left( \begin{matrix} a \\ b\end{matrix} \right) + \left( \begin{matrix} O_x \\ O_y \end{matrix} \right)$

So if you know $(x,y)$ and are looking for $(a,b)$, simply invert this operation

$\left( \begin{matrix} a \\ b\end{matrix} \right) = \left( \begin{matrix} A_x & B_x\\ A_y & B_y\end{matrix} \right)^{-1} \left( \begin{matrix} x-O_x \\ y-O_y \end{matrix} \right) $

Not that this 2x2 matrix inversion is expensive, but it can be performed once ahead of time to save a few electrons. Obviously, this matrix is singular if $\mathbf{A}$ and $\mathbf{B}$ are parallel.