计算均方，绝对偏差和自定义相似性度量-Python/NumPy [英] Compute mean squared, absolute deviation and custom similarity measure - Python/NumPy
问题描述
我有一个较大的图像作为2D数组(假设它是500 x 1000像素的灰度图像).我有一个小图像(假设它是15 x 15像素).我想将小图像滑到大图像上，并针对小图像的给定位置，计算小图像与大图像下半部分之间的相似度.
I have a large image as an 2D array (let's assume that it is a 500 by 1000 pixels gray scale image). And I have one small image (let's say is it 15 by 15 pixels). I would like to slide the small image over the large one and for a given position of the small image I would like to calculate a measure of similarity between the small image and the underling part of the big image.
我想灵活选择相似度.例如，我可能想计算均方偏差或均值绝对偏差或其他某种结果(只是某些操作需要两个大小相同的矩阵并返回一个实数).
结果应为2D数组.我想高效地执行此操作(这意味着我想避免循环).
The result should be a 2D array. I want to do this operation efficiently (which means that I want to avoid loops).
为了衡量相似度，我计划使用两个图像颜色之间的平方偏差.但是，正如我已经提到的，如果可以更改度量(例如使用颜色之间的相关性)，那将是很好的.
As a measure of similarity I plan to use the square deviations between the colors of the two images. However, as I have mentioned, it would be nice if I can change the measure (for example to use correlation between the colors).
在numpy中有一个函数可以做到吗?
Is there a function in numpy that can do it?
推荐答案
导入模块
首先，让我们导入将在本文中列出的各种方法中使用的所有相关模块/功能-
Import modules
To start off, let's import all the relevant modules/functions that would be used across the various approaches listed in this 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
A. MAD，MSD基于视图的方法
用于计算mean absolute deviation的基于图像的方法:
A. Views based approach for MAD, MSD
Skimage based approaches for computing mean absolute deviation :
使用scikit-image获取 np.mean 平均计算-
Using scikit-image for getting slided 4D array of views and then np.mean for the average calculations -
def skimage_views_MAD_v1(img, tmpl):
return np.abs(view_as_windows(img, tmpl.shape) - tmpl).mean(axis=(2,3))
使用scikit图像获取4D滑动视图数组，然后 np.einsum 进行平方平均值计算-
Using scikit-image for getting slided 4D array of views and then np.einsum for the squared-average calculations -
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)
用于计算mean squared deviation的基于图像的方法:
Skimage based approaches for computing mean squared deviation :
使用类似的技术，对于mean squared deviation-
Using the similar techniques, we would have two approaches for 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)
B.基于卷积的MSD方法
基于卷积的计算方法mean squared deviations:
B. Convolution based approach for MSD
Convolution based approaches for computing mean squared deviations :
Convolution可用于以一些调整的方式来计算均方差.对于偏差平方和的情况，我们在每个窗口内执行逐元素减法，然后对每个减法求平方，然后对所有这些求和.
Convolution could be used to used to compute mean squared deviations in a bit of tweaked way. For the case of sum of squared deviations, within each window we are performing elementwise subtractions and then squaring each subtraction and then summing up all those.
让我们考虑一下一维示例-
Let's consider a 1D example for a closer look -
a : [a1, a2, a3, a4, a5, a6, a7, a8] # Image array
b : [b1, b2, b3] # Template array
对于第一个窗口操作，我们将:
For the first window operation, we would have :
(a1-b1)**2 + (a2-b2)**2 + (a3-b3)**2
让我们使用(a-b)**2公式:
(a - b)**2 = a**2 - 2*a*b +b**2
因此，我们将拥有第一个窗口:
Thus, we would have for the first window :
(a1**2 - 2*a1*b1 +b1**2) + (a2**2 - 2*a2*b2 +b2**2) + (a3**2 - 2*a3*b3 +b3**2)
类似地，对于第二个窗口:
Similarly, for the second window :
(a2**2 - 2*a2*b1 +b1**2) + (a3**2 - 2*a3*b2 +b2**2) + (a4**2 - 2*a4*b3 +b3**2)
以此类推.
因此，我们将这些计算分为三部分-
So, we have three parts to these computations -
-
a的平方和在滑动窗口中的总和.
Squaring of a's and summing of those in sliding windows.
b的平方并求和.在所有窗口中保持不变.
Squaring of b's and summing of those. This stays the same among all windows.
最后一个难题是:(2*a1*b1, 2*a2*b2, 2*a3*b3)，(2*a2*b1, 2*a3*b2, 2*a4*b3)等等，对于第一个，第二个，依此类推.根据2D与a的卷积和b的翻转版本来计算. > convolution 可以从另一个方向以滑动的方式运行内核，并计算逐元素的乘法运算，并对每个窗口内的元素求和，因此需要此处的翻转.
Final piece of puzzle is : (2*a1*b1, 2*a2*b2, 2*a3*b3), (2*a2*b1, 2*a3*b2, 2*a4*b3) and so on for the first, second and so on windows. This could be computed by a 2D convolution with a and flipped version of b as per the definition of convolution that runs the kernel from the other direction in a sliding and computes elementwise multiplication and sums elements within each window and hence the flipped needed here.
将这些想法扩展到2D情况，并使用Scipy's convolve2d和uniform_filter，我们将有另外两种方法来计算mean squared deviations，就像这样-
Extending these ideas to 2D case and using Scipy's convolve2d and uniform_filter, we would have two more approaches to compute mean squared deviations, like so -
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.基于Skimage的NCC方法
基于Skimage的计算方法normalized cross-correlation:
C. Skimage based approach for NCC
Skimage based approaches for computing normalized cross-correlation :
def skimage_match_template(img, tmpl):
return match_template(img, tmpl)
请注意，由于这些是互相关值，因此图像和模板之间的紧密度将以较高的输出值为特征.
Please note that, since these are cross-correlation values, the closeness between image and template would be characterized by a high output value.
OpenCV提供了多种 template-matching 进行模板分类的方法-
OpenCV offers various template-matching methods of classifying templates -
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视图，并沿最后两个轴执行自定义相似性度量，如NumPy ufuncs.
E. Views based approach : Custom function
We could use 4D views as shown earlier in this post and perform the custom similarity measures as NumPy ufuncs along the last two axes.
因此，我们将滑动窗口作为以前使用的4D数组获得，就像这样-
Thus, we get the sliding windows as a 4D array as used previously, like so -
img_4D = view_as_windows(img, tmpl.shape)
请注意，作为输入图像的视图，它不再占用内存.但是后面的操作将根据这些操作本身进行复制.比较操作将减少内存占用量(确切地说，在Linux上要少8倍).
Please note that being a view into the input image, it won't cost anymore on memory. But the later operations would make copy depending on those operations themselves. A comparison operation would result in much lesser memory footprint (8 times lesser on Linux to be exact).
然后，我们在img_4D和tmpl之间执行预期的操作，这在线性映射的操作中将导致在 NumPy ufuncs 之一这里.我们需要在img_sub的最后两个轴上使用此ufunc.同样，许多功能可以让我们一次在多个轴上工作.例如，前面我们一次沿最后两个轴使用了mean计算.否则，我们需要一个接一个地沿着这两个轴进行操作.
Then, we perform the intended operation between img_4D and tmpl, which in a linearly mapped operation would result in another 4D array following 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 could be used here. We need to use this ufunc along the last two axes on img_sub. Again, many ufuncs allow us to work on more than one axis at a time. For example, earlier we used mean computation along last two axes in one go. Otherwise, we need to operate along those two axes one after another.
示例
作为有关如何使用的示例，让我们考虑一个自定义函数:
As an example on how to use let's consider a custom function :
mean((img_W**tmpl)*tmpl - 2*img*tmpl**2)
在这里，我们以img_W作为img的滑动窗口，而tmpl通常是在img的高度和宽度上滑动的模板.
Here, we have img_W as the sliding window from img and tmpl as usual is the template that is slided across the height and width of img.
通过两个嵌套循环，我们将:
Implemented with two nested loops, we would have :
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，我们将获得向量化的解决方案:
Now, using view_as_windows, we would have a vectorized solution :
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
基准化:均方差
大小合适的数据集:
Benchmarking : Mean squared deviation
Decent sized datasets :
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方法之外的其他方法，这些方法会留下明显的内存占用.
For bigger datasizes, depending on the system RAM available, we might have to fall back to methods other than views method that leaves noticeable memory footprint.
让我们用其他方法测试更大的数据大小-
Let's test out on bigger datasizes with the remaining approaches -
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，那将是最好的方法.
Summary
When working with known similarity measures , i.e. one of the six methods listed with OpenCV based template matching method and if we have access to OpenCV, that would be the best one.
在没有OpenCV的情况下，对于均方差之类的特殊情况，我们可以利用卷积来获得相当有效的方法.
Without OpenCV, for a special case like mean squared deviation, we could make use of convolution to have decently efficient approaches.
对于通用/自定义函数以及适当大小的数据大小，我们可以研究
4D视图以获取有效的矢量化解决方案.对于大数据量，我们可能希望使用一个循环并使用3D视图而不是4D来减少内存占用.对于很大的循环，您可能不得不退回到两个嵌套循环.For generic/custom functions and with decent sized datasizes, we could look into
4Dviews to have efficient vectorized solutions. For large datasizes, we might want to use one loop and use3Dviews instead of4Dwith the intention of reducing memory footprint. For really large ones, you might have to fall back to two nested loops.这篇关于计算均方，绝对偏差和自定义相似性度量-Python/NumPy的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
