Самый быстрый способ выполнить вычисления для каждого подмассива NXN в двумерном массиве

Question

У меня есть двумерный массив, который представляет изображение в градациях серого.Мне нужно извлечь каждый N x N вложенный массив в этом массиве с указанным перекрытием между вложенными массивами и вычислить такое свойство, как среднее значение, стандартное отклонение или медиана.

Приведенный ниже код выполняет эту задачу, но довольно медленно, потому что он использует Python для циклов.Любые идеи о том, как векторизовать этот расчет или иным образом ускорить его?

import numpy as np

img = np.random.randn(100, 100)
N = 4
step = 2

h, w = img.shape
out = []
for i in range(0, h - N, step):
    outr = []
    for j in range(0, w - N, step):
        outr.append(np.mean(img[i:i+N, j:j+N]))
    out.append(outr)
out = np.array(out)

Answer 1

Для среднего и стандартного отклонения существует быстрое решение на основе cumsum.

Здесь приведены временные параметры для изображения 500x200, окна 30x20 и размеров шагов 5 и 3. Для сравнения я использую skimage.util.view_as_windows сNumPy среднее и стандартное.

mn + sd using cumsum     1.1531693299184553 ms
mn using view_as_windows 3.495307120028883 ms
sd using view_as_windows 21.855629019846674 ms

Код:

import numpy as np
from math import gcd
from timeit import timeit

def wsum2d(A, winsz, stepsz, canoverwriteA=False):
    M, N = A.shape
    m, n = winsz
    i, j = stepsz
    for X, x, s in ((M, m, i), (N, n, j)):
        g = gcd(x, s)
        if g > 1:
            X //= g
            x //= g
            s //= g
            A = A[:X*g].reshape(X, g, -1).sum(axis=1)
        elif not canoverwriteA:
            A = A.copy()
        canoverwriteA = True
        A[x:] -= A[:-x]
        A = A.cumsum(axis=0)[x-1::s]
        A = A.T
    return A

def w2dmnsd(A, winsz, stepsz):
    # combine A and A*A into a complex, so overheads apply only once
    M21 = wsum2d(A*(A+1j), winsz, stepsz, True)
    M2, mean_ = M21.real / np.prod(winsz), M21.imag / np.prod(winsz)
    sd = np.sqrt(M2 - mean_*mean_)
    return mean_, sd

# test
np.random.seed(0)
A = np.random.random((500, 200))
wsz = (30, 20)
stpsz = (5, 3)
mn, sd = w2dmnsd(A, wsz, stpsz)
from skimage.util import view_as_windows
Av = view_as_windows(A, wsz, stpsz) # this emits a warning on my system
assert np.allclose(mn, np.mean(Av, axis=(2, 3)))
assert np.allclose(sd, np.std(Av, axis=(2, 3)))
from timeit import repeat

print('mn + sd using cumsum    ', min(repeat(lambda: w2dmnsd(A, wsz, stpsz), number=100))*10, 'ms')
print('mn using view_as_windows', min(repeat(lambda: np.mean(Av, axis=(2, 3)), number=100))*10, 'ms')
print('sd using view_as_windows', min(repeat(lambda: np.std(Av, axis=(2, 3)), number=100))*10, 'ms')

Answer 2

Общий случай может быть решен с использованием scipy.ndimage.generic_filter:

import numpy as np

from scipy.ndimage import generic_filter

img = np.random.randn(100, 100)

N = 4
filtered = generic_filter(img.astype(np.float), np.std, size=N)

step = 2
output = filtered[::step, ::step]

Однако, на самом деле, это может работать не намного быстрее, чем простой цикл for.

Для применения среднего и медианного фильтра вы можете использовать skimage.rank.mean и skimage.rank.median соответственно, что должно быть быстрее.Существует также scipy.ndimage.median_filter.В противном случае среднее значение также может быть эффективно вычислено путем простой свертки с массивом (N, N) со значениями 1./N^2.Для стандартного отклонения вам, вероятно, придется прикусить пулю и использовать generic_filter, если размер вашего шага не больше или равен N.

Answer 3

Если Numba является опцией, единственное, что нужно сделать, это избежать добавления списка (он работает и с добавками в список, но медленнее . Чтобы использовать и параллелизацию, немного переписали реализацию, чтобыизбегайте шага в пределах диапазона, который не поддерживается при использовании parfor.

Пример

@nb.njit(error_model='numpy',parallel=True)
def calc_p(img,N,step):
  h,w=img.shape

  i_w=(h - N)//step
  j_w=(w - N)//step
  out = np.empty((i_w,j_w))
  for i in nb.prange(0, i_w):
      for j in range(0, j_w):
          out[i,j]=np.std(img[i*step:i*step+N, j*step:j*step+N])
  return out

def calc_n(img,N,step):
  h, w = img.shape
  out = []
  for i in range(0, h - N, step):
      outr = []
      for j in range(0, w - N, step):
          outr.append(np.std(img[i:i+N, j:j+N]))
      out.append(outr)
  return(np.array(out))

Сроки

Всевремя без компиляции составляет около 0,5 с (первый вызов функции исключен из времени).

#Data
img = np.random.randn(100, 100)
N = 4
step = 2

calc_n :17ms
calc_p :0.033ms

Поскольку это на самом деле скользящее среднее, есть дополнительные возможности для улучшения, если N получитбольше.

Answer 4

Вы можете использовать scikit-image block_reduce:

Таким образом, ваш код становится:

import numpy as np
import skimage.measure

N = 4

# Your main array
a = np.arange(9).reshape(3,3)

mean = skimage.measure.block_reduce(a, (N,N), np.mean) 
std_dev = skimage.measure.block_reduce(a, (N,N), np.std)
median = skimage.measure.block_reduce(a, (N,N), np.median)

Однако приведенный выше код работает только для шагов / шаговразмером 1.

Для среднего вы можете использовать среднее пул, который доступен в любом современном пакете ML.Что касается срединного и стандартного отклонения, то это правильный подход.