计算均方,绝对偏差和自定义相似性度量 - Python/NumPy

Question

我有一个大图像作为2D数组(让我们假设它是一个500 x 1000像素的灰度图像).我有一个小图像(比方说是15 x 15像素).我想将小图像滑过大图像,对于小图像的给定位置,我想计算小图像和大图像的下部分之间的相似性度量.

我想灵活地选择一种相似度.例如,我可能想要计算均方偏差或平均绝对偏差或其他东西(只是一些采用相同大小的两个矩阵并返回实数的操作).

结果应该是2D数组.我想有效地执行此操作(这意味着我想避免循环).

作为相似性的度量,我计划使用两个图像的颜色之间的平方偏差.但是,正如我所提到的,如果我可以更改度量(例如使用颜色之间的相关性),那将是很好的.

numpy中有功能吗？

Answer 1

导入模块

首先,让我们导入将在本文中列出的各种方法中使用的所有相关模块/功能 -

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的基于视图的方法

基于Skimage的计算方法mean absolute deviation:

使用scikit-image获取slided 4D array of views然后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)

B.基于卷积的MSD方法

基于卷积的计算方法mean squared deviations:

Convolution可以用来计算一些调整方式的均方偏差.对于平方偏差之和的情况,在每个窗口内我们执行元素减法,然后对每个减法进行平方,然后对所有这些进行求和.

让我们仔细考虑一维示例 -

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)

等等.

所以,我们对这些计算有三个部分 -

• 平滑和滑动窗口中的那些.

• b的平方和那些的总和.所有窗户都保持不变.

• 最后一块拼图是:(2*a1*b1, 2*a2*b2, 2*a3*b3),(2*a2*b1, 2*a3*b2, 2*a4*b3)等等用于第一,第二等的窗户.这可以通过2D卷积a和翻转版本来计算,b根据定义,convolution在滑动中从另一个方向运行内核并计算元素乘法并对每个窗口内的元素求和,因此需要翻转.

扩展这些想法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.基于Skimage的NCC方法

基于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前面所示的视图,并沿着最后两个轴执行自定义相似性度量为NumPy ufuncs.

因此,我们将滑动窗口视为之前使用的4D阵列,如此 -

img_4D = view_as_windows(img, tmpl.shape)

请注意,作为输入图像的视图,它不会再花费在内存上.但是后来的操作会根据这些操作本身进行复制.比较操作会导致更少的内存占用(确切地说是Linux上的8倍).

然后,我们在img_4D和之间执行预期的操作tmpl,在线性映射操作中将导致另一个4D数组跟随broadcasting.我们称之为img_sub.接下来,最有可能的是,我们会进行一些减少操作来为我们提供2D输出.同样,在大多数情况下,其中一个NumPy ufuncs可以在这里使用.我们需要在最后两个轴上使用这个ufunc img_sub.同样,许多ufunc允许我们一次在多个轴上工作.例如,早先我们mean一次性沿着最后两个轴使用计算.否则,我们需要一个接一个地沿着这两个轴操作.

例

作为如何使用let的例子考虑一个自定义函数:

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]: # 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,对于像均方偏差这样的特殊情况,我们可以利用卷积来获得相当有效的方法.

• 对于通用/自定义函数和具有相当大的数据量,我们可以查看4D视图以获得高效的矢量化解决方案.对于大型数据集,我们可能希望使用一个循环并使用3D视图而不是4D为了减少内存占用.对于非常大的,您可能必须回退到两个嵌套循环.