高效地生成python中点的点阵

Question

帮助我更快地编写代码：我的python代码需要生成位于边界矩形内的点的二维点阵。我将一些代码（如下所示）汇集在一起​​生成此格。然而，这个函数被多次调用，并且已经成为我应用程序中的一个严重瓶颈。

我相信有一个更快的方法来做到这一点，可能涉及numpy数组而不是列表。任何建议更快，更优雅的方式来做到这一点？

功能描述： 我有两个2D矢量，v1和v2。这些载体define a lattice。就我而言，我的矢量定义了一个几乎但不完全是六角形的格子。我想生成一些边界矩形中的这个点阵上的所有二维点的集合。在我的情况下，矩形的其中一个角位于（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() 

Answer 1

如果你想要矢量化整个事物，生成一个正方形格子然后剪切它。然后切掉落在你盒子外面的边缘。

这是我想出来的。还有很多可以改进的地方，但这是基本的想法。

def generate_lattice(image_shape, lattice_vectors) : 
    center_pix = numpy.array(image_shape) // 2 
    # Get the lower limit on the cell size. 
    dx_cell = max(abs(lattice_vectors[0][0]), abs(lattice_vectors[1][0])) 
    dy_cell = max(abs(lattice_vectors[0][1]), abs(lattice_vectors[1][1])) 
    # Get an over estimate of how many cells across and up. 
    nx = image_shape[0]//dx_cell 
    ny = image_shape[1]//dy_cell 
    # Generate a square lattice, with too many points. 
    # Here I generate a factor of 4 more points than I need, which ensures 
    # coverage for highly sheared lattices. If your lattice is not highly 
    # sheared, than you can generate fewer points. 
    x_sq = np.arange(-nx, nx, dtype=float) 
    y_sq = np.arange(-ny, nx, dtype=float) 
    x_sq.shape = x_sq.shape + (1,) 
    y_sq.shape = (1,) + y_sq.shape 
    # Now shear the whole thing using the lattice vectors 
    x_lattice = lattice_vectors[0][0]*x_sq + lattice_vectors[1][0]*y_sq 
    y_lattice = lattice_vectors[0][1]*x_sq + lattice_vectors[1][1]*y_sq 
    # Trim to fit in box. 
    mask = ((x_lattice < image_shape[0]/2.0) 
      & (x_lattice > -image_shape[0]/2.0)) 
    mask = mask & ((y_lattice < image_shape[1]/2.0) 
        & (y_lattice > -image_shape[1]/2.0)) 
    x_lattice = x_lattice[mask] 
    y_lattice = y_lattice[mask] 
    # Translate to the centre pix. 
    x_lattice += center_pix[0] 
    y_lattice += center_pix[1] 
    # Make output compatible with original version. 
    out = np.empty((len(x_lattice), 2), dtype=float) 
    out[:, 0] = y_lattice 
    out[:, 1] = x_lattice 
    return out 

Answer 2

由于lower_bounds和upper_bounds只有2个元素组成的数组，numpy的可能不是这里的正确选择。尝试使用基本的Python的东西来代替

if all(lower_bounds < lp) and all(lp < upper_bounds): 

：

if lower1 < lp and lower2 < lp and lp < upper1 and lp < upper2: 

根据timeit，第二种方法要快得多：

>>> timeit.timeit('all(lower < lp)', 'import numpy\nlp=4\nlower = numpy.array((1,5))') 
3.7948939800262451 
>>> timeit.timeit('lower1 < 4 and lower2 < 4', 'lp = 4\nlower1, lower2 = 1,5') 
0.074192047119140625 

从我的经验，只要您不需要处理n-diminsional数据，如果你不需要双精度浮点数，使用基本的Python数据类型和构造inst通常会更快ead of numpy，在这种情况下有点过载 - 请看this other question。

---

另一个小的改进可能是只计算一次range(-num_vectors, num_vectors)然后重新使用它。此外，您可能希望使用product iterator而不是嵌套循环 - 尽管我不认为这些更改会对性能产生重大影响。

Answer 3

也许你可以用这个替换两个for循环。

i,j = numpy.mgrid[-num_vectors:num_vectors, -num_vectors:num_vectors] 
numel = num_vectors ** 2; 
i = i.reshape(numel, 1) 
j = j.reshape(numel, 1) 
lp = i * lattice_vectors[0] + j * lattice_vectors[1] + center_pix 
valid = numpy.all(lower_bounds < lp, 1) and numpy.all(lp < upper_bounds, 1) 
lattice_points = lp[valid] 

可能会有一些小错误，但您明白了。

编辑

我做编辑的 “numpy.all（lower_bounds ..）” 考虑到正确的尺寸。

Answer 4

我得到了比2X加速更通过重复添加而不是乘法免去您的lp计算。 xrange优化似乎是无足轻重的（尽管它可能不会受伤）;重复添加似乎比乘法更有效。将这与上面提到的其他优化结合起来应该会给你更多的加速。但是，当然，最好的你可以得到一个恒定的因素加速，因为你的输出的大小是二次的，就像你的原始代码一样。

lv0, lv1 = lattice_vectors[0], lattice_vectors[1] 
lv0_i = -num_vectors * lv0 + center_pix 
lv1_orig = -num_vectors * lv1 
lv1_j = lv1_orig 
for i in xrange(-num_vectors, num_vectors): 
    for j in xrange(-num_vectors, num_vectors): 
     lp = lv0_i + lv1_j 
     if all(lower_bounds < lp) and all(lp < upper_bounds): 
      lattice_points.append(lp) 
     lv1_j += lv1 
    lv0_i += lv0 
    lv1_j = lv1_orig 

定时器结果：

>>> t = Timer("generate_lattice(image_shape, lattice_vectors, orig=True)", "from __main__ import generate_lattice, lattice_vectors, image_shape") 
>>> print t.timeit(number=50) 
5.20121788979 
>>> t = Timer("generate_lattice(image_shape, lattice_vectors, orig=False)", "from __main__ import generate_lattice, lattice_vectors, image_shape") 
>>> print t.timeit(number=50) 
2.17463898659