Najszybszy splot 2D lub filtr obrazu w języku Python

Question

Kilku użytkowników pytało o szybkość lub zużycie pamięci w nawijaczach obrazów w trybie numpy lub scipy [1, 2, 3, 4]. Z odpowiedzi i mojego doświadczenia z Numpy, wierzę, że może to być poważnym brakiem numpy w porównaniu do Matlab lub IDL.

Żadna z dotychczasowych odpowiedzi nie odniosła się do ogólnego pytania, więc oto jest: "Jaka jest najszybsza metoda obliczania splotów 2D w Pythonie?" Wspólne moduły python to uczciwa gra: numpy, scipy i PIL (inne?). W celu przeprowadzenia trudnego porównania chciałbym zaproponować następujące zasady:

1. Macierze wejściowe mają odpowiednio 2048 x 2048 i 32 x 32.
2. Zmienność pojedyncza lub podwójna zmiennoprzecinkowa są dopuszczalne.
3. Czas spędzony na przekształceniu macierzy wejściowej w odpowiedni format nie jest zliczany - wystarczy krok splotu.
4. Wymiana matrycy wejściowy ze swoim wyjściem jest dopuszczalna (czy jakiekolwiek wsparcie biblioteka Pythona, który?)
5. bezpośrednie DLL wzywa do wspólnych bibliotek C są w porządku - Lapack lub scalapack
6. PyCUDA jest prosto. Używanie niestandardowego sprzętu GPU nie jest sprawiedliwe.

Answer 1

To naprawdę zależy od tego, co chcesz zrobić ... Wiele czasu, nie potrzebujesz w pełni generycznego (czytaj: wolniej) splot 2D ... (tzn. Jeśli filtr jest rozdzielny, używać dwóch 1D zwoje zamiast ... to dlaczego różne scipy.ndimage.gaussian, scipy.ndimage.uniform, są znacznie szybsze niż samo zaimplementowany jako generycznych zwojów nD)

w każdym razie, jako punkt odniesienia.

t = timeit.timeit(stmt='ndimage.convolve(x, y, output=x)', number=1, 
setup=""" 
import numpy as np 
from scipy import ndimage 
x = np.random.random((2048, 2048)).astype(np.float32) 
y = np.random.random((32, 32)).astype(np.float32) 
""") 
print t 

Zajmuje to 6,9 s na moim komputerze ...

Porównaj to z fftconvolve

t = timeit.timeit(stmt="signal.fftconvolve(x, y, mode='same')", number=1, 
setup=""" 
import numpy as np 
from scipy import signal 
x = np.random.random((2048, 2048)).astype(np.float32) 
y = np.random.random((32, 32)).astype(np.float32) 
""") 
print t 

Trwa to około 10,8 sek. Jednak przy różnych wielkościach wejściowych używanie fftów do splatania może być znacznie szybsze (Chociaż nie mogę wymyślić dobrego przykładu, w tej chwili ...).

Answer 2

Na moim komputerze, okrągły splot ręcznie wykonane przy użyciu FFT wydaje się być na czczo:

import numpy 
x = numpy.random.random((2048, 2048)).astype(numpy.float32) 
y = numpy.random.random((32, 32)).astype(numpy.float32) 
z = numpy.fft.irfft2(numpy.fft.rfft2(x) * numpy.fft.rfft2(y, x.shape)) 

Należy pamiętać, że może to traktować obszary blisko krawędzi inaczej niż inne sposoby, bo jest to okrągły splot.

Answer 3

Zrobiłem też kilka eksperymentów z tym. Domyślam się, że splot SciPy nie korzysta z biblioteki BLAS do przyspieszania obliczeń. Używając BLAS, udało mi się zakodować splot 2D, który był porównywalny z szybkością do MATLAB-a. To więcej pracy, ale najlepiej jest przekodować splot w C++.

Tutaj jest ciasna część pętli (proszę wybaczyć odniesienie do tablicy weird(), jest to moja klasa wygody dla macierzy MATLAB) Kluczową częścią jest to, że nie robisz iteracji po obrazie, iterujesz po filtruj i pozwól BLAS na iterację obrazu, ponieważ zazwyczaj obraz jest znacznie większy niż filtr.

for(int n = 0; n < filt.numCols; n++) 
    { 
    for(int m = 0; m < filt.numRows; m++) 
    { 
     const double filt_val = filt(filt.numRows-1-m,filt.numCols-1-n); 
     for (int i =0; i < diffN; i++) 
     { 
     double *out_ptr = &outImage(0,i); 
     const double *im_ptr = &image(m,i+n); 
     cblas_daxpy(diffM,filt_val,im_ptr, 1, out_ptr,1); 

     } 
    } 
}