Wykorzystanie kroków do wydajnego filtra średniej ruchomej

Nie znam Pythona na tyle, aby napisać do tego Kod, ale dwa najlepsze sposoby na przyspieszenie splotów to albo oddzielenie filtra, albo użycie transformaty Fouriera.

Separated filter: splot to O (M*N), gdzie M I N to odpowiednio liczba pikseli obrazu i filtra. Ponieważ średnie filtrowanie z jądrem 3 na 3 jest równoważne filtrowaniu najpierw z jądrem 3 na 1, a następnie z jądrem 1 na 3, można uzyskać (3+3)/(3*3) = ~30% poprawę prędkości przez kolejne splot z dwoma jądrami 1-d (to oczywiście staje się lepsze, gdy jądro się powiększa). Oczywiście nadal możesz używać trików z krokami.

transformata Fouriera : conv(A,B) jest równoważne ifft(fft(A)*fft(B)), tzn. splot w przestrzeni bezpośredniej staje się mnożeniem w przestrzeni Fouriera, gdzie A jest twoim obrazem, a {[4] } filtrem. Ponieważ mnożenie transformat Fouriera wymaga, aby a i B były tej samej wielkości, B jest tablicą size(A) z jądrem w samym centrum obrazu i zerami wszędzie indziej. Aby umieścić jądro 3 na 3 na środku tablicy, być może będziesz musiał użyć pad A do nieparzystego rozmiaru. W zależności od implementacji transformaty Fouriera może to być dużo szybsze niż splot (a jeśli zastosujesz ten sam filtr wiele razy, możesz wstępnie obliczyć fft(B), oszczędzając kolejne 30% czasu obliczeń).

Jedna rzecz, którą jestem pewien, że należy naprawić, to Twoja tablica widoków b.

Ma kilka elementów z nieprzydzielonej pamięci, więc będziesz miał awarie.

Biorąc pod uwagę nowy opis algorytmu, pierwszą rzeczą, która wymaga naprawy, jest fakt, że stąpasz poza alokacją a:

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

Update

Ponieważ wciąż nie do końca rozumiem metodę i wydaje się, że są prostsze sposoby rozwiązania problemu, zamierzam umieścić to tutaj:

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

...co wydaje się prostym podejściem. Jedyną operacją zewnętrzną jest to, że przydziela i wypełnia B tylko raz. Wszystkie dodawanie, dzielenie i indeksowanie musi odbywać się niezależnie. Jeśli wykonujesz 16 pasm, musisz przydzielić B tylko raz, jeśli chcesz zapisać obraz. Nawet jeśli to nie pomoże, może to wyjaśnić, dlaczego nie rozumiem problemu, lub przynajmniej służyć jako punkt odniesienia do czasu przyspieszenia innych metod. To działa w 2.6 s na moim laptopie na macierzy 5K x 5K float64, z których 0.5 jest tworzeniem B

Zobaczmy:

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)

Teraz, jakiego rodzaju poprawy wydajności można się spodziewać?

Update:
Przede wszystkim Ostrzeżenie: kod w obecnym stanie nie dostosowuje się właściwie do kształtu 'jądra'. Jednak to nie jest moim głównym zmartwieniem w tej chwili (w każdym razie pomysł jest już gotowy, jak prawidłowo się dostosować).

Właśnie wybrałem nowy kształt 4D a intuicyjnie, dla mnie naprawdę sensowne jest myślenie o centrum jądra 2D, które będzie wyśrodkowane do każdej pozycji siatki oryginalnego 2D A.

Ale to kształtowanie 4D może nie być tak naprawdę "najlepsze". Myślę, że prawdziwym problemem jest tutaj wydajność sumowania. Należy być w stanie znaleźć "najlepszy porządek" (z 4D A) w celu pełnego wykorzystania architektury pamięci podręcznej maszyn. Jednak kolejność ta może nie być taka sama dla "małych" tablic, które rodzaj "współpracuje" z pamięci podręcznej maszyn i tych większych, które nie (przynajmniej nie tak prosty sposób).

Aktualizacja 2:
Oto nieco zmodyfikowana wersja mf. Oczywiście lepiej najpierw przekształcić tablicę 3D, a następnie zamiast sumować po prostu zrób produkt dot (ma to tę zaletę, że jądro może być dowolne). Jednak nadal jest trochę 3x wolniejszy (na moim komputerze) niż zaktualizowana funkcja Paulsa.

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)