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

Ответ 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 post to numpy-discussion.
    <http://www.mail-archive.com/[email protected]/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, во всяком случае...

Ответ 2

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

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

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

Ответ 3

Одна вещь, которую я уверен, должна быть исправлена, - это ваш массив представлений 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 секунды на моем ноутбуке на 5k x 5k массиве float64, из которых 0,5 - создание b

Ответ 4

Давайте посмотрим:

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

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)

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

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

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

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

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

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)