00001
00002
00003 #include <math.h>
00004 #include <sys/wait.h>
00005 #include <signal.h>
00006 #include "misc.h"
00007 #include "config.h"
00008 #include "FourierTransform.h"
00009
00010 #ifdef HAVE_LIBFFTW
00011 #include "fftwlib.h"
00012 #endif
00013
00014 void FourierTransform::_fft(float *data, int nn, int isign)
00015 {
00016 int n, mmax, m, j, istep, i;
00017 float wtemp, wr, wpr, wpi, wi, theta;
00018 float tempr, tempi;
00019
00020 n = nn<<1;
00021 j = 1;
00022 for (i=1;i<n;i+=2) {
00023 if (j>i) {
00024 SWAP(data[j],data[i]);
00025 SWAP(data[j+1],data[i+1]);
00026 }
00027 m = n>>1;
00028 while((m>=2) && (j>m)) {
00029 j -= m;
00030 m >>= 1;
00031 }
00032 j+=m;
00033 }
00034
00035 mmax = 2;
00036 while( n > mmax) {
00037 istep = mmax<<1;
00038 theta = isign*(6.28318530/mmax);
00039 wtemp = sin(0.5*theta);
00040 wpr = -2.0*wtemp*wtemp;
00041 wpi = sin(theta);
00042 wr = 1.0;
00043 wi = 0.0;
00044 for (m=1;m<mmax;m+=2) {
00045 for (i=m;i<=n;i+=istep) {
00046 j=i+mmax;
00047 tempr = wr*data[j]-wi*data[j+1];
00048 tempi = wr*data[j+1]+wi*data[j];
00049 data[j] = data[i]-tempr;
00050 data[j+1] = data[i+1]-tempi;
00051 data[i] += tempr;
00052 data[i+1] += tempi;
00053 }
00054 wr = (wtemp=wr)*wpr-wi*wpi+wr;
00055 wi = wi*wpr+wtemp*wpi+wi;
00056 }
00057 mmax=istep;
00058 }
00059 }
00060
00061 FourierTransform::FourierTransform(FloatImage &input)
00062 {
00063 #ifdef HAVE_LIBFFTW
00064 input.makeEvenSized();
00065 real_part = input;
00066 imag_part.makeImage(input.getWidth(),input.getHeight());
00067 imag_part.clearImage(0);
00068 #else
00069 input.makeEvenSized();
00070 real_part = input;
00071 int w = real_part.getWidth();
00072 int h = real_part.getHeight();
00073 int p2_w = 1; int p2_h = 1;
00074
00075 for(int i=0;i<16;i++) {
00076 if (w <= p2_w) break;
00077 p2_w = p2_w<<1;
00078 }
00079 for(int i=0;i<16;i++) {
00080 if (h <= p2_h) break;
00081 p2_h = p2_h<<1;
00082 }
00083 real_part.pad((p2_w - w)/2, (p2_h - h)/2, ZEROS);
00084 imag_part.makeImage(real_part.getWidth(),real_part.getHeight());
00085 imag_part.clearImage(0);
00086 #endif
00087 }
00088
00089 FourierTransform::FourierTransform(FloatImage &input_real, FloatImage &input_imag)
00090 {
00091 if ((input_real.getWidth() != input_imag.getWidth()) ||
00092 (input_real.getHeight() != input_imag.getHeight())) {
00093 ERROR("Real and imaginary input images have different sizes");
00094 }
00095 else {
00096 #ifdef HAVE_LIBFFTW
00097 input_real.makeEvenSized();
00098 input_imag.makeEvenSized();
00099 real_part = input_real;
00100 imag_part = input_imag;
00101 #else
00102 int w = input_real.getWidth();
00103 int h = input_imag.getHeight();
00104 int p2_w = 1; int p2_h = 1;
00105
00106 int flag_w = 0, flag_h =0;
00107 for(int i=0;i<16;i++) {
00108 if (w == p2_w) {
00109 flag_w = 1;
00110 break;
00111 }
00112 p2_w = p2_w<<1;
00113 }
00114 if (!flag_w)
00115 ERROR("Image width not a power of 2; results may be incorrect");
00116 for(int i=0;i<16;i++) {
00117 if (h == p2_h) {
00118 flag_h = 1;
00119 break;
00120 }
00121 p2_h = p2_h<<1;
00122 }
00123 if (!flag_h)
00124 ERROR("Image height not a power of 2; results may be incorrect");
00125 real_part = input_real;
00126 imag_part = input_imag;
00127 #endif
00128 }
00129 }
00130
00131 void FourierTransform::fft(void)
00132 {
00133
00134 #ifdef HAVE_LIBFFTW
00135 FloatImage temp_real = real_part;
00136 FloatImage temp_imag = imag_part;
00137 _fft2d(temp_real.getPixels(), temp_imag.getPixels(),
00138 real_part.getWidth() , real_part.getHeight(),
00139 real_part.getPixels(), imag_part.getPixels());
00140 #else
00141 int w = real_part.getWidth();
00142 int h = real_part.getHeight();
00143
00144 float * temp = (float *) malloc(2*w* sizeof(float));
00145 for (int i=0;i<h;i++) {
00146 for(int j=0;j<w;j++) {
00147 temp[2*j] = real_part[i][j];
00148 temp[2*j+1] = imag_part[i][j];
00149 }
00150 _fft(temp-1, w, 1);
00151 for(int j=0;j<w;j++) {
00152 real_part[i][j] = temp[2*j];
00153 imag_part[i][j] = temp[2*j+1];
00154 }
00155 }
00156 free(temp);
00157 temp = (float *) malloc(2*h * sizeof(float));
00158 for(int j=0;j<w;j++) {
00159 for (int i=0;i<h;i++) {
00160 temp[2*i] = real_part[i][j];
00161 temp[2*i+1] = imag_part[i][j];
00162 }
00163 _fft(temp-1, h, 1);
00164 for(int i=0;i<h;i++) {
00165 real_part[i][j] = temp[2*i];
00166 imag_part[i][j] = temp[2*i+1];
00167 }
00168 }
00169 free(temp);
00170 for(int j=0;j<w;j++)
00171 for(int i=0;i<h;i++) {
00172 real_part[i][j] /= w*h;
00173 imag_part[i][j] /= w*h;
00174 }
00175 HALF_SHIFT(real_part.getPixels(),w,h);
00176 HALF_SHIFT(imag_part.getPixels(),w,h);
00177 #endif
00178 }
00179
00180 void FourierTransform::ifft(void)
00181 {
00182 #ifdef HAVE_LIBFFTW
00183 FloatImage temp_real = real_part;
00184 FloatImage temp_imag = imag_part;
00185 _ifft2d(temp_real.getPixels(), temp_imag.getPixels(),
00186 real_part.getWidth() , real_part.getHeight(),
00187 real_part.getPixels(), imag_part.getPixels());
00188 #else
00189 int w = real_part.getWidth();
00190 int h = real_part.getHeight();
00191
00192 HALF_SHIFT(real_part.getPixels(),w,h);
00193 HALF_SHIFT(imag_part.getPixels(),w,h);
00194 float *temp = (float *) malloc(2*w * sizeof(float));
00195 for (int i=0;i<h;i++) {
00196 for(int j=0;j<w;j++) {
00197 temp[2*j] = real_part[i][j];
00198 temp[2*j+1] = imag_part[i][j];
00199 }
00200 _fft(temp-1, w, -1);
00201 for(int j=0;j<w;j++) {
00202 real_part[i][j] = temp[2*j];
00203 imag_part[i][j] = temp[2*j+1];
00204 }
00205 }
00206 free(temp);
00207 temp = (float *) malloc(2*h * sizeof(float));
00208 for(int j=0;j<w;j++) {
00209 for (int i=0;i<h;i++) {
00210 temp[2*i] = real_part[i][j];
00211 temp[2*i+1] = imag_part[i][j];
00212 }
00213 _fft(temp-1, h, -1);
00214 for(int i=0;i<h;i++) {
00215 real_part[i][j] = temp[2*i];
00216 imag_part[i][j] = temp[2*i+1];
00217 }
00218 }
00219 free(temp);
00220 #endif
00221 }
00222
00223 void FourierTransform::showMagnitude(int log_values)
00224 {
00225 FloatImage magnitude(real_part.getWidth(), real_part.getHeight());
00226
00227 for (int j=0;j<real_part.getWidth();j++)
00228 for (int i=0;i<real_part.getHeight();i++)
00229 magnitude[i][j] = sqrt(real_part[i][j]*real_part[i][j] +
00230 imag_part[i][j]*imag_part[i][j]);
00231 if (log_values) {
00232 for (int j=0;j<real_part.getWidth();j++)
00233 for (int i=0;i<real_part.getHeight();i++)
00234 magnitude[i][j] = log(magnitude[i][j]);
00235 }
00236 if(pid_magnitude) {
00237 kill(pid_magnitude, SIGUSR1);
00238 waitpid(pid_magnitude, NULL, 0);
00239 pid_magnitude = magnitude.showImage();
00240 }
00241 else
00242 pid_magnitude = magnitude.showImage();
00243 }
00244
00245 void FourierTransform::showPhase(int log_values)
00246 {
00247 FloatImage phase(imag_part.getWidth(), imag_part.getHeight());
00248 for (int j=0;j<imag_part.getWidth();j++)
00249 for (int i=0;i<imag_part.getHeight();i++)
00250 phase[i][j] = atan(imag_part[i][j]/real_part[i][j]);
00251 if (log_values) {
00252 for (int j=0;j<imag_part.getWidth();j++)
00253 for (int i=0;i<imag_part.getHeight();i++)
00254 phase[i][j] = log(phase[i][j]);
00255 }
00256 if(pid_phase) {
00257 kill(pid_phase, SIGUSR1);
00258 waitpid(pid_phase, NULL, 0);
00259 pid_phase = phase.showImage();
00260 }
00261 else
00262 pid_phase = phase.showImage();
00263 }
00264
00265 void FourierTransform::closeWindows()
00266 {
00267 if(pid_magnitude) {
00268 kill(pid_magnitude, SIGUSR1);
00269 waitpid(pid_magnitude, NULL, 0);
00270 pid_magnitude = 0;
00271 }
00272 if(pid_phase) {
00273 kill(pid_phase, SIGUSR1);
00274 waitpid(pid_phase, NULL, 0);
00275 pid_phase = 0;
00276 }
00277 }
00278
00279 void fft(FloatImage &input, FloatImage &output_real, FloatImage &output_imag)
00280 {
00281 FourierTransform a(input);
00282 a.fft();
00283 output_real = a.getReal();
00284 output_imag = a.getImag();
00285 }
00286
00287 void ifft(FloatImage &input_real, FloatImage &input_imag, FloatImage &output)
00288 {
00289 FourierTransform a(input_real, input_imag);
00290 a.ifft();
00291 output = a.getReal();
00292 }
00293