code: def block(shape, block_shape):
      
      shape = np.array(shape)
      block_shape = np.array(block_shape)
      i,j = np.mgrid[0:shape[0],0:shape[1]]
      ijmax = (shape.astype(float)/block_shape.astype(float)).astype(int)
      ijmax = np.maximum(ijmax, 1)
      multiplier = ijmax.astype(float) / shape.astype(float)
      i = (i * multiplier[0]).astype(int)
      j = (j * multiplier[1]).astype(int)
      labels = i * ijmax[1] + j
      indexes = np.array(list(range(np.product(ijmax))))
      return labels, indexes
docs: Create a labels image that divides the image into blocks
      
      shape - the shape of the image to be blocked
      block_shape - the shape of one block
      
      returns a labels matrix and the indexes of all labels generated
      
      The idea here is to block-process an image by using SciPy label
      routines. This routine divides the image into blocks of a configurable
      dimension. The caller then calls scipy.ndimage functions to process
      each block as a labeled image. The block values can then be applied
      to the image via indexing. For instance:
      
      labels, indexes = block(image.shape, (60,60))
      minima = scind.minimum(image, labels, indexes)
      img2 = image - minima[labels]