Эффективно сгенерируйте решетку точек в python

Помогите ускорить мой код : Мой код на Python должен генерировать двумерную решетку точек, которые попадают в ограничивающий прямоугольник. Я собрал код (показанный ниже), который генерирует эту решетку. Однако эта функция вызывается много раз, и она стала серьезным узким местом в моем приложении.

Я уверен, что есть более быстрый способ сделать это, вероятно, с использованием массивов numpy вместо списков. Есть ли предложения по более быстрому и элегантному способу сделать это?

Описание функции : У меня есть два 2D-вектора, v1 и v2. Эти векторы определяют решетку . В моем случае мои векторы определяют решетку, которая почти, но не совсем гексагональна. Я хочу создать набор всех 2D-точек на этой решетке, которые находятся в некотором ограничивающем прямоугольнике. В моем случае один из углов прямоугольника находится в точке (0, 0), а другие углы имеют положительные координаты.

Пример : Если бы дальний угол моего ограничивающего прямоугольника находился в точке (3, 3), а мои векторы решетки были:

v1 = (1.2, 0.1)
v2 = (0.2, 1.1)

Я бы хотел, чтобы моя функция возвращала точки:

(1.2, 0.1) #v1
(2.4, 0.2) #2*v1
(0.2, 1.1) #v2
(0.4, 2.2) #2*v2
(1.4, 1.2) #v1 + v2
(2.6, 1.3) #2*v1 + v2
(1.6, 2.3) #v1 + 2*v2
(2.8, 2.4) #2*v1 + 2*v2

Меня не интересуют крайние случаи; например, не имеет значения, возвращает ли функция (0, 0).

Медленный способ, которым я сейчас это делаю :

import numpy, pylab

def generate_lattice( #Help me speed up this function, please!
    image_shape, lattice_vectors, center_pix='image', edge_buffer=2):

    ##Preprocessing. Not much of a bottleneck:
    if center_pix == 'image':
        center_pix = numpy.array(image_shape) // 2
    else: ##Express the center pixel in terms of the lattice vectors
        center_pix = numpy.array(center_pix) - (numpy.array(image_shape) // 2)
        lattice_components = numpy.linalg.solve(
            numpy.vstack(lattice_vectors[:2]).T,
            center_pix)
        lattice_components -= lattice_components // 1
        center_pix = (lattice_vectors[0] * lattice_components[0] +
                      lattice_vectors[1] * lattice_components[1] +
                      numpy.array(image_shape)//2)
    num_vectors = int( ##Estimate how many lattice points we need
        max(image_shape) / numpy.sqrt(lattice_vectors[0]**2).sum())
    lattice_points = []
    lower_bounds = numpy.array((edge_buffer, edge_buffer))
    upper_bounds = numpy.array(image_shape) - edge_buffer

    ##SLOW LOOP HERE. 'num_vectors' is often quite large.
    for i in range(-num_vectors, num_vectors):
        for j in range(-num_vectors, num_vectors):
            lp = i * lattice_vectors[0] + j * lattice_vectors[1] + center_pix
            if all(lower_bounds < lp) and all(lp < upper_bounds):
                lattice_points.append(lp)
    return lattice_points


##Test the function and display the output.
##No optimization needed past this point.
lattice_vectors = [
    numpy.array([-40., -1.]),
    numpy.array([ 18., -37.])]
image_shape = (1000, 1000)
spots = generate_lattice(image_shape, lattice_vectors)

fig=pylab.figure()
pylab.plot([p[1] for p in spots], [p[0] for p in spots], '.')
pylab.axis('equal')
fig.show()