使用矩陣乘法的numpy模板匹配

Question

我想通過沿圖像移動模板來匹配模板和二進制圖像（只有黑色和白色）。然後將模板和圖像之間的最小距離返回到確實發生該最小距離的相應位置。例如：

IMG：

0 1 0 
0 0 1 
0 1 1 

模板：

0 1 
1 1 

這個模板圖像最擅長的位置（1,1）相匹配，然後在距離爲0。到目前爲止，事情不是很難，我已經有一些代碼可以做到這一點。

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] 

但對於這個已花費大約4.5秒，這是太慢我所需要的尺寸（250×100之間的圖像間500×200像素和模板）的圖像。而且我知道使用矩陣乘法可以更快地完成同樣的事情（在matlab中，我相信這可以使用im2col和repmat來完成）。任何人都可以解釋我如何在Python/numpy中做到這一點？

btw。我知道有一個opencv matchTemplate函數完全符合我的需求，但由於我稍後可能需要稍後修改代碼，所以我更願意使用我完全理解並可以改變的解決方案。

謝謝！

編輯：如果有人能解釋我如何在不到0.2秒的時間內做到這一點，這也將是偉大的。我對源代碼進行了簡短的介紹，但這些事情對我來說總是很複雜。

EDIT2：用Cython代碼

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

有可能開始使用純numpy的/ SciPy的魔法此做一個奇特的方式。但是，在將來可能會更容易（並且在將來查看代碼時更容易理解）放入Cython中以完成此操作。有一個很好的教程，用於將Cython與numpy集成在http://docs.cython.org/src/tutorial/numpy.html。

編輯： 我用你的Cython代碼做了一個快速測試，它用一個100x200模板運行了一個500x400 img〜15秒。經過一些調整（消除了numpy方法調用和numpy邊界檢查）後，我在3秒鐘內完成了調整。這對你來說可能還不夠，但它表明了可能性。

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)] 

Answer 2

通過卷積（可以是強力或FFT）做一個你想做的事情的一種可能方式。矩陣乘法AFAIK不起作用。您需要將數據與模板進行卷積。並找到最大值（您還需要進行一些縮放以使其正常工作）。

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.]])