Finding the best matching block/patch in Python

Question

I wish to locate the closest matching NxN block within a WxW window centred at location (x,y) of a larger 2D array. The code below works fine but is very slow for my needs as I need to run this operation many times.

Is there a better way to do this?

Here N = 3, W = 15, x =15, y = 15 and (bestx, besty) is the centre of the best matching block

import numpy as np

## Generate some test data
CurPatch = np.random.randint(20, size=(3, 3))
Data = np.random.randint(20,size=(30,30))

# Current Location 
x,y = 15,15
# Initialise Best Match
bestcost = 999.0
bestx = 0;besty=0

for Wy in xrange(-7,8):
    for Wx in xrange(-7,8):
            Ywj,Ywi = y+Wy,x+Wx 

            cost = 0.0
            for py in xrange(3):
                for px in xrange(3):
                    cost += abs(Data[Ywj+py-1,Ywi+px-1] - CurPatch[py,px]) 
                    if cost >= bestcost:
                       break

            if cost < bestcost:
                bestcost = cost
                besty,bestx = Ywj,Ywi

print besty,bestx

Answer 1

The following outlines a couple variations of a numpy solution. For the problem sizes in your post, they all clock faster than 100us, so you should get your 30,744 calls in about 3 seconds.

import numpy as np

N = 3

W = 15
x, y = 15, 15

data = np.random.randint(20, size=(30, 30))
current_patch = np.random.randint(20, size=(N, N))

# isolate the patch
x_start = x - W//2
y_start = y - W//2

patch = data[x_start:x_start+W, y_start:y_start+W]

# take a windowed view of the array
from numpy.lib.stride_tricks import as_strided
shape = tuple(np.subtract(patch.shape, N-1)) + (N, N)
windowed_patch = as_strided(patch, shape=shape, strides=patch.strides*2)

# this works, but creates a large intermediate array
cost = np.abs(windowed_patch - current_patch).sum(axis=(-1, -2))

# this is much more memory efficient, but uses squared differences,
# and the fact that (a - b)^2 = a^2 + b^2 - 2*a*b
patch2 = patch*patch
ssqd =  as_strided(patch2, shape=shape,
                   strides=patch2.strides*2).sum(axis=(-1, -2),
                                                 dtype=np.double)
ssqd += np.sum(current_patch * current_patch, dtype=np.double)
ssqd -= 2 * np.einsum('ijkl,kl->ij', windowed_patch, current_patch,
                      dtype=np.double)

# for comparison with the other method
cost2 = windowed_patch - current_patch
cost2 = np.sum(cost2*cost2, axis=(-1, -2))

# with any method, to find the location of the max in cost:
best_X, best_y = np.unravel_index(np.argmax(cost), cost.shape)

Answer 2

Using sum, you can make your code more concise and slightly more efficient :

cost = sum(abs(Data[Wy+py-1,Wx+px-1] - CurPatch[py,px]) for py in xrange(3) for px in xrange(3))

---

You can avoid manipulation of Ywj,Ywi by sliding the ranges you are using :

for Wy in xrange(y-7,y+8):
    for Wx in xrange(x-7,x+8):

---

By doing this, I got tiny tiny performance improvements.

Then, there are a few details of style (mostly about whitespaces) you could change to make your code more pythonic.