计算均方、绝对偏差和自定义相似性度量 - Python/NumPy
Compute mean squared, absolute deviation and custom similarity measure - Python/NumPy
我有一个二维数组的大图像(假设它是一个 500 x 1000 像素的灰度图像)。我有一张小图片(假设它是 15 x 15 像素)。我想将小图像滑到大图像上,对于小图像的给定位置,我想计算小图像与大图像下层部分之间的相似性度量。
我想灵活地选择相似性度量。例如,我可能想计算均方差或平均绝对差或其他东西(只是一些需要两个相同大小的矩阵和 returns 实数的操作)。
结果应该是一个二维数组。我想高效地执行此操作(这意味着我想避免循环)。
作为相似度的衡量标准,我计划使用两个图像的颜色之间的平方差。但是,正如我所提到的,如果我可以更改度量(例如使用颜色之间的相关性),那就太好了。
numpy中有没有函数可以做到?
您指的是互相关操作吗?
但是,如果您严格地想要检查平方偏差的相似性,您可以使用 skimage 中的模板匹配,它使用更快的互相关实现。这里的例子:http://scikit-image.org/docs/dev/auto_examples/plot_template.html
否则,您可以使用 correlate2d 来实现,如下所示:
1. 对零均值信号执行互相关(意味着 signals/images 应该以零为中心)
2. 检查局部最大值 scipy.signal.argrelmax 或(如果您认为只有一个匹配项)使用 np.argmax
查找全局最大值
这是一个示例(从文档中提取),如果需要,您可以将 np.argmax 替换为 signal.argrelmax
from scipy import signal
from scipy import misc
lena = misc.lena() - misc.lena().mean()
template = np.copy(lena[235:295, 310:370]) # right eye
template -= template.mean()
lena = lena + np.random.randn(*lena.shape) * 50 # add noise
corr = signal.correlate2d(lena, template, boundary='symm', mode='same')
y, x = np.unravel_index(np.argmax(corr), corr.shape) # find the match
来源:
https://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.signal.correlate2d.html
导入模块
首先,让我们导入所有相关的 modules/functions,这些 modules/functions 将用于此 post -
中列出的各种方法
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 的计算方法 mean absolute deviation
:
使用 scikit-image 获取 slided 4D array of views
and then 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 的计算方法 mean squared deviation
:
使用类似的技术,我们将有两种方法 mean squared deviation
-
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)
乙。 MSD
基于卷积的方法
基于卷积的计算方法mean squared deviations
:
Convolution
可用于以稍微调整的方式计算均方差。对于偏差平方和的情况,在每个 window 中,我们执行元素减法,然后对每个减法求平方,然后对所有这些求和。
让我们考虑一个一维的例子来仔细看看 -
a : [a1, a2, a3, a4, a5, a6, a7, a8] # Image array
b : [b1, b2, b3] # Template array
对于第一个 window 操作,我们将有:
(a1-b1)**2 + (a2-b2)**2 + (a3-b3)**2
让我们使用 (a-b)**2
公式:
(a - b)**2 = a**2 - 2*a*b +b**2
因此,第一个 window :
(a1**2 - 2*a1*b1 +b1**2) + (a2**2 - 2*a2*b2 +b2**2) + (a3**2 - 2*a3*b3 +b3**2)
同样,对于第二个 window :
(a2**2 - 2*a2*b1 +b1**2) + (a3**2 - 2*a3*b2 +b2**2) + (a4**2 - 2*a4*b3 +b3**2)
等等。
因此,我们将这些计算分为三个部分 -
a 的平方和滑动 windows.
b 的平方和它们的总和。这在所有 windows.
中保持不变
最后一块拼图是:(2*a1*b1, 2*a2*b2, 2*a3*b3)
,(2*a2*b1, 2*a3*b2, 2*a4*b3)
,依此类推第一块,第二块,依此类推windows。这可以通过 2D
与 a
的卷积和 b
的翻转版本来计算,根据 convolution
的定义,它从滑动的另一个方向运行内核并计算每个 window 中的元素乘法和求和元素,因此这里需要翻转。
将这些想法扩展到 2D
案例并使用 Scipy's convolve2d
和 uniform_filter
,我们将有另外两种方法来计算 mean squared deviations
,就像这样 -
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
C。 NCC 的基于 Skimage 的方法
基于 Skimage 的计算方法 normalized cross-correlation
:
def skimage_match_template(img, tmpl):
return match_template(img, tmpl)
请注意,由于这些是互相关值,图像和模板之间的接近度将以高输出值为特征。
D.用于各种相似性度量的基于 OpenCV 的方法
OpenCV 提供了多种template-matching
模板分类方法 -
def opencv_generic(img, tmpl, method_string ='SQDIFF'):
# Methods :
# 'CCOEFF' : Correlation coefficient
# 'CCOEFF_NORMED' : Correlation coefficient normalized
# 'CCORR' : Cross-correlation
# 'CCORR_NORMED' : Cross-correlation normalized
# 'SQDIFF' : Squared differences
# 'SQDIFF_NORMED' : Squared differences normalized
method = eval('cv2.TM_' + method_string)
return cv2m(img.astype('uint8'),tmpl.astype('uint8'),method)
E.基于视图的方法:自定义函数
我们可以使用 4D
前面显示的视图 post 并沿着最后两个轴执行自定义相似性度量作为 NumPy ufunc。
因此,我们得到滑动 windows 作为一个 4D 数组,如前所述 -
img_4D = view_as_windows(img, tmpl.shape)
请注意,作为输入图像的视图,它不会再占用内存。但是后面的操作会根据这些操作本身进行复制。比较操作将导致更少的内存占用(确切地说 Linux 减少 8 倍)。
然后,我们在 img_4D
和 tmpl
之间执行预期的操作,这在线性映射操作中会导致在 broadcasting
. Let's call it img_sub
. Next up, most probably, we would have some reduction operation to give us a 2D
output. Again, in most cases, one of the NumPy ufuncs
之后的另一个 4D 数组可以在这里使用。我们需要沿着 img_sub
上的最后两个轴使用这个 ufunc。同样,许多 ufunc 允许我们一次处理多个轴。例如,之前我们使用 mean
沿最后两个轴一次性计算。否则,我们需要沿着这两个轴依次操作。
例子
作为如何使用的示例,让我们考虑一个自定义函数:
mean((img_W**tmpl)*tmpl - 2*img*tmpl**2)
这里,我们有img_W
作为img
的滑动window和tmpl
像往常一样是在[=的高度和宽度上滑动的模板60=].
用两个嵌套循环实现,我们将有:
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]: # Sample image(a) and template(b)
...: 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]: # Inputs
...: 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,我们可能不得不回退到 views
方法以外的方法,这会留下明显的内存占用。
让我们用剩下的方法测试更大的数据量 -
In [97]: # Inputs
...: 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,对于像均方差这样的特殊情况,我们可以利用卷积来获得相当有效的方法。
对于 generic/custom 函数和适当大小的数据,我们可以查看 4D
视图以获得高效的矢量化解决方案。对于大数据量,我们可能希望使用一个循环并使用 3D
视图而不是 4D
以减少内存占用。对于非常大的循环,您可能不得不回退到两个嵌套循环。
我有一个二维数组的大图像(假设它是一个 500 x 1000 像素的灰度图像)。我有一张小图片(假设它是 15 x 15 像素)。我想将小图像滑到大图像上,对于小图像的给定位置,我想计算小图像与大图像下层部分之间的相似性度量。
我想灵活地选择相似性度量。例如,我可能想计算均方差或平均绝对差或其他东西(只是一些需要两个相同大小的矩阵和 returns 实数的操作)。
结果应该是一个二维数组。我想高效地执行此操作(这意味着我想避免循环)。
作为相似度的衡量标准,我计划使用两个图像的颜色之间的平方差。但是,正如我所提到的,如果我可以更改度量(例如使用颜色之间的相关性),那就太好了。
numpy中有没有函数可以做到?
您指的是互相关操作吗?
但是,如果您严格地想要检查平方偏差的相似性,您可以使用 skimage 中的模板匹配,它使用更快的互相关实现。这里的例子:http://scikit-image.org/docs/dev/auto_examples/plot_template.html
否则,您可以使用 correlate2d 来实现,如下所示: 1. 对零均值信号执行互相关(意味着 signals/images 应该以零为中心) 2. 检查局部最大值 scipy.signal.argrelmax 或(如果您认为只有一个匹配项)使用 np.argmax
查找全局最大值这是一个示例(从文档中提取),如果需要,您可以将 np.argmax 替换为 signal.argrelmax
from scipy import signal
from scipy import misc
lena = misc.lena() - misc.lena().mean()
template = np.copy(lena[235:295, 310:370]) # right eye
template -= template.mean()
lena = lena + np.random.randn(*lena.shape) * 50 # add noise
corr = signal.correlate2d(lena, template, boundary='symm', mode='same')
y, x = np.unravel_index(np.argmax(corr), corr.shape) # find the match
来源:
https://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.signal.correlate2d.html
导入模块
首先,让我们导入所有相关的 modules/functions,这些 modules/functions 将用于此 post -
中列出的各种方法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 的计算方法 mean absolute deviation
:
使用 scikit-image 获取 slided 4D array of views
and then 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 的计算方法 mean squared deviation
:
使用类似的技术,我们将有两种方法 mean squared deviation
-
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)
乙。 MSD
基于卷积的方法基于卷积的计算方法mean squared deviations
:
Convolution
可用于以稍微调整的方式计算均方差。对于偏差平方和的情况,在每个 window 中,我们执行元素减法,然后对每个减法求平方,然后对所有这些求和。
让我们考虑一个一维的例子来仔细看看 -
a : [a1, a2, a3, a4, a5, a6, a7, a8] # Image array
b : [b1, b2, b3] # Template array
对于第一个 window 操作,我们将有:
(a1-b1)**2 + (a2-b2)**2 + (a3-b3)**2
让我们使用 (a-b)**2
公式:
(a - b)**2 = a**2 - 2*a*b +b**2
因此,第一个 window :
(a1**2 - 2*a1*b1 +b1**2) + (a2**2 - 2*a2*b2 +b2**2) + (a3**2 - 2*a3*b3 +b3**2)
同样,对于第二个 window :
(a2**2 - 2*a2*b1 +b1**2) + (a3**2 - 2*a3*b2 +b2**2) + (a4**2 - 2*a4*b3 +b3**2)
等等。
因此,我们将这些计算分为三个部分 -
a 的平方和滑动 windows.
b 的平方和它们的总和。这在所有 windows.
中保持不变最后一块拼图是:
(2*a1*b1, 2*a2*b2, 2*a3*b3)
,(2*a2*b1, 2*a3*b2, 2*a4*b3)
,依此类推第一块,第二块,依此类推windows。这可以通过2D
与a
的卷积和b
的翻转版本来计算,根据convolution
的定义,它从滑动的另一个方向运行内核并计算每个 window 中的元素乘法和求和元素,因此这里需要翻转。
将这些想法扩展到 2D
案例并使用 Scipy's convolve2d
和 uniform_filter
,我们将有另外两种方法来计算 mean squared deviations
,就像这样 -
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
C。 NCC 的基于 Skimage 的方法
基于 Skimage 的计算方法 normalized cross-correlation
:
def skimage_match_template(img, tmpl):
return match_template(img, tmpl)
请注意,由于这些是互相关值,图像和模板之间的接近度将以高输出值为特征。
D.用于各种相似性度量的基于 OpenCV 的方法
OpenCV 提供了多种template-matching
模板分类方法 -
def opencv_generic(img, tmpl, method_string ='SQDIFF'):
# Methods :
# 'CCOEFF' : Correlation coefficient
# 'CCOEFF_NORMED' : Correlation coefficient normalized
# 'CCORR' : Cross-correlation
# 'CCORR_NORMED' : Cross-correlation normalized
# 'SQDIFF' : Squared differences
# 'SQDIFF_NORMED' : Squared differences normalized
method = eval('cv2.TM_' + method_string)
return cv2m(img.astype('uint8'),tmpl.astype('uint8'),method)
E.基于视图的方法:自定义函数
我们可以使用 4D
前面显示的视图 post 并沿着最后两个轴执行自定义相似性度量作为 NumPy ufunc。
因此,我们得到滑动 windows 作为一个 4D 数组,如前所述 -
img_4D = view_as_windows(img, tmpl.shape)
请注意,作为输入图像的视图,它不会再占用内存。但是后面的操作会根据这些操作本身进行复制。比较操作将导致更少的内存占用(确切地说 Linux 减少 8 倍)。
然后,我们在 img_4D
和 tmpl
之间执行预期的操作,这在线性映射操作中会导致在 broadcasting
. Let's call it img_sub
. Next up, most probably, we would have some reduction operation to give us a 2D
output. Again, in most cases, one of the NumPy ufuncs
之后的另一个 4D 数组可以在这里使用。我们需要沿着 img_sub
上的最后两个轴使用这个 ufunc。同样,许多 ufunc 允许我们一次处理多个轴。例如,之前我们使用 mean
沿最后两个轴一次性计算。否则,我们需要沿着这两个轴依次操作。
例子
作为如何使用的示例,让我们考虑一个自定义函数:
mean((img_W**tmpl)*tmpl - 2*img*tmpl**2)
这里,我们有img_W
作为img
的滑动window和tmpl
像往常一样是在[=的高度和宽度上滑动的模板60=].
用两个嵌套循环实现,我们将有:
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]: # Sample image(a) and template(b)
...: 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]: # Inputs
...: 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,我们可能不得不回退到 views
方法以外的方法,这会留下明显的内存占用。
让我们用剩下的方法测试更大的数据量 -
In [97]: # Inputs
...: 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,对于像均方差这样的特殊情况,我们可以利用卷积来获得相当有效的方法。
对于 generic/custom 函数和适当大小的数据,我们可以查看
4D
视图以获得高效的矢量化解决方案。对于大数据量,我们可能希望使用一个循环并使用3D
视图而不是4D
以减少内存占用。对于非常大的循环,您可能不得不回退到两个嵌套循环。