00001 #include "sadie.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 void
00040 FILTER (IMAGE * in,
00041 double (*flt) (double, double, double, double),
00042
00043 double xradius,
00044 double yradius,
00045 double option,
00046
00047
00048
00049
00050 IMAGE ** out
00051
00052
00053 )
00054 {
00055 register short i, j, k;
00056 char msg[SLEN];
00057 double fx, fy, ampl, pinc, psum;
00058 IMAGE *scratch = 0;
00059
00060 if (TIMES)
00061 TIMING (T_FILTER);
00062 if (NAMES)
00063 {
00064 MESSAGE ('I', "");
00065 MESSAGE ('I', "FILTER");
00066 MESSAGE ('I', "");
00067 sprintf (msg, " Input image: %s", in->text);
00068 MESSAGE ('I', msg);
00069 if (flt == FLTBOX)
00070 {
00071 MESSAGE ('I', " Box filter");
00072 }
00073 else if (flt == FLTCONE)
00074 {
00075 MESSAGE ('I', " Cone filter");
00076 }
00077 else if (flt == FLTCYL)
00078 {
00079 MESSAGE ('I', " Cylinder filter");
00080 }
00081 else if (flt == FLTEXP)
00082 {
00083 MESSAGE ('I', " Exponential filter");
00084 }
00085 else if (flt == FLTGAUSS)
00086 {
00087 MESSAGE ('I', " Gaussian filter");
00088 }
00089 sprintf (msg, " Filter radius: cycles/line: %8.4f", yradius);
00090 MESSAGE ('I', msg);
00091 sprintf (msg, " cycles/pixel: %8.4f", xradius);
00092 MESSAGE ('I', msg);
00093 if (option == 0.0)
00094 {
00095 MESSAGE ('I', " Low-pass filter");
00096 }
00097 else if (option == 1.0)
00098 {
00099 MESSAGE ('I', " High-pass filter");
00100 }
00101 else
00102 {
00103 sprintf (msg, " Parametric high-boost filter, c = %5.2f", option);
00104 MESSAGE ('I', msg);
00105 }
00106 MESSAGE ('I', " ...............");
00107 }
00108
00109
00110 if (!CHECKIMG (in))
00111 {
00112 MESSAGE ('E', " Can't identify image.");
00113 goto the_end;
00114 }
00115 else if (in->nlin != 1L << rnd (log10 ((double) in->nlin) / log10 (2.0)))
00116 {
00117 MESSAGE ('E', " Number of lines must be a power of two.");
00118 goto the_end;
00119 }
00120 else if (in->npix != 1L << rnd (log10 ((double) in->npix) / log10 (2.0)))
00121 {
00122 MESSAGE ('E', " Number of pixels/line must be a power of two.");
00123 goto the_end;
00124 }
00125 else if (!flt)
00126 {
00127 MESSAGE ('E', " Can't identify filter function.");
00128 goto the_end;
00129 }
00130
00131
00132 FFT2D (in, FORWARD, &scratch);
00133 if (ABORT || !scratch)
00134 goto the_end;
00135 sprintf (scratch->text, "Scratch");
00136
00137
00138 if (LINES && !PROGRESS (psum = 0.0))
00139 goto the_end;
00140 pinc = 1.0 / (double) scratch->nbnd / (double) (scratch->nlin / 2 - 1);
00141
00142
00143 for (i = 0; i < scratch->nbnd; i++)
00144 {
00145 for (j = 1; j < scratch->nlin / 2; j++)
00146 {
00147 fy = 0.5 - (double) j / (double) scratch->nlin;
00148 for (k = 2; k < scratch->npix / 2 - 1; k += 2)
00149 {
00150 fx = 0.5 - (double) k / (double) (scratch->npix);
00151 ampl = (*flt) (fx, fy, xradius, yradius);
00152 if (option != 0.0)
00153 ampl = option - ampl;
00154
00155
00156 scratch->data[i][j][k] *= ampl;
00157 scratch->data[i][j][scratch->npix - k] *= ampl;
00158
00159
00160 scratch->data[i][scratch->nlin - j][k] =
00161 scratch->data[i][j][scratch->npix - k];
00162 scratch->data[i][scratch->nlin - j][scratch->npix - k] =
00163 scratch->data[i][j][k];
00164
00165
00166 scratch->data[i][j][k + 1] *= ampl;
00167 scratch->data[i][j][scratch->npix - k + 1] *= ampl;
00168
00169
00170 scratch->data[i][scratch->nlin - j][k + 1] =
00171 -scratch->data[i][j][scratch->npix - k + 1];
00172 scratch->data[i][scratch->nlin - j][scratch->npix - k + 1] =
00173 -scratch->data[i][j][k + 1];
00174 }
00175 if (LINES && !PROGRESS (psum += pinc))
00176 goto the_end;
00177 }
00178 }
00179
00180
00181 for (i = 0; i < scratch->nbnd; i++)
00182 {
00183 for (j = 1; j < scratch->nlin; j++)
00184 {
00185 fx = 0.0;
00186 fy = 0.5 - (double) j / (double) scratch->nlin;
00187 ampl = (*flt) (fx, fy, xradius, yradius);
00188 if (option != 0.0)
00189 ampl = option - ampl;
00190 scratch->data[i][j][scratch->npix / 2] *= ampl;
00191 scratch->data[i][j][scratch->npix / 2 + 1] *= ampl;
00192 }
00193 }
00194 for (i = 0; i < scratch->nbnd; i++)
00195 {
00196 for (k = 2; k < scratch->npix; k += 2)
00197 {
00198 fx = 0.5 - (double) k / (double) scratch->npix;
00199 fy = 0.0;
00200 ampl = (*flt) (fx, fy, xradius, yradius);
00201 if (option != 0.0)
00202 ampl = option - ampl;
00203 scratch->data[i][scratch->nlin / 2][k] *= ampl;
00204 scratch->data[i][scratch->nlin / 2][k + 1] *= ampl;
00205 }
00206 }
00207
00208
00209 for (i = 0; i < scratch->nbnd; i++)
00210 {
00211 for (j = 0; j < scratch->nlin; j++)
00212 {
00213 fx = 0.5;
00214 fy = 0.5 - (double) j / (double) scratch->nlin;
00215 ampl = (*flt) (fx, fy, xradius, yradius);
00216 if (option != 0.0)
00217 ampl = option - ampl;
00218 scratch->data[i][j][0] *= ampl;
00219 scratch->data[i][j][1] *= ampl;
00220 }
00221 }
00222 for (i = 0; i < scratch->nbnd; i++)
00223 {
00224 for (k = 2; k < scratch->npix; k += 2)
00225 {
00226 fx = 0.5 - (double) k / (double) scratch->npix;
00227 fy = 0.5;
00228 ampl = (*flt) (fx, fy, xradius, yradius);
00229 if (option != 0.0)
00230 ampl = option - ampl;
00231 scratch->data[i][0][k] *= ampl;
00232 scratch->data[i][0][k + 1] *= ampl;
00233 }
00234 }
00235
00236
00237 FFT2D (scratch, INVERSE, out);
00238
00239 the_end:
00240 if (scratch)
00241 RELMEM (scratch);
00242 if (TIMES)
00243 TIMING (T_EXIT);
00244 }
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255 double
00256 FLTBOX (double fx,
00257 double fy,
00258 double xradius,
00259 double yradius
00260
00261 )
00262 {
00263 double ampl;
00264
00265 if (TIMES)
00266 TIMING (T_FLTBOX);
00267
00268 ampl = fabs (fx / xradius) <= 1.0 && fabs (fy / yradius) <= 1.0 ? 1.0 : 0.0;
00269
00270 the_end:
00271 if (TIMES)
00272 TIMING (T_EXIT);
00273 return (ampl);
00274 }
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285 double
00286 FLTCONE (double fx,
00287 double fy,
00288 double xradius,
00289 double yradius
00290
00291 )
00292 {
00293 double ampl;
00294
00295 if (TIMES)
00296 TIMING (T_FLTCONE);
00297
00298 ampl =
00299 max (0.0,
00300 1.0 - sqrt ((fx * fx) / (xradius * xradius) +
00301 (fy * fy) / (yradius * yradius)));
00302
00303 the_end:
00304 if (TIMES)
00305 TIMING (T_EXIT);
00306 return (ampl);
00307 }
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318 double
00319 FLTCYL (double fx,
00320 double fy,
00321 double xradius,
00322 double yradius
00323
00324 )
00325 {
00326 double ampl;
00327
00328 if (TIMES)
00329 TIMING (T_FLTCYL);
00330
00331 ampl =
00332 sqrt ((fx * fx) / (xradius * xradius) +
00333 (fy * fy) / (yradius * yradius)) <= 1.0 ? 1.0 : 0.0;
00334
00335 the_end:
00336 if (TIMES)
00337 TIMING (T_EXIT);
00338 return (ampl);
00339 }
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351 double
00352 FLTEXP (double fx,
00353 double fy,
00354 double xradius,
00355 double yradius
00356
00357 )
00358 {
00359 double ampl;
00360
00361 if (TIMES)
00362 TIMING (T_FLTEXP);
00363
00364 ampl =
00365 exp (-sqrt
00366 ((fx * fx) / (xradius * xradius) + (fy * fy) / (yradius * yradius)));
00367
00368 the_end:
00369 if (TIMES)
00370 TIMING (T_EXIT);
00371 return (ampl);
00372 }
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383 double
00384 FLTGAUSS (double fx,
00385 double fy,
00386 double xradius,
00387 double yradius
00388
00389 )
00390 {
00391 double ampl;
00392
00393 if (TIMES)
00394 TIMING (T_FLTGAUSS);
00395
00396 ampl =
00397 exp (-
00398 ((fx * fx) / (xradius * xradius) + (fy * fy) / (yradius * yradius)));
00399
00400 the_end:
00401 if (TIMES)
00402 TIMING (T_EXIT);
00403 return (ampl);
00404 }