Вычислить средний квадрат, абсолютное отклонение и пользовательскую меру сходства - Python/NumPy

Question

У меня есть большое изображение в виде 2D-массива (предположим, что это изображение размером 500 на 1000 пикселей в оттенках серого). И у меня есть одно маленькое изображение (скажем, 15 на 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

---

А. Взгляд основанный на MAD, MSD

Skimage основанные подходы к вычислениям mean absolute deviation :

Использование scikit-изображения для получения 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-изображения для получения скользящего 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 может быть использован для вычисления средних квадратов отклонений немного измененным способом. Для случая суммы квадратов отклонений в каждом окне мы выполняем поэлементное вычитание, а затем возводим в квадрат каждое вычитание и затем суммируем все эти.

Давайте рассмотрим 1D пример для более детального изучения -

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)

и так далее.

Итак, у нас есть три части к этим вычислениям -

• Возведение в квадрат и суммирование в раздвижных окнах.

• Возведение в квадрат б и суммирование тех. Это остается одинаковым среди всех окон.

• Последний кусок головоломки: (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 представлений, как показано ранее в этом посте, и выполняйте пользовательские измерения сходства как ufuncs NumPy вдоль двух последних осей.

Таким образом, мы получаем скользящие окна в виде 4D-массива, как и раньше, например:

img_4D = view_as_windows(img, tmpl.shape)

Обратите внимание, что просмотр входного изображения не потребует больше памяти. Но более поздние операции сделают копию в зависимости от самих этих операций. Операция сравнения приведет к гораздо меньшему использованию памяти (если быть точным, в 8 раз меньше в Linux).

Затем мы выполняем предполагаемую операцию между img_4D а также tmpl, что в линейно отображенной операции приведет к следующему массиву 4D broadcasting, Давайте назовем это img_sub, Далее, скорее всего, у нас будет некоторая операция сокращения, чтобы дать нам 2D выход. Опять же, в большинстве случаев, один из NumPy ufuncs может быть использован здесь. Нам нужно использовать этот ufunc вдоль двух последних осей на img_sub, Опять же, многие функции позволяют нам работать более чем на одной оси одновременно. Например, ранее мы использовали 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]: # 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

Для больших размеров данных, в зависимости от доступной оперативной памяти системы, нам, возможно, придется прибегнуть к методам, отличным от 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 с намерением уменьшить объем памяти. Для действительно больших, вам, возможно, придется вернуться к двум вложенным циклам.

Answer 2

Вы имеете в виду операцию кросс-корреляции?

Однако, если вы строго хотите проверить сходство с квадратом отклонения, вы можете использовать сопоставление с шаблоном в Skimage, который использует более быструю реализацию взаимной корреляции. Пример здесь: http://scikit-image.org/docs/dev/auto_examples/plot_template.html

В противном случае вы можете использовать correlate2d для достижения этой цели следующим образом: 1. Выполните взаимную корреляцию для сигнала с нулевым средним (то есть оба сигнала / изображения должны быть центрированы относительно нуля) 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

https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.argrelmax.html