# Finding the best matching block/patch in Python

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

• What constitutes the 'closet matching block'? – BeetDemGuise Jun 19 '14 at 11:49
• Here I am using the sum of absolute difference. Pixel by pixel comparison in the NxN blocks. – user3752615 Jun 19 '14 at 12:30
• Could you give an approximation of your current running time and give us a target to achieve ? It would help to better compare future reviews. – Marc-Andre Jun 19 '14 at 12:32
• Do you consider using numpy an option? – Jaime Jun 19 '14 at 13:57
• @Marc-Andre This function within my program takes around 8 minutes to run on my machine. Note that this function is called 30,744 times. I've used a C++ implementation that is able to achieve this in mere seconds. – user3752615 Jun 19 '14 at 14:10

## 2 Answers

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)


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.