00001 #include "sadie.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 void
00014 FFT1D (short option,
00015
00016
00017 short size,
00018 double *sine,
00019 double *cosn,
00020 double *array
00021
00022 )
00023 {
00024 register short n, m, l, k, j, i;
00025 register double temp, e_real, e_imag, o_real, o_imag;
00026
00027 if (TIMES)
00028 TIMING (T_FFT1D);
00029
00030
00031 for (i = j = 0; j < size; j += 1, i += k)
00032 {
00033 if (i > j)
00034 {
00035 temp = array[2 * i];
00036 array[2 * i] = array[2 * j];
00037 array[2 * j] = temp;
00038 temp = array[2 * i + 1];
00039 array[2 * i + 1] = array[2 * j + 1];
00040 array[2 * j + 1] = temp;
00041 }
00042 for (k = size / 2; 1 <= k && k <= i; i -= k, k /= 2);
00043 }
00044
00045
00046 for (i = 2; i <= size; i *= 2)
00047 {
00048 for (j = 0; j < size / i; j++)
00049 {
00050 for (k = 0; k < i / 2; k++)
00051 {
00052 l = 2 * (j * i + k);
00053 m = 2 * (j * i + i / 2 + k);
00054 n = (short) ((long) k * (long) size / (long) i);
00055 e_real = array[l];
00056 e_imag = array[l + 1];
00057 o_real = array[m] * cosn[n] - array[m + 1] * sine[n];
00058 o_imag = array[m] * sine[n] + array[m + 1] * cosn[n];
00059 array[l] = (e_real + o_real) / 2.0;
00060 array[l + 1] = (e_imag + o_imag) / 2.0;
00061 array[m] = (e_real - o_real) / 2.0;
00062 array[m + 1] = (e_imag - o_imag) / 2.0;
00063 }
00064 }
00065 }
00066
00067 if (option == INVERSE)
00068 for (i = 0; i < 2 * size; array[i++] *= (double) size);
00069
00070 the_end:
00071 if (TIMES)
00072 TIMING (T_EXIT);
00073 }
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091 void
00092 FFT2D (IMAGE * in,
00093 short option,
00094
00095
00096 IMAGE ** out
00097
00098 )
00099 {
00100 register short i, j, k;
00101 char msg[SLEN];
00102 double *sine = 0, *cosn = 0, *array = 0, pinc, psum;
00103 PIXEL temp;
00104 IMAGE *img;
00105
00106 if (TIMES)
00107 TIMING (T_FFT2D);
00108 if (NAMES)
00109 {
00110 MESSAGE ('I', "");
00111 MESSAGE ('I', "FFT2D");
00112 MESSAGE ('I', "");
00113 sprintf (msg, " Input image: %s", in->text);
00114 MESSAGE ('I', msg);
00115 if (option == FORWARD)
00116 {
00117 MESSAGE ('I', " Forward FFT");
00118 }
00119 else if (option == INVERSE)
00120 {
00121 MESSAGE ('I', " Inverse FFT");
00122 }
00123 MESSAGE ('I', " ...............");
00124 }
00125
00126
00127 if (!CHECKIMG (in))
00128 {
00129 MESSAGE ('E', " Can't identify image.");
00130 goto the_end;
00131 }
00132 else if (in->nlin != 1L << rnd (log10 ((double) in->nlin) / log10 (2.0)))
00133 {
00134 MESSAGE ('E', " Number of lines must be a power of two.");
00135 goto the_end;
00136 }
00137 else if (in->npix != 1L << rnd (log10 ((double) in->npix) / log10 (2.0)))
00138 {
00139 MESSAGE ('E', " Number of pixels/line must be a power of two.");
00140 goto the_end;
00141 }
00142 else if (option != FORWARD && option != INVERSE)
00143 {
00144 MESSAGE ('E',
00145 " Transformation option must be either FORWARD or INVERSE.");
00146 goto the_end;
00147 }
00148
00149
00150 if (!CHECKIMG (*out))
00151 GETMEM (in->nbnd, in->nlin,
00152 (short) ((option == FORWARD ? 2.0 : 0.5) * in->npix), out);
00153 if (!*out)
00154 goto the_end;
00155 img = option == FORWARD ? *out : in;
00156
00157
00158 i = max (img->nlin, img->npix / 2);
00159 if (!(sine = (double *) malloc (i * sizeof (double))))
00160 {
00161 MESSAGE ('E', " Memory request failed.");
00162 goto the_end;
00163 }
00164 if (!(cosn = (double *) malloc (i * sizeof (double))))
00165 {
00166 MESSAGE ('E', " Memory request failed.");
00167 goto the_end;
00168 }
00169 for (j = 0; j < i; j++)
00170 {
00171 sine[j] = (double) option *sin ((2.0 * PI * (double) j) / (double) i);
00172 cosn[j] = cos ((2.0 * PI * (double) j) / (double) i);
00173 }
00174
00175
00176 if (!(array = (double *) malloc (2 * i * sizeof (double))))
00177 {
00178 MESSAGE ('E', " Memory request failed.");
00179 goto the_end;
00180 }
00181
00182
00183 if (LINES && !PROGRESS (psum = 0.0))
00184 goto the_end;
00185 pinc = 1.0 / (double) img->nbnd / (double) (img->nlin + img->npix / 2);
00186
00187
00188 for (i = 0; i < img->nbnd; i++)
00189 {
00190
00191
00192 if (option == FORWARD)
00193 {
00194 for (j = 0; j < (*out)->nlin; j += 1)
00195 {
00196 for (k = 0; k < (*out)->npix; k += 2)
00197 {
00198 (*out)->data[i][j][k] = in->data[i][j][k / 2];
00199 (*out)->data[i][j][k + 1] = (PIXEL) 0;
00200 }
00201 }
00202 }
00203
00204
00205 for (j = 0; j < img->nlin / 2; j++)
00206 {
00207 for (k = 0; k < img->npix / 2; k++)
00208 {
00209 temp = img->data[i][j][k];
00210 img->data[i][j][k] =
00211 img->data[i][img->nlin / 2 + j][img->npix / 2 + k];
00212 img->data[i][img->nlin / 2 + j][img->npix / 2 + k] = temp;
00213 temp = img->data[i][img->nlin / 2 + j][k];
00214 img->data[i][img->nlin / 2 + j][k] =
00215 img->data[i][j][img->npix / 2 + k];
00216 img->data[i][j][img->npix / 2 + k] = temp;
00217 }
00218 }
00219
00220
00221 for (j = 0; j < img->nlin; j++)
00222 {
00223 for (k = 0; k < img->npix; k++)
00224 {
00225 array[k] = (double) img->data[i][j][k];
00226 }
00227 FFT1D (option, img->npix / 2, sine, cosn, array);
00228 for (k = 0; k < img->npix; k++)
00229 {
00230 img->data[i][j][k] = (PIXEL) array[k];
00231 }
00232 if (LINES && !PROGRESS (psum += pinc))
00233 goto the_end;
00234 }
00235
00236
00237 for (k = 0; k < img->npix; k += 2)
00238 {
00239 for (j = 0; j < img->nlin; j += 1)
00240 {
00241 array[2 * j] = (double) img->data[i][j][k];
00242 array[2 * j + 1] = (double) img->data[i][j][k + 1];
00243 }
00244 FFT1D (option, img->nlin, sine, cosn, array);
00245 for (j = 0; j < img->nlin; j += 1)
00246 {
00247 img->data[i][j][k] = (PIXEL) array[2 * j];
00248 img->data[i][j][k + 1] = (PIXEL) array[2 * j + 1];
00249 }
00250 if (LINES && !PROGRESS (psum += pinc))
00251 goto the_end;
00252 }
00253
00254
00255 for (j = 0; j < img->nlin / 2; j++)
00256 {
00257 for (k = 0; k < img->npix / 2; k++)
00258 {
00259 temp = img->data[i][j][k];
00260 img->data[i][j][k] =
00261 img->data[i][img->nlin / 2 + j][img->npix / 2 + k];
00262 img->data[i][img->nlin / 2 + j][img->npix / 2 + k] = temp;
00263 temp = img->data[i][img->nlin / 2 + j][k];
00264 img->data[i][img->nlin / 2 + j][k] =
00265 img->data[i][j][img->npix / 2 + k];
00266 img->data[i][j][img->npix / 2 + k] = temp;
00267 }
00268 }
00269
00270
00271 if (option == INVERSE)
00272 {
00273 for (j = 0; j < in->nlin; j += 1)
00274 {
00275 for (k = 0; k < in->npix; k += 2)
00276 {
00277 (*out)->data[i][j][k / 2] = in->data[i][j][k];
00278 }
00279 }
00280 }
00281 }
00282
00283 the_end:
00284 if (sine)
00285 free (sine);
00286 if (cosn)
00287 free (cosn);
00288 if (array)
00289 free (array);
00290 if (TIMES)
00291 TIMING (T_EXIT);
00292 }