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

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

• I think rchilton1980's suggestion of scanline algorithms in computer graphics is the right direction, but for reference in mathematics this could also be called "enumerating lattice points inside a triangle/polygon" (mathoverflow.net/search?q=enumerating+lattice+points). Lattice points is also easier to google than grid points. – Kirill Oct 6 '17 at 0:46

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.

• "Scanline conversion" sounds like exactly the term I needed... will try to implement this and let you know how it goes. Leaving the question open for a bit longer to see if any other ideas come along. – Floris Oct 6 '17 at 0:45

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

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.

• This is very close to the method I currently have implemented - in essence calculating $a$ and $b$ which are the actual numbers I need. – Floris Oct 6 '17 at 12:07

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.

• The problem is, I need to do this lots of times, and A and B are different each time so there is no value in precomputing (as far as I can see). Efficiency is the key - I am not sure your answer provides it. Did I miss something? – Floris Nov 2 '17 at 4:51
• I thought you had a single set of $\mathbf{A}$ and $\mathbf{B}$, and a lot of grid points. – LedHead Nov 2 '17 at 4:54
• Or rather for every set of $\mathbf{A}$, $\mathbf{B}$ vectors, you have a bunch of grid points to test. This way allows you to do preprocessing once for each $\mathbf{A}$, $\mathbf{B}$ and then all you have to do is a [2x2]*[2xN] matrix-matrix product, where N is the number of candidate grid points. This operation should be very fast. – LedHead Nov 2 '17 at 5:05