Do przodu FFT obraz i do tyłu FFT obraz, aby uzyskać ten sam wynik

Question

Próbuję FFT obrazu przy użyciu biblioteki z http://www.fftw.org/, dzięki czemu można zrobić splot w dziedzinie częstotliwości. Ale nie mogę wymyślić, jak to działa. Aby zrozumieć, jak to zrobić, próbuję przesłać dalej FFT obraz jako tablicę pikseli, a następnie FFT do tyłu, aby uzyskać tę samą tablicę pikseli. Oto, co zrobić:

fftw_plan planR, planG, planB; 
fftw_complex *inR, *inG, *inB, *outR, *outG, *outB, *resultR, *resultG, *resultB; 

//Allocate arrays. 
inR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width); 
inG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width); 
inB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width); 

outR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width); 
outG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width); 
outB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width); 

resultR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width); 
resultG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width); 
resultB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width); 

//Fill in arrays with the pixelcolors. 
for (int y = 0; y < height; y++) { 
    for (int x = 0; x < width; x++) { 
     int currentIndex = ((y * width) + (x)) * 3; 
     inR[y * width + x][0] = pixelColors[currentIndex]; 
     inG[y * width + x][0] = pixelColors[currentIndex + 1]; 
     inB[y * width + x][0] = pixelColors[currentIndex + 2]; 
    } 
} 

//Forward plans. 
planR = fftw_plan_dft_2d(width, width, inR, outR, FFTW_FORWARD, FFTW_MEASURE); 
planG = fftw_plan_dft_2d(width, width, inG, outG, FFTW_FORWARD, FFTW_MEASURE); 
planB = fftw_plan_dft_2d(width, width, inB, outB, FFTW_FORWARD, FFTW_MEASURE); 

//Forward FFT. 
fftw_execute(planR); 
fftw_execute(planG); 
fftw_execute(planB); 

//Backward plans. 
planR = fftw_plan_dft_2d(width, width, outR, resultR, FFTW_BACKWARD, FFTW_MEASURE); 
planG = fftw_plan_dft_2d(width, width, outG, resultG, FFTW_BACKWARD, FFTW_MEASURE); 
planB = fftw_plan_dft_2d(width, width, outB, resultB, FFTW_BACKWARD, FFTW_MEASURE); 

//Backward fft 
fftw_execute(planR); 
fftw_execute(planG); 
fftw_execute(planB); 

//Overwrite the pixelcolors with the result. 
for (int y = 0; y < height; y++) { 
    for (int x = 0; x < width; x++) { 
     int currentIndex = ((y * width) + (x)) * 3; 
     pixelColors[currentIndex] = resultR[y * width + x][0]; 
     pixelColors[currentIndex + 1] = resultG[y * width + x][0]; 
     pixelColors[currentIndex + 2] = resultB[y * width + x][0]; 
    } 
} 

Czy ktoś mógłby mi pokazać przykład sposobu przekazania FFT obraz, a następnie do tyłu FFT obraz, używając FFTW aby uzyskać ten sam wynik? Patrzyłem na wiele przykładów pokazujących, jak używać FFTW do FFT, ale nie mogę zrozumieć, jak to ma zastosowanie w mojej sytuacji, w której mam szereg pikseli reprezentujących obraz.

Answer 1

Jedna ważna rzecz, którą należy zwrócić, kiedy wykonujemy do przodu FFT, a następnie odwrotną FFT, jest to, że zwykle skutkuje to współczynnikiem skalowania N, który jest stosowany do końcowego wyniku, tj. Wynikowe wartości pikseli obrazu muszą być podzielone przez N w aby dopasować oryginalne wartości pikseli. (N to rozmiar FFT). Więc prawdopodobnie pętla wyjście powinno wyglądać mniej więcej tak:

//Overwrite the pixelcolors with the result. 
for (int y = 0; y < height; y++) { 
    for (int x = 0; x < width; x++) { 
     int currentIndex = ((y * width) + (x)) * 3; 
     pixelColors[currentIndex] = resultR[y * width + x][0]/(width * height); 
     pixelColors[currentIndex + 1] = resultG[y * width + x][0]/(width * height); 
     pixelColors[currentIndex + 2] = resultB[y * width + x][0]/(width * height); 
    } 
} 

Należy również pamiętać, że prawdopodobnie chcesz zrobić prawdziwy do kompleksu FFT następuje przez złożony-to- prawdziwe IFFT (nieco bardziej wydajne pod względem pamięci i wydajności). Na razie jednak wygląda na to, że robisz skomplikowane i złożone w obu kierunkach, co jest w porządku, ale nie wypełniasz poprawnie tablic wejściowych. Jeśli masz zamiar trzymać się z kompleksu do kompleksu następnie prawdopodobnie chcesz zmienić pętli wejściowego coś takiego:

//Fill in arrays with the pixelcolors. 
for (int y = 0; y < height; y++) { 
    for (int x = 0; x < width; x++) { 
     int currentIndex = ((y * width) + (x)) * 3; 
     inR[y * width + x][0] = (double)pixelColors[currentIndex]; 
     inR[y * width + x][1] = 0.0; 
     inG[y * width + x][0] = (double)pixelColors[currentIndex + 1]; 
     inG[y * width + x][1] = 0.0; 
     inB[y * width + x][0] = (double)pixelColors[currentIndex + 2]; 
     inB[y * width + x][1] = 0.0; 
    } 
} 

czyli wartości pikseli iść do rzeczywistych części kompleksu wartości wejściowego i Części urojonych należy wyzerować.

Jeszcze jedna rzecz do zapamiętania: kiedy w końcu uda się to zrobić, okaże się, że wydajność jest straszna - przygotowanie planu w stosunku do czasu potrzebnego na faktyczną FFT zajmuje dużo czasu. Chodzi o to, że plan tworzy się tylko raz, ale używa się go do wykonywania wielu FFT. Więc chcesz oddzielić tworzenie planu od rzeczywistego kodu FFT i umieścić go w procedurze inicjalizacji lub konstruktorze lub czymkolwiek innym.

Answer 2

Ale jeśli używasz funkcji realToComplex lub ComplexToRealFunction, zwróć uwagę, że obraz będzie przechowywany w macierzy o wymiarach [wysokość x (szerokość/2 +1)] i jeśli chcesz wykonać obliczenia pośrednie w w dziedzinie częstotliwości, staną się nieco trudniejsze ...