python - Fast Way to slice image into overlapping patches and merge patches to image

Question

Trying to slice a grayscale image of size 100x100 into patches of size 39x39 which are overlapping, with a stride-size of 1. That means that the next patch which starts one pixel to the right/or below is only different to the previous patch in one additional column/or row.

Rough Outline of the code: First compute the indices for each patch, so to be able to construct the 2D array of patches from an image and to be able to construct an image from patches:

patches = imgFlat[ind]

'patches' is a 2D array with each column containing a patch in vector form.

Those patches are processed, each patch individually and afterwards are merged together to an image again, with the precomputed indices.

img = np.sum(patchesWithColFlat[ind],axis=2)

As patches overlap it is necessary at the end to multiply the img with the precomputed weights:

imgOut = weights*imgOut

My code is really slow and speed is a critical issue, as this should be done on ca. 10^8 patches.

The functions get_indices_for_un_patchify and weights_unpatchify can be precomputed once, so speed is only an issue for patchify and unpatchify.

Thanks for any tipps.

Carlos

import numpy as np
import scipy
import collections
import random as rand


def get_indices_for_un_patchify(sImg,sP,step):
    ''' creates indices for fast patchifying and unpatchifying

    INPUTS:
      sx    image size
      sp    patch size
      step  offset between two patches (default == [1,1])

      OUTPUTS:
       patchInd             collection with indices
       patchInd.img2patch   patchifying indices
                            patch = img(patchInd.img2patch);
       patchInd.patch2img   unpatchifying indices

    NOTE: * for unpatchifying necessary to add a 0 column to the patch matrix
          * matrices are constructed row by row, as normally there are less rows than columns in the 
            patchMtx
     '''
    lImg = np.prod(sImg)
    indImg = np.reshape(range(lImg), sImg)

    # no. of patches which fit into the image
    sB = (sImg - sP + step) / step

    lb              = np.prod(sB)
    lp              = np.prod(sP)
    indImg2Patch    = np.zeros([lp, lb])
    indPatch        = np.reshape(range(lp*lb), [lp, lb])

    indPatch2Img = np.ones([sImg[0],sImg[1],lp])*(lp*lb+1)

    # default value should be last column
    iRow   = 0;
    for jCol in range(sP[1]):
        for jRow in range(sP[0]):
            tmp1 = np.array(range(0, sImg[0]-sP[0]+1, step[0]))
            tmp2 = np.array(range(0, sImg[1]-sP[1]+1, step[1]))
            sel1                    = jRow  + tmp1
            sel2                    = jCol  + tmp2
            tmpIndImg2Patch = indImg[sel1,:]          
            # do not know how to combine following 2 lines in python
            tmpIndImg2Patch = tmpIndImg2Patch[:,sel2]
            indImg2Patch[iRow, :]   = tmpIndImg2Patch.flatten()

            # next line not nice, but do not know how to implement it better
            indPatch2Img[min(sel1):max(sel1)+1, min(sel2):max(sel2)+1, iRow] = np.reshape(indPatch[iRow, :, np.newaxis], sB)
            iRow                    += 1

    pInd = collections.namedtuple
    pInd.patch2img = indPatch2Img
    pInd.img2patch = indImg2Patch

    return pInd

def weights_unpatchify(sImg,pInd):
    weights = 1./unpatchify(patchify(np.ones(sImg), pInd), pInd)
    return weights

# @profile
def patchify(img,pInd):
    imgFlat = img.flat
   # imgFlat = img.flatten()
    ind = pInd.img2patch.tolist()
    patches = imgFlat[ind]

    return patches

# @profile
def unpatchify(patches,pInd):
    # add a row of zeros to the patches matrix    
    h,w = patches.shape
    patchesWithCol = np.zeros([h+1,w])
    patchesWithCol[:-1,:] = patches
    patchesWithColFlat = patchesWithCol.flat
   #  patchesWithColFlat = patchesWithCol.flatten()
    ind = pInd.patch2img.tolist()
    img = np.sum(patchesWithColFlat[ind],axis=2)
    return img

I call those functions, here e.g. with a random image

if __name__ =='__main__':
    img = np.random.randint(255,size=[100,100])
    sImg = img.shape
    sP = np.array([39,39])  # size of patch
    step = np.array([1,1])  # sliding window step size
    pInd = get_indices_for_un_patchify(sImg,sP,step)
    patches = patchify(img,pInd)
    imgOut = unpatchify(patches,pInd)
    weights = weights_unpatchify(sImg,pInd)
    imgOut = weights*imgOut

    print 'Difference of img and imgOut = %.7f' %sum(img.flatten() - imgOut.flatten())

Answer 1

An efficient way to "patchify" an array, that is, to get an array of windows to the original array is to create a view with custom strides, the number of bytes to jump to the following element. It can be helpful to think of a numpy array as a (glorified) chunk of memory, and then strides are a way to map indices to memory address.

For example, in

a = np.arange(10).reshape(2, 5)

a.itemsize equals 4 (ie, 4 bytes or 32 bits for each element) and a.strides is (20, 4) (5 elements, 1 element) so that a[1,2] refers to the element which is 1*20 + 2*4 bytes (or 1*5 + 2 elements) after the first one:

0 1 2 3 4
5 6 7 x x

In fact, the elements are placed in memory one after another, 0 1 2 3 4 5 6 7 x x but the strides let us index it as a 2D array.

Building on that concept, we can rewrite patchify as follows

def patchify(img, patch_shape):
    img = np.ascontiguousarray(img)  # won't make a copy if not needed
    X, Y = img.shape
    x, y = patch_shape
    shape = ((X-x+1), (Y-y+1), x, y) # number of patches, patch_shape
    # The right strides can be thought by:
    # 1) Thinking of `img` as a chunk of memory in C order
    # 2) Asking how many items through that chunk of memory are needed when indices
    #    i,j,k,l are incremented by one
    strides = img.itemsize*np.array([Y, 1, Y, 1])
    return np.lib.stride_tricks.as_strided(img, shape=shape, strides=strides)

This function returns a view of img, so no memory is allocated and it runs in just a few tens of microseconds. The output shape is not exactly what you want, and in fact it must be copied to be able to get that shape.

One has to be careful when dealing with views of an array that are much larger than the base array, because operations can trigger a copy which would need to allocate a lot of memory. In your case, since the arrays isn't too big and there aren't that many patches, it should be ok.

Finally, we can ravel the array of patches a bit:

patches = patchify(img, (39,39))
contiguous_patches = np.ascontiguousarray(patches)
contiguous_patches.shape = (-1, 39**2)

This doesn't reproduce the output of your patchify function because you kind of develop the patches in Fortran order. I recommend you to use this instead because

1. It leads to more natural indexing later on (ie, the first patch is patches[0] instead of patches[:, 0] for your solution).

2. It's also easier in numpy to use C ordering everywhere because you need less typing (you avoid stuff like order='F', arrays are created in C order by default...).

"Hints" in case you insist: strides = img.itemsize * np.array([1, Y, Y, 1]), use .reshape(..., order='F') on contiguous_patches and finally transpose it .T