Использование шагов для эффективного фильтра скользящих средних

Question

Недавно я узнал об успехах в ответе на этот пост, и мне было интересно, как я мог бы использовать их для вычисления фильтра скользящего среднего более эффективно, чем то, что я предложил в этом посте (используя фильтры свертки).

Это то, что я до сих пор. Он просматривает исходный массив, затем скручивает его на необходимое количество и суммирует значения ядра, чтобы вычислить среднее значение. Я знаю, что края не обрабатываются правильно, но я могу позаботиться об этом позже... Есть ли лучший и более быстрый способ? Цель состоит в том, чтобы отфильтровать большие массивы с плавающей запятой размером до 5000x5000 x 16 слоев, задача, которая scipy.ndimage.filters.convolve довольно медленно

Обратите внимание, что я ищу подключение с 8 соседями, то есть фильтр 3x3 занимает в среднем 9 пикселей (8 вокруг фокусного пикселя) и присваивает это значение пикселю в новом изображении.

import numpy, scipy

filtsize = 3
a = numpy.arange(100).reshape((10,10))
b = numpy.lib.stride_tricks.as_strided(a, shape=(a.size,filtsize), strides=(a.itemsize, a.itemsize))
for i in range(0, filtsize-1):
    if i > 0:
        b += numpy.roll(b, -(pow(filtsize,2)+1)*i, 0)
filtered = (numpy.sum(b, 1) / pow(filtsize,2)).reshape((a.shape[0],a.shape[1]))
scipy.misc.imsave("average.jpg", filtered)

РЕДАКТИРОВАТЬ разъяснение о том, как я вижу это работает:

Текущий код:

1. используйте stride_tricks для генерации массива вроде [[0,1,2],[1,2,3],[2,3,4]...], который соответствует верхней строке ядра фильтра.
2. Прокрутите вдоль вертикальной оси, чтобы получить средний ряд ядра [[10,11,12],[11,12,13],[13,14,15]...] и добавить его в массив, который я получил в 1)
3. Повторите, чтобы получить нижний ряд ядра [[20,21,22],[21,22,23],[22,23,24]...]. На этом этапе я беру сумму каждой строки и делю ее на количество элементов в фильтре, давая мне среднее значение для каждого пикселя (сдвинутое на 1 строку и 1 столбец и с некоторыми странностями по краям, но я могу позаботься об этом позже).

Я надеялся на лучшее использование stride_tricks для непосредственного получения 9 значений или суммы элементов ядра для всего массива, или чтобы кто-то смог убедить меня в другом более эффективном методе...

Answer 1

Что бы это ни стоило, вот как вы можете сделать это, используя "причудливые" уловки. Я собирался опубликовать это вчера, но отвлекся на реальную работу!:)

У @Paul и @eat есть хорошие реализации, использующие различные другие способы сделать это. Просто чтобы продолжить вещи из предыдущего вопроса, я решил опубликовать N-мерный эквивалент.

Вы не сможете значительно победить scipy.ndimage функции для>1D массивов, однако. (scipy.ndimage.uniform_filter должен бить scipy.ndimage.convolve, хоть)

Более того, если вы пытаетесь получить многомерное движущееся окно, вы рискуете взорвать использование памяти всякий раз, когда случайно сделаете копию массива. Хотя исходный "скользящий" массив - это просто просмотр памяти вашего исходного массива, любые промежуточные шаги, которые копируют массив, сделают копию, которая на порядки больше, чем ваш исходный массив (т. Е. Предположим, что вы работаете с исходный массив 100x100... Вид на него (для фильтра размером (3,3)) будет 98x98x3x3, но будет использовать ту же память, что и оригинал. Однако любые копии будут использовать объем памяти, который заполнен массивом 98x98x3x3 было бы!!)

По сути, использование безумных пошаговых трюков отлично подходит для случаев, когда вы хотите векторизовать операции с движущимся окном на одной оси ndarray. Это позволяет очень легко вычислять такие вещи, как движущееся стандартное отклонение и т. Д. С минимальными издержками. Когда вы хотите начать делать это по нескольким осям, это возможно, но обычно вам лучше использовать более специализированные функции. (Такие как scipy.ndimage, так далее)

Во всяком случае, вот как вы это делаете:

import numpy as np

def rolling_window_lastaxis(a, window):
    """Directly taken from Erik Rigtorp's post to numpy-discussion.
    <http://www.mail-archive.com/[e-mail]/msg29450.html>"""
    if window < 1:
       raise ValueError, "`window` must be at least 1."
    if window > a.shape[-1]:
       raise ValueError, "`window` is too long."
    shape = a.shape[:-1] + (a.shape[-1] - window + 1, window)
    strides = a.strides + (a.strides[-1],)
    return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)

def rolling_window(a, window):
    if not hasattr(window, '__iter__'):
        return rolling_window_lastaxis(a, window)
    for i, win in enumerate(window):
        if win > 1:
            a = a.swapaxes(i, -1)
            a = rolling_window_lastaxis(a, win)
            a = a.swapaxes(-2, i)
    return a

filtsize = (3, 3)
a = np.zeros((10,10), dtype=np.float)
a[5:7,5] = 1

b = rolling_window(a, filtsize)
blurred = b.mean(axis=-1).mean(axis=-1)

Итак, что мы получаем, когда делаем b = rolling_window(a, filtsize) является массивом 8x8x3x3, это фактически представление той же памяти, что и исходный массив 10x10. Мы могли бы так же легко использовать фильтры разных размеров по разным осям или работать только по выбранным осям N-мерного массива (т.е. filtsize = (0,3,0,3) на 4-мерном массиве даст нам 6-мерное представление).

Затем мы можем несколько раз применить произвольную функцию к последней оси, чтобы эффективно вычислять вещи в движущемся окне.

Однако, поскольку мы храним временные массивы, которые намного больше, чем наш исходный массив на каждом шаге mean (или же std или что угодно), это совсем не эффективно памяти! Это также не будет очень быстро, либо.

Эквивалент для ndimage просто:

blurred = scipy.ndimage.uniform_filter(a, filtsize, output=a)

Это будет обрабатывать различные граничные условия, выполнять "размывание" на месте, не требуя временной копии массива, и будет очень быстрым. Уловки шага - хороший способ применить функцию к движущемуся окну вдоль одной оси, но они не являются хорошим способом сделать это вдоль нескольких осей, обычно....

Просто мои 0,02 доллара, во всяком случае...

Answer 2

Я недостаточно знаком с Python, чтобы написать код для этого, но два лучших способа ускорить свертки - это либо разделить фильтр, либо использовать преобразование Фурье.

Разделенный фильтр: свертка - это O(M*N), где M и N - количество пикселей в изображении и фильтре соответственно. Поскольку средняя фильтрация с ядром 3 на 3 эквивалентна фильтрации сначала с ядром 3 на 1, а затем с ядром 1 на 3, вы можете получить (3+3)/(3*3) = ~30% улучшение скорости за счет последовательной свертки с двумя 1-d ядрами (это, очевидно, становится лучше, когда ядро ​​становится больше). Конечно, вы все еще можете использовать трюки с шагами.

Преобразование Фурье: conv(A,B) эквивалентно ifft(fft(A)*fft(B))то есть свертка в прямом пространстве становится умножением в пространстве Фурье, где A ваше изображение и B это ваш фильтр Поскольку (поэлементное) умножение преобразований Фурье требует, чтобы A и B имели одинаковый размер, B представляет собой массив size(A) с вашим ядром в самом центре изображения и нулями повсюду. Чтобы разместить ядро ​​3 на 3 в центре массива, вам может понадобиться A в нечетный размер. В зависимости от вашей реализации преобразования Фурье, это может быть намного быстрее, чем свертка (и если вы применяете один и тот же фильтр несколько раз, вы можете предварительно вычислить fft(B), экономя еще 30% вычислительного времени).

Answer 3

Посмотрим:

Это не совсем понятно из вашего вопроса, но я предполагаю, что вы хотите значительно улучшить усреднение такого рода.

import numpy as np
from numpy.lib import stride_tricks as st

def mf(A, k_shape= (3, 3)):
    m= A.shape[0]- 2
    n= A.shape[1]- 2
    strides= A.strides+ A.strides
    new_shape= (m, n, k_shape[0], k_shape[1])
    A= st.as_strided(A, shape= new_shape, strides= strides)
    return np.sum(np.sum(A, -1), -1)/ np.prod(k_shape)

if __name__ == '__main__':
    A= np.arange(100).reshape((10, 10))
    print mf(A)

Теперь, какого улучшения производительности вы бы ожидали?

Обновить:
Прежде всего, предупреждение: код в его текущем состоянии не адаптируется должным образом к форме "ядра". Однако это не моя главная задача сейчас (во всяком случае, идея уже есть, как правильно адаптироваться).

Я только что выбрал новую форму 4D A интуитивно, для меня действительно имеет смысл подумать о двумерном "ядре", центр которого будет центрирован в каждой позиции сетки исходного 2D A.

Но это 4D формирование не может быть "лучшим". Я думаю, что настоящая проблема здесь заключается в производительности суммирования. Нужно уметь находить порядок "лучшего порядка" (из 4D A), чтобы полностью использовать архитектуру кэширования ваших машин. Однако этот порядок может не совпадать с "маленькими" массивами, которые "взаимодействуют" с кэшем вашей машины, и теми, которые не работают (по крайней мере, не так просто).

Обновление 2:
Вот слегка измененная версия mf, Очевидно, что сначала лучше изменить форму на 3D-массив, а затем вместо суммирования просто делать точечное произведение (это имеет преимущество, так что ядро ​​может быть произвольным). Однако это все еще в 3 раза медленнее (на моей машине), чем обновленная функция Паулса.

def mf(A):
    k_shape= (3, 3)
    k= np.prod(k_shape)
    m= A.shape[0]- 2
    n= A.shape[1]- 2
    strides= A.strides* 2
    new_shape= (m, n)+ k_shape
    A= st.as_strided(A, shape= new_shape, strides= strides)
    w= np.ones(k)/ k
    return np.dot(A.reshape((m, n, -1)), w)

Answer 4

Я уверен, что нужно исправить одну вещь - это ваш массив представлений b,

В нем есть несколько элементов из нераспределенной памяти, поэтому вы получите сбои.

Учитывая ваше новое описание вашего алгоритма, первое, что нужно исправить, это то, что вы шагаете за пределы распределения a:

bshape = (a.size-filtsize+1, filtsize)
bstrides = (a.itemsize, a.itemsize)
b = numpy.lib.stride_tricks.as_strided(a, shape=bshape, strides=bstrides)

---

Обновить

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

A = numpy.arange(100).reshape((10,10))

shifts = [(-1,-1),(-1,0),(-1,1),(0,-1),(0,1),(1,-1),(1,0),(1,1)]
B = A[1:-1, 1:-1].copy()
for dx,dy in shifts:
    xstop = -1+dx or None
    ystop = -1+dy or None
    B += A[1+dx:xstop, 1+dy:ystop]
B /= 9

... что кажется простым подходом. Единственная посторонняя операция состоит в том, что она имеет B только однажды. Все сложение, деление и индексация должны выполняться независимо. Если вы делаете 16 групп, вам все еще нужно выделить B один раз, если вы хотите сохранить изображение. Даже если это не поможет, это может прояснить, почему я не понимаю проблему, или, по крайней мере, служить ориентиром для определения времени ускорения других методов. Это выполняется за 2,6 с на моем ноутбуке на массиве float64 5 x 5 000, 0,5 из которых - создание B