title: Efficiently generate a lattice of points in python
body: <p><strong>Help make my code faster</strong>: My python code needs to generate a 2D lattice of points that fall inside a bounding rectangle. I kludged together some code (shown below) that generates this lattice. However, this function is called many many times and has become a serious bottleneck in my application.</p>
  <p>I'm sure there's a faster way to do this, probably involving numpy arrays instead of lists. Any suggestions for a faster, more elegant way to do this?</p>
  <p><strong>Description of the function</strong>:
  I have two 2D vectors, v1 and v2. These vectors <a href="http://en.wikipedia.org/wiki/Bravais_lattice#Bravais_lattices_in_at_most_2_dimensions" rel="nofollow noreferrer">define a lattice</a>. In my case, my vectors define a lattice that is almost, but not quite, hexagonal. I want to generate the set of all 2D points on this lattice that are in some bounding rectangle. In my case, one of the corners of the rectangle is at (0, 0), and the other corners are at positive coordinates.</p>
  <p><strong>Example</strong>:
  If the far corner of my bounding rectangle was at (3, 3), and my lattice vectors were:</p>
  <pre><code>v1 = (1.2, 0.1)
  v2 = (0.2, 1.1)
  </code></pre>
  <p>I'd want my function to return the points:</p>
  <pre><code>(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
  </code></pre>
  <p>I'm not concerned with edge cases; it's doesn't matter if the function returns (0, 0), for example.</p>
  <p><strong>The slow way I'm currently doing it</strong>:</p>
  <pre><code>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 &lt; lp) and all(lp &lt; 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()
  </code></pre>