00001
00002
00003 #include <stdlib.h>
00004 #include <stdio.h>
00005 #include <string.h>
00006 #include <math.h>
00007 #include <config.h>
00008
00009 #ifdef HAVE_LIBFFTW
00010
00011 #include <fftw.h>
00012 #include <rfftw.h>
00013 #include "fftwlib.h"
00014
00015
00016
00017
00018
00019 void _convolve(float *input_image_1,
00020 float *input_image_2,
00021 int size_x, int size_y,
00022 float *convolved_image)
00023 {
00024
00025 fftw_real a[size_x][2*(size_y/2+1)], b[size_x][2*(size_y/2+1)],c[size_x][size_y];
00026 fftw_complex *A, *B, C[size_x][size_y/2+1];
00027 fftwnd_plan p, pinv;
00028 fftw_real scale = 1.0 / (size_x * size_y);
00029 int i, j, ij;
00030
00031 p = rfftw2d_create_plan(size_x, size_y,
00032 FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE | FFTW_IN_PLACE);
00033 pinv = rfftw2d_create_plan(size_x, size_y,
00034 FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
00035
00036 half_shift(input_image_1,size_x,size_y);
00037
00038
00039 A = (fftw_complex *) &a[0][0];
00040 B = (fftw_complex *) &b[0][0];
00041
00042 for(i=0; i<size_x; ++i)
00043 for(j=0; j<size_y; ++j)
00044 {
00045 a[i][j] = input_image_1[j*size_x+i];
00046 b[i][j] = input_image_2[j*size_x+i];
00047 }
00048
00049 rfftwnd_one_real_to_complex(p,&a[0][0],NULL);
00050 rfftwnd_one_real_to_complex(p,&b[0][0],NULL);
00051
00052 for(i=0;i<size_x;++i)
00053 for(j=0;j<(size_y/2+1);++j)
00054 {
00055 ij = i*(size_y/2+1)+j;
00056 C[i][j].re = (A[ij].re * B[ij].re -
00057 A[ij].im * B[ij].im)*scale;
00058
00059 C[i][j].im = (A[ij].re * B[ij].im +
00060 A[ij].im * B[ij].re)*scale;
00061 }
00062
00063 rfftwnd_one_complex_to_real(pinv,&C[0][0],&c[0][0]);
00064
00065 for(i=0; i<size_x; ++i)
00066 for(j=0; j<size_y; ++j)
00067 convolved_image[j*size_x+i] = (float)c[i][j];
00068
00069 half_shift(input_image_1,size_x,size_y);
00070
00071 rfftwnd_destroy_plan(p);
00072 rfftwnd_destroy_plan(pinv);
00073 }
00074
00075 void _fft2d(float *real_input_image, float *imag_input_image,
00076 int size_x, int size_y,
00077 float *fourier_real_part, float *fourier_imag_part)
00078 {
00079
00080 fftw_complex a[size_x*size_y];
00081 fftwnd_plan p;
00082 fftw_real scale;
00083 int i, j, ij;
00084
00085 scale = 1.0 / (size_x * size_y);
00086
00087 p = fftw2d_create_plan(size_x, size_y, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_IN_PLACE);
00088
00089 half_shift(real_input_image,size_x,size_y);
00090 if (imag_input_image != NULL)
00091 half_shift(imag_input_image,size_x,size_y);
00092
00093 for(j=0; j<size_y; ++j)
00094 for(i=0; i<size_x; ++i)
00095 {
00096 ij = i*size_y+j;
00097 a[ij].re = real_input_image[j*size_x+i];
00098 if (imag_input_image != NULL)
00099 a[ij].im = imag_input_image[j*size_x+i];
00100 else
00101 a[ij].im = 0.0;
00102 }
00103
00104 fftwnd_one(p,&a[0],NULL);
00105
00106 for(j=0; j<size_y; ++j)
00107 for(i=0; i<size_x; ++i)
00108 {
00109 ij = i*size_y+j;
00110 fourier_real_part[j*size_x+i] = (float) a[ij].re;
00111 fourier_imag_part[j*size_x+i] = (float) a[ij].im;
00112 }
00113
00114 half_shift(fourier_real_part,size_x,size_y);
00115 half_shift(fourier_imag_part,size_x,size_y);
00116 fftwnd_destroy_plan(p);
00117 }
00118
00119 void _ifft2d(float *fourier_real_image, float *fourier_imag_image,
00120 int size_x, int size_y, float *real_image, float *imag_image)
00121 {
00122
00123
00124 fftw_complex a[size_x*size_y];
00125 fftwnd_plan p;
00126 fftw_real scale;
00127 int i, j, ij;
00128
00129 scale = 1.0 / (size_x * size_y);
00130
00131 p = fftw2d_create_plan(size_x, size_y, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_IN_PLACE);
00132
00133 half_shift(fourier_real_image,size_x,size_y);
00134 half_shift(fourier_imag_image,size_x,size_y);
00135
00136 for(j=0; j<size_y; ++j)
00137 for(i=0; i<size_x; ++i)
00138 {
00139 ij = i*size_y+j;
00140 a[ij].re = (double) fourier_real_image[j*size_x+i];
00141 a[ij].im = (double) fourier_imag_image[j*size_x+i];
00142 }
00143
00144 fftwnd_one(p,&a[0],NULL);
00145
00146 for(j=0; j<size_y; ++j)
00147 for(i=0; i<size_x; ++i)
00148 {
00149 ij = i*size_y+j;
00150 real_image[j*size_x+i] = (float) a[ij].re * scale;
00151 imag_image[j*size_x+i] = (float) a[ij].im * scale;
00152 }
00153
00154 half_shift(real_image,size_x,size_y);
00155 half_shift(imag_image,size_x,size_y);
00156 fftwnd_destroy_plan(p);
00157 }
00158
00159 #endif