Adelante FFT una imagen y FFT hacia atrás una imagen para obtener el mismo resultado

Question

Estoy tratando de FFT una imagen usando la biblioteca de http://www.fftw.org/ para que pueda hacer una convolución en el dominio de la frecuencia. Pero no puedo encontrar la manera de hacerlo funcionar. Para entender cómo hacer esto, estoy tratando de enviar FFT una imagen como una matriz de pixelcolors y luego hacia atrás FFT para obtener la misma matriz de pixelcolors. Esto es lo que hago:

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]; 
    } 
} 

Podría alguien por favor me muestran un ejemplo de cómo enviar una imagen FFT FFT y luego hacia atrás la imagen utilizando FFTW para obtener el mismo resultado? He estado viendo muchos ejemplos que muestran cómo usar FFTW para FFT, pero no puedo entender cómo se aplica a mi situación donde tengo una matriz de pixelcolors que representa una Imagen.

Answer 1

Una cosa importante a tener en cuenta cuando se avanza FFT seguida de FFT inversa es que normalmente se aplica un factor de escala de N al resultado final, es decir, los valores de píxeles de la imagen resultante deberán dividirse por N en para que coincida con los valores de píxeles originales. (Siendo N el tamaño de la FFT.) De modo que su lazo de salida debe probablemente ser algo como esto:

//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); 
    } 
} 

También tenga en cuenta que es probable que desee hacer una FFT-real a complejo seguido de un complejo-a- IFFT real (algo más eficiente en términos de memoria y rendimiento). Por ahora, aunque parece que estás haciendo complejo a complejo en ambas direcciones, lo cual está bien, pero no estás llenando tus matrices de entrada correctamente. Si usted va a seguir con complejo a complejo entonces es probable que desee cambiar su bucle de entrada a algo como esto:

//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; 
    } 
} 

es decir, los valores de los píxeles en los sitios reales de los valores de entrada complejas y la las partes imaginarias necesitan ser puestas a cero.

Una cosa más a tener en cuenta: cuando finalmente funcione, encontrará que el rendimiento es terrible; lleva mucho tiempo crear un plan relativo al tiempo que lleva la FFT real. La idea es crear el plan solo una vez, pero usarlo para realizar muchas FFT. Así que querrá separar la creación del plan del código FFT real y ponerlo en una rutina de inicialización o constructor o lo que sea.

Answer 2

Pero si usa la función RealToComplex o ComplexToReal, preste atención al hecho de que la imagen se almacenará en una matriz de dimensiones [altura x (ancho/2 +1)] y si desea hacer algunos cálculos intermedios en el dominio de la frecuencia, se volverán un poco más difíciles ...