Скользящая статистика по Image Python

Причина, по которой uniform_filter() намного быстрее, чем generic_filter(), связана с Python — для generic_filter() Python вызывается для каждого пикселя, а для uniform_filter() все изображение обрабатывается в собственном коде. (Я обнаружил, что boxFilter() OpenCV даже быстрее, чем uniform_filter(), см. мой ответ на вопрос "дисперсии окна".)

В оставшейся части этого ответа я покажу, как выполнить вычисление перекоса с использованием uniform_filter(), что значительно ускоряет версию на основе generic_filter(), например:

import scipy.ndimage as ndimage, scipy.stats as st
ndimage.generic_filter(img, st.skew, size=(1,5))

st.skew() SciPy (см., например, v0.17.0) вычисляет перекос как

m3 / m2**1.5

где m3 = E[(X-m)**3] (третий центральный момент), m2 = E[(X-m)**2] (дисперсия) и m = E[X] (среднее значение).

Чтобы использовать uniform_filter(), нужно записать это в терминах необработанных моментов, таких как m3p = E[X**3] и m2p = E[X**2] (обычно используется штрих, чтобы отличить необработанный момент от центрального):

m3 = E[(X-m)**3] = ... = m3p - 3*m*m2p + 2*m**3
m2 = E[(X-m)**2] = ... = m2p - m*m

(В случае, если мой "..." пропускает слишком много, этот ответ имеет полный вывод для m2.) Тогда можно реализовать skew() с помощью uniform_filter() (или boxFilter() для некоторого дополнительного ускорения):

def winSkew(img, wsize):
  imgS = img*img
  m, m2p, m3p = (ndimage.uniform_filter(x, wsize) for x in (img, imgS, imgS*img))
  mS = m*m
  return (m3p-3*m*m2p+2*mS*m)/(m2p-mS)**1.5

По сравнению с generic_filter(), winSkew() обеспечивает 654-кратное ускорение в следующем примере на моей машине:

In [185]: img = np.random.randint(0, 256, (500,500)).astype(np.float)

In [186]: %timeit ndimage.generic_filter(img, st.skew, size=(1,5))
1 loops, best of 3: 14.2 s per loop

In [188]: %timeit winSkew(img, (1,5))
10 loops, best of 3: 21.7 ms per loop

И два расчета дают практически идентичные результаты:

In [190]: np.allclose(winSkew(img, (1,5)), ndimage.generic_filter(img, st.skew, size=(1,5)))
Out[190]: True

Таким же образом можно получить код для расчета эксцесса.

Проблема в том, что generic_filter() не может предположить, что ваш фильтр разделим по осям x или y. Таким образом, он должен работать как настоящий 2D-фильтр, а не как последовательность из двух 1D-фильтров, поэтому время выполнения будет намного медленнее.

Средний фильтр и эквивалентен (я думаю) uniform_filter(), который, если вы читали документацию, реализован в виде серии из двух одномерных однородных фильтров.

Я сравнил время с помощью этого блока кода:

import numpy as np
from scipy import ndimage as ndi
from scipy import misc
baboonfile = '/Users/curt/Downloads/BaboonRGB.jpg'  #local download of http://read.pudn.com/downloads169/sourcecode/graph/texture_mapping/776733/Picture/BaboonRGB__.jpg
im = misc.imread(baboonfile)

meanfilt2D = ndi.generic_filter(im, np.mean, size=[3, 3, 1])
%timeit meanfilt2D = ndi.generic_filter(im, np.mean, size=[3, 3, 1])
print meanfilt2D.shape

meanfiltU = ndi.uniform_filter(im, size=[3, 3, 1])
%timeit meanfiltU = ndi.uniform_filter(im, size=[3, 3, 1])
print meanfiltU.shape

Вывод этого блока был:

1 loops, best of 3: 5.22 s per loop
(512, 512, 3)
100 loops, best of 3: 11.8 ms per loop
(512, 512, 3)

так что настоящий двухмерный generic_filter() занимает 5 секунд для маленького изображения, а двухпроходный 1D uniform_filter() занимает всего миллисекунды. (Примечание: разностное изображение meanfilt2D-meanfiltU не было тождественно нулю, но максимальный элемент был равен 2; я думаю, что различия вызваны округлением и неточным типом данных (uint8), используемым для im.)

Для дисперсии и других фильтров вы должны увидеть этот старый пост о переполнении стека, который отвечает на очень связанный вопрос .