Matching numpy patterns using matrix multiplications

Question

I am trying to match a pattern with a binary image (black and white only) by moving the pattern along the image. And return the minimum distance between the template and the image with the corresponding position at which this minimum distance actually occurred. For instance:

IMG:

0 1 0 0 0 1 0 1 1 

template:

 0 1 1 1 

This pattern matches the best image in position (1,1), and then it will be 0. Now everything is not so complicated, and I already have code that does the trick.

 def match_template(img, template): mindist = float('inf') idx = (-1,-1) for y in xrange(img.shape[1]-template.shape[1]+1): for x in xrange(img.shape[0]-template.shape[0]+1): #calculate Euclidean distance dist = np.sqrt(np.sum(np.square(template - img[x:x+template.shape[0],y:y+template.shape[1]]))) if dist < mindist: mindist = dist idx = (x,y) return [mindist, idx] 

But for images of the required size (an image between 500 x 200 pixels and a template of 250 x 100) it takes about 4.5 seconds, which is too slow. And I know that the same thing can be done much faster using matrix multiplications (in Matlab, I believe that this can be done with im2col and repmat). Can someone explain me how to do this in python / numpy?

by the way. I know that there is an opencv matchTemplate function that does exactly what I need, but since I may need to change the code a bit later, I would prefer a solution that I fully understand and can change.

Thanks!

edit: If anyone can explain to me how opencv does this in less than 0.2 seconds, that would also be great. I looked briefly at the source code, but these things somehow always look quite complicated for me.

edit2: Cython code

 import numpy as np cimport numpy as np DTYPE = np.int ctypedef np.int_t DTYPE_t def match_template(np.ndarray img, np.ndarray template): cdef float mindist = float('inf') cdef int x_coord = -1 cdef int y_coord = -1 cdef float dist cdef unsigned int x, y cdef int img_width = img.shape[0] cdef int img_height = img.shape[1] cdef int template_width = template.shape[0] cdef int template_height = template.shape[1] cdef int range_x = img_width-template_width+1 cdef int range_y = img_height-template_height+1 for y from 0 <= y < range_y: for x from 0 <= x < range_x: dist = np.sqrt(np.sum(np.square(template - img[ x:<unsigned int>(x+template_width), y:<unsigned int>(y+template_height) ]))) #calculate euclidean distance if dist < mindist: mindist = dist x_coord = x y_coord = y return [mindist, (x_coord,y_coord)] img = np.asarray(img, dtype=DTYPE) template = np.asarray(template, dtype=DTYPE) match_template(img, template) 

Answer 1

One possible way to do what you want is convolution (which can be brute force or FFT). AFAIK matrix multiplications will not work. You need to drill your data using a template. And find the maximum (you will also need to scale a little for it to work correctly).

 xs=np.array([[0,1,0],[0,0,1],[0,1,1]])*1. ys=np.array([[0,1],[1,1]])*1. print scipy.ndimage.convolve(xs,ys,mode='constant',cval=np.inf) >>> array([[ 1., 1., inf], [ 0., 2., inf], [ inf, inf, inf]]) print scipy.signal.fftconvolve(xs,ys,mode='valid') >>> array([[ 1., 1.], [ 0., 2.]]) 

Answer 2

A bizarre way to do this may occur, using pure numpy / scipy magic. But it might be easier (and more understandable when you look at the code in the future) to just get into Cython to do this. There's a good tutorial on integrating Cython with numpy at http://docs.cython.org/src/tutorial/numpy.html .

EDIT: I did a quick test with your Cython code, and it ran ~ 15 seconds for a 500x400 img with a 100x200 pattern. After some tweaks (excluding calls to the numpy method and checking for number limits), I got it in less than 3 seconds. This may not be enough for you, but it shows the possibility.

 import numpy as np cimport numpy as np cimport cython from libc.math cimport sqrt DTYPE = np.int ctypedef np.int_t DTYPE_t @cython.boundscheck(False) def match_template(np.ndarray[DTYPE_t, ndim=2] img, np.ndarray[DTYPE_t, ndim=2] template): cdef float mindist = float('inf') cdef int x_coord = -1 cdef int y_coord = -1 cdef float dist cdef unsigned int x, y cdef int img_width = img.shape[0] cdef int img_height = img.shape[1] cdef int template_width = template.shape[0] cdef int template_height = template.shape[1] cdef int range_x = img_width-template_width+1 cdef int range_y = img_height-template_height+1 cdef DTYPE_t total cdef int delta cdef unsigned int j, k, j_plus, k_plus for y from 0 <= y < range_y: for x from 0 <= x < range_x: #dist = np.sqrt(np.sum(np.square(template - img[ x:<unsigned int>(x+template_width), y:<unsigned int>(y+template_height) ]))) #calculate euclidean distance # Do the same operations, but in plain C total = 0 for j from 0 <= j < template_width: j_plus = <unsigned int>x + j for k from 0 <= k < template_height: k_plus = <unsigned int>y + k delta = template[j, k] - img[j_plus, k_plus] total += delta*delta dist = sqrt(total) if dist < mindist: mindist = dist x_coord = x y_coord = y return [mindist, (x_coord,y_coord)]