导入模块
首先,让我们导入此帖子中列出的各种方法中将使用的所有相关模块/函数 -
from skimage.util import view_as_windows
from skimage.feature import match_template
import cv2
from cv2 import matchTemplate as cv2m
from scipy.ndimage.filters import uniform_filter as unif2d
from scipy.signal import convolve2d as conv2
基于视图的方法来计算MAD和MSD
使用Skimage进行计算平均绝对偏差的方法:
使用scikit-image获取滑动4D视图数组
,然后使用np.mean
来进行平均值计算 -
def skimage_views_MAD_v1(img, tmpl):
return np.abs(view_as_windows(img, tmpl.shape) - tmpl).mean(axis=(2,3))
使用scikit-image获取滑动4D数组的视图,然后使用np.einsum
进行平方平均计算 -
def skimage_views_MAD_v2(img, tmpl):
subs = np.abs(view_as_windows(img, tmpl.shape) - tmpl)
return np.einsum('ijkl->ij',subs)/float(tmpl.size)
使用Skimage计算均方差的方法:
使用相似技术,我们有两种计算均方差的方法 -
def skimage_views_MSD_v1(img, tmpl):
return ((view_as_windows(img, tmpl.shape) - tmpl)**2).mean(axis=(2,3))
def skimage_views_MSD_v2(img, tmpl):
subs = view_as_windows(img, tmpl.shape) - tmpl
return np.einsum('ijkl,ijkl->ij',subs, subs)/float(tmpl.size)
B. 基于卷积的MSD方法
用于计算均方偏差的基于卷积的方法:
卷积
可以被用来在略微调整的方式下计算均方偏差。对于平方偏差和,我们在每个窗口内执行元素级别的减法,然后对每个减法的结果进行平方并将它们全部加起来。
让我们考虑一个一维示例来更仔细地看一下 -
a : [a1, a2, a3, a4, a5, a6, a7, a8] # Image array
b : [b1, b2, b3] # Template array
对于第一个窗口操作,我们会有:
(a1-b1)**2 + (a2-b2)**2 + (a3-b3)**2
让我们使用(a-b)**2
公式:
(a - b)**2 = a**2 - 2*a*b +b**2
因此,我们对于第一个窗口将会有:
(a1**2 - 2*a1*b1 +b1**2) + (a2**2 - 2*a2*b2 +b2**2) + (a3**2 - 2*a3*b3 +b3**2)
同样地,对于第二个窗口:
(a2**2 - 2*a2*b1 +b1**2) + (a3**2 - 2*a3*b2 +b2**2) + (a4**2 - 2*a4*b3 +b3**2)
等等。
因此,我们需要对这些计算进行三个部分的划分 -
a 的平方并在滑动窗口中求和。
b 的平方并求和。这在所有窗口中保持不变。
拼图的最后一部分是:(2*a1*b1, 2*a2*b2, 2*a3*b3)
, (2*a2*b1, 2*a3*b2, 2*a4*b3)
以此类推,对于第一个、第二个等窗口。这可以通过使用 a
和 b
的翻转版本来进行 2D
卷积计算,根据 卷积
的定义运行内核从另一个方向滑动,并在每个窗口内计算元素乘法并求和元素,因此需要进行翻转。
将这些思想扩展到 2D
情况,并使用 Scipy的convolve2d
和 uniform_filter
,我们将有另外两种方法来计算 均方差
,如下所示 -
def convolution_MSD(img, tmpl):
n = tmpl.shape[0]
sums = conv2(img**2,np.ones((n,n)),'valid')
out = sums + (tmpl**2).sum() -2*conv2(img,tmpl[::-1,::-1],'valid')
return out/(n*n)
def uniform_filter_MSD(img, tmpl):
n = tmpl.shape[0]
hWSZ = (n-1)//2
sums = unif2d(img.astype(float)**2,size=(n,n))[hWSZ:-hWSZ,hWSZ:-hWSZ]
out = sums + (tmpl**2).mean() -2*conv2(img,tmpl[::-1,::-1],'valid')/float(n*n)
return out
基于Skimage的NCC方法
使用Skimage计算标准化互相关(Normalized Cross-Correlation)的方法:
def skimage_match_template(img, tmpl):
return match_template(img, tmpl)
请注意,由于这些是交叉相关值,因此图像和模板之间的接近程度将以高输出值来表征。
使用基于OpenCV的方法进行各种相似度测量
OpenCV提供了各种分类模板的方法 -template-matching
。
def opencv_generic(img, tmpl, method_string ='SQDIFF'):
method = eval('cv2.TM_' + method_string)
return cv2m(img.astype('uint8'),tmpl.astype('uint8'),method)
E. 基于视图的方法:自定义函数
我们可以像本文前面所示使用 4D 视图,并沿着最后两个轴执行自定义相似度测量作为 NumPy ufuncs。 因此,我们得到了滑动窗口作为先前使用的 4D 数组,如下所示 -
img_4D = view_as_windows(img, tmpl.shape)
请注意,作为对输入图像的视图,它不会在内存上造成额外的开销。但是,后续的操作会根据这些操作本身进行复制。比较操作将导致更小的内存占用(在Linux上减少8倍)。
然后,我们在
img_4D
和
tmpl
之间执行预期的操作,这在线性映射操作中将导致另一个4D数组,遵循
broadcasting
。我们称之为
img_sub
。接下来,最可能的情况是我们会进行一些减少操作,以给出一个
2D
输出。同样,在大多数情况下,可以使用其中一个
NumPy ufunc
。我们需要在
img_sub
的最后两个轴上使用此ufunc。同样,许多ufunc允许我们一次处理多个轴。例如,之前我们一次性沿着最后两个轴使用了
mean
计算。否则,我们需要依次沿着这两个轴进行操作。
示例
作为如何使用的示例,考虑以下自定义函数:
mean((img_W**tmpl)*tmpl - 2*img*tmpl**2)
这里,我们有 img_W
作为从 img
滑动的窗口,tmpl
通常是在 img
的高度和宽度上滑动的模板。
通过两个嵌套循环实现,我们会得到:
def func1(a,b):
m1,n1 = a.shape
m2,n2 = b.shape
mo,no = m1-m2+1, n1-n2+1
out = np.empty((mo,no))
for i in range(mo):
for j in range(no):
out[i,j] = ((a[i:i+m2,j:j+n2]**2)*b - 2*a[i:i+m2,j:j+n2]*(b**2)).mean()
return out
现在,使用 view_as_windows
,我们可以得到一个矢量化的解决方案:
def func2(a,b):
a4D = view_as_windows(img, tmpl.shape)
return ((a4D**2)*b - 2*a4D*(b**2)).mean(axis=(2,3))
运行时测试 -
In [89]:
...: a = np.random.randint(4,9,(50,100))
...: b = np.random.randint(2,9,(15,15))
...:
In [90]: %timeit func1(a,b)
1 loops, best of 3: 147 ms per loop
In [91]: %timeit func2(a,b)
100 loops, best of 3: 17.8 ms per loop
基准测试: 均方偏差
适量大小的数据集:
In [94]:
...: img = np.random.randint(0,255,(50,100))
...: tmpl = np.random.randint(0,255,(15,15))
...:
In [95]: out1 = skimage_views_MSD_v1(img, tmpl)
...: out2 = skimage_views_MSD_v2(img, tmpl)
...: out3 = convolution_MSD(img, tmpl)
...: out4 = uniform_filter_MSD(img, tmpl)
...: out5 = opencv_generic(img, tmpl, 'SQDIFF')/tmpl.size
...:
...: print np.allclose(out1, out2)
...: print np.allclose(out1, out3)
...: print np.allclose(out1, out4)
...: print np.allclose(out1, out5)
...:
True
True
True
True
In [96]: %timeit skimage_views_MSD_v1(img, tmpl)
...: %timeit skimage_views_MSD_v2(img, tmpl)
...: %timeit convolution_MSD(img, tmpl)
...: %timeit uniform_filter_MSD(img, tmpl)
...: %timeit opencv_generic(img, tmpl, 'SQDIFF')/tmpl.size
...:
100 loops, best of 3: 8.49 ms per loop
100 loops, best of 3: 3.87 ms per loop
100 loops, best of 3: 5.96 ms per loop
100 loops, best of 3: 3.25 ms per loop
10000 loops, best of 3: 201 µs per loop
对于更大的数据集,取决于系统可用的RAM,我们可能需要使用其他方法而不是视图
方法,该方法会留下明显的内存印记。
让我们在剩下的方法上测试更大的数据大小 -
In [97]:
...: img = np.random.randint(0,255,(500,1000))
...: tmpl = np.random.randint(0,255,(15,15))
...:
In [98]: %timeit convolution_MSD(img, tmpl)
...: %timeit uniform_filter_MSD(img, tmpl)
...: %timeit opencv_generic(img, tmpl, 'SQDIFF')/tmpl.size
...:
1 loops, best of 3: 910 ms per loop
1 loops, best of 3: 483 ms per loop
100 loops, best of 3: 16.1 ms per loop
总结
当使用已知相似度量方法时,即使用OpenCV模板匹配方法列出的六种方法之一,并且我们可以访问OpenCV,那么这将是最佳选择。
如果没有OpenCV,对于像均方差这样的特殊情况,我们可以利用卷积来得到相当有效的方法。
对于通用/自定义函数和具有相当大的数据大小,我们可以查看4D
视图以获得高效的向量化解决方案。对于大型数据大小,我们可能需要使用一个循环并使用3D
视图而不是4D
,以减少内存占用。对于非常大的数据大小,您可能必须退回到两个嵌套循环。
如何检查一个二维NumPy数组是否包含特定的值模式?
. - Divakarskimage.util.view_as_windows
获得一个3或4维数组(取决于是否拆分了RBG值),如Divakar的解决方案,然后广播您的比较函数。或者您可以循环(对于高效的搜索功能可能会更快,但这完全取决于相似性度量,因此不能推广)。 - Daniel F