6
\$\begingroup\$

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
\$\endgroup\$
6
  • \$\begingroup\$ What constitutes the 'closet matching block'? \$\endgroup\$ – BeetDemGuise Jun 19 '14 at 11:49
  • \$\begingroup\$ Here I am using the sum of absolute difference. Pixel by pixel comparison in the NxN blocks. \$\endgroup\$ – user3752615 Jun 19 '14 at 12:30
  • \$\begingroup\$ 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. \$\endgroup\$ – Marc-Andre Jun 19 '14 at 12:32
  • \$\begingroup\$ Do you consider using numpy an option? \$\endgroup\$ – Jaime Jun 19 '14 at 13:57
  • \$\begingroup\$ @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. \$\endgroup\$ – user3752615 Jun 19 '14 at 14:10
3
\$\begingroup\$

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)
\$\endgroup\$
0
\$\begingroup\$

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.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.