Covering a 2D plotting area with lattice points

My goal is to cover the plotting area with lattice points.

In this example we are working in 2D. We call the set Λ ⊆ R^2 a lattice if there exists a basis B ⊆ R^2 with Λ = Λ(B). The set Λ(B) is a set of all integer linear combinations of the basis vectors, so Λ(B) = {xb1 + yb2 | x,y integers}.

In other words, we get a grid by calculating the integer combinations of the given basis vectors. For the standard basis E=[[1,0]^T, [0,1]^T] we can write the point [8,4]^T as 8 * [1,0]^T + 4 * [0,1]^T where both 8 and 4 are integers. The set of all integer combinations (hence the lattice Λ) then looks like this:

If we change the basis, this will result into a different lattice. Here is an example for b1=[2,3]^T, b2=[4,5]^T:

In order to produce these images I am using the following Python code:

import numpy as np
import matplotlib.pyplot as plt
from typing import List, Tuple


def plotLattice(ax: plt.Axes, basis_vectors: List[np.ndarray], 
                   ldown: np.ndarray, rup: np.ndarray, color: str, 
                   linewidth: float, alpha: float) -> List[np.ndarray]:
    """
    Draws a two-dimensional lattice.

    Args:
        ax: The Matplotlib Axes instance to plot on.
        basis_vectors: A list of two NumPy arrays representing the basis vectors of the lattice.
        ldown: A NumPy array representing the lower left corner of the rectangular area to draw the lattice in.
        rup: A NumPy array representing the upper right corner of the rectangular area to draw the lattice in.
        color: A string representing the color of the lattice points and basis vectors.
        linewidth: A float representing the linewidth of the lattice points.
        alpha: A float representing the alpha value of the lattice points.

    Returns:
        A list of NumPy arrays representing the lattice points.
    """
    
    # get the basis vectors
    b1, b2 = np.array(basis_vectors[0]), np.array(basis_vectors[1])
    # list to hold the lattice points
    points = []
    
    # upper bounds for the for loops
    xmax, ymax = 0, 0
    if b1[0] == 0:
        xmax = np.floor(rup[0] / abs(b2[0]))
    elif b2[0] == 0:
        xmax = np.floor(rup[0] / abs(b1[0]))
    else:
        xmax = np.floor(rup[0] / min(abs(b1[0]),abs(b2[0])))
    
    if b1[1] == 0:
        ymax = np.floor(rup[1] / abs(b2[1]))
    elif b2[1] == 0:
        ymax = np.floor(rup[1] / abs(b1[1]))
    else:
        ymax = np.floor(rup[1] / min(abs(b1[1]),abs(b2[1])))
    
    # increase the bounds by 1
    xmax = int(xmax) + 1
    ymax = int(ymax) + 1
    
    # get the lower bounds for the for loops
    xmin, ymin = -int(xmax), -int(ymax)
    
    for i in range(xmin, int(xmax)):
        for j in range(ymin, int(ymax)):
            # make the linear combination
            p = i * b1 + j * b2
            # if the point is within the plotting area, plot it and add the point to the list
            if ldown[0] <= p[0] <= rup[0] and ldown[1] <= p[1] <= rup[1]:
                ax.scatter(p[0], p[1], color=color, linewidths=linewidth, alpha=alpha)
                points.append(p)

    # plot basis vectors
    ax.quiver(0, 0, b1[0], b1[1], color=color, scale_units='xy', scale=1, alpha=1)
    ax.quiver(0, 0, b2[0], b2[1], color=color, scale_units='xy', scale=1, alpha=1)

    return points


if __name__ == '__main__':
    # pick a basis
    b1 = np.array([2, 3])
    b2 = np.array([-4, 5])
    basis_vectors = [b1, b2]
    
    # define the plotting area
    ldown = np.array([-15, -15])
    rup = np.array([15, 15])

    fig, ax = plt.subplots()
    points = plotLattice(ax, basis_vectors, ldown, rup, 'blue', 3, 0.25)

    # resize the plotting window
    mngr = plt.get_current_fig_manager()
    mngr.resize(960, 1080)
    
    # tune axis
    ax.set_aspect('equal')
    ax.grid(True, which='both')
    ax.set_xlim(ldown[0], rup[0])
    ax.set_ylim(ldown[1], rup[1])
    
    # show the plot
    plt.show()

And now we get to the problem. For the basis vectors b1=[1,1], b2=[1,2] the code does not cover the plotting area:

We can increase the problem by choosing some not nicely orthogonal vectors:

So, the problem arises every time when the vectors are getting closer to each other, hence when the dot product is big. Now consider the example I picked before:

My approach was to take the minimum values of the absolute coordinate values and to estimate how much points I can have along one axis. This is also the approach you can see in the code. Taking the minimum of the x coordinates of b1=[1,1] and b2=[1,2] we get 1. Taking the minimum of the y coordinates we get 1. My plotting area is defined by the square which is given by the points ldown=[-15,-15] and rup=[15,15]. Hence I know that I can have maximum floor(rup[0]/1) = 15 points along the x-axis, and maximum floor(rup[1]/1) = 15 along the y-axis. Including the zero point it results in 31 points for each axis, so that I expect to see 31*31 = 961 points on the plot. Hence, I think that I am done and set xmax=15, xmin=-15, ymax=15, ymin=-15.

But this gives me the result presented above. So, the calculation is wrong. Then I say, "Ok, I know that the point [15,-15] has to be in the plot". Hence I can solve the system Bx = [15,-15]^T. This results into the vector x=[45,-30]. Now I can set xmax=45, ymin=-30. Doing the same for the point [-15,15] gives me the vector x=[-45,30]. So I can set xmin=-45, ymin=-30. The resulting plot is:

This looks almost well, but you can notice that the points [15,-15] and [-15,15] are missing in the plot. Hence I have to enlarge the bounds by 1 by setting xmax=46, xmin=-46, ymax=31, ymin=-31. After that, the whole area is covered.

So, the drawback of this mechanism, is that I cheated a bit. Here, I just knew that the point [15,-15] must be on the plot. I could solve the equations system and determine the bounds for the for loop. Occasionally, this point was also the most distant point from the origin, so that I knew that covering it I should automatically cover the whole plotting plane. However, there are basis vectors for which I cannot determine such point, and we can just pick one of the previous plots:

Here, my approach would say that we can have min(2,4) = 2 points along the x-axis and min(3,5)=3 points along the y-axis. But I cannot simply say that the point [14,-9]=[7*2,-3*3] is on the plot (because it is not). Moreover, I cannot even say which of the points [12,-12], [12,-15], [14,-12],[14-15] are part of the plot, and which are not. Knowing the plot I see that [12,-15] and [14,-12] must be in the plot. Without knowing that I do not even know for which point I have to solve the Bx=b system.

Choosing different basis or a different (not origin-centered) plotting area makes the problem surprisingly complex for me, - even though we are acting in a 2D plane only.

So, now when the problem is described, I can formulate it as follows: Given the points rup and ldown of the plotting area, a basis b1, b2, define the bounds xmax, xmin, ymax, ymin for the for loops, so that the whole plotting area gets covered by the lattice points.

I am not even asking the code to be efficient at the moment, however a solution of the type xmax = sys.maxsize or xmax = 100 * rup[0] do not count.

What would be your approach?

Here is an outline of a solution.

Your problem is that the 2 vectors from which you build your grid are not necessarily the closest points in your grid. So we want to discover the closest points. After we know that, we can build the grid by taking each point on the grid, and adding its closest points, then theirs and so on until we've covered the whole plot area.

So how do we find the closest points in the grid?

What we're going to do is keep adding points to the queue, prioritizing those close to the origin. We keep going until every point is either in a list of minimal points, or is the sum or difference of two closer points.

To prioritize we can use a priority queue.

Now let's take your example (4, 1) and (4, 2) as the starting vectors and see how this works. We'll have a queue upcoming of points we'll look for. I'll just write it in sorted order to make it clear.

We'll have a list of examined points.

We'll have a set of not minimal points.

We start with something like:

upcoming = [Point(4, 1), Point(4, 2)]
examined = []
not_minimal = set([])

Now we take out the first point, add it to examined, add the sum or difference of it and all examined points either to upcoming or not_minimal depending on whether it is nearer or farther from the points we added.

And now our code is something like this (untested)

while 0 < len(upcoming):
    point = heapq.heappop(upcoming)
    if point not in not_minimal: # Might have found another way to get it.
        examined.append(point)
        for point2 in examined:
            for point3 in [
                point + point2,
                point - point2,
                -point + point2,
                -point - point2
            ]):
                if max(distance(point), distance(point2)) < distance(point3):
                    not_minimal.add(point3)
                else:
                    heapq.heappush(upcoming, point3)
minimal_points = [point for point in examined if point not in not_minimal]

And now you will wind up with the nearest points to the origin. Now start with a point you want in your final answer, and build up the grid around it. You may, depending on the grid, need to go outside your plotting area by max((distance(point) for point in minimal_points)) to discover some of the corner points in your grid. But you should miss none.

The initial discovery of finite points takes a finite time bounded above by time O(m log(m)) where m is how many vectors are in a fixed circle/sphere/whatever around the origin of radius double your larger vector.

The discovery of the full grid takes time O(n * k) where k is the number of minimal points and n is the number of points in your grid that you need to find.

This algorithm should work in any number of dimensions.

---

And for fun I coded it up. This just returns the points that you want, instead of drawing them. Fixing that is easy. It also will handle more than 2 dimensions. Drawing that will take more work.

import numpy as np
import heapq

def distance_to_box (ldown, lup, point):
    deltas = []
    for i in range(len(ldown)):
        (x_min, x_max) = sorted([ldown[i], lup[i]])
        if point[i] < x_min:
            deltas.append(x_min - point[i])
        elif x_max < point[i]:
            deltas.append(point[i] - x_max)
        else:
            deltas.append(0)
    da = np.array(deltas)
    return np.linalg.norm(da)

def grid (basis_vectors, ldown, lup, epsilon = 0.1**8):
    upcoming = []
    i = 0
    for v in basis_vectors:
        point = np.array(v)
        i += 1
        heapq.heappush(upcoming, (np.dot(point, point), i, point))

    minimal = {}
    not_minimal = {}
    while len(upcoming):
        (d, _, p) = heapq.heappop(upcoming)
        if str(p) not in not_minimal:
            minimal[str(p)] = p
            to_delete = []
            for s2, p2 in minimal.items():
                if s2 in not_minimal:
                    to_delete.append(s2)
                else:
                    d2 = np.linalg.norm(p2)
                    for p3 in [p + p2, p - p2, -p + p2, -p - p2]:
                        d3 = np.linalg.norm(p3)
                        if max(d, d2) + epsilon < d3:
                            not_minimal[str(p3)] = p3
                        elif d3 < epsilon:
                            pass # ignore 0 vector
                        else:
                            i += 1
                            heapq.heappush(upcoming, (d3, i, p3))
            for s2 in to_delete:
                minimal.pop(s2)
    directions = minimal.values()

    start = np.array([0 for _ in ldown])
    improved = True
    while improved:
        improved = False
        for direction in directions:
            step = direction
            # Does walking in this direction help?
            while np.linalg.norm(start + step - ldown) < np.linalg.norm(start - ldown):
                improved = True
                start = start + step
                step = step + step

    todo = [start]
    seen = {}
    max_d = max([np.linalg.norm(d) for d in directions]) + epsilon
    answer = []
    while len(todo):
        p = todo.pop()
        d = distance_to_box(ldown, lup, p)
        if d < max_d:
            if str(p) not in seen:
                for direction in directions:
                    todo.append(p + direction)
                if d < epsilon:
                    answer.append(p)
                seen[str(p)] = p
    return answer



for p in grid([(4, 1), (4, 2)], (-4, -4), (4, 4)):
    print(p)