00001 #include "sadie.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014 static void
00015 FHT1D (short option,
00016
00017
00018 short size,
00019 short p,
00020 short *a,
00021 short *pwr,
00022 double *sine,
00023 double *tang,
00024 double *array
00025
00026 )
00027 {
00028 register short o, n, m, l, k, j, i;
00029 register double t, u, v, w, x, y;
00030
00031 if (TIMES)
00032 TIMING (T_FHT1D);
00033
00034 a[0] = 0;
00035 a[1] = 1;
00036 for (i = 2; i <= p / 2 + p % 2; i++)
00037 {
00038 for (j = 0; j < pwr[i - 1]; j++)
00039 {
00040 a[j + pwr[i - 1]] = (a[j] *= 2) + 1;
00041 }
00042 }
00043 for (i = 1; i < pwr[p / 2]; i++)
00044 {
00045 k = i;
00046 l = m = pwr[p / 2] * a[i];
00047 t = array[k];
00048 array[k] = array[m];
00049 array[m] = t;
00050 for (j = 1; j <= a[i] - 1; j++)
00051 {
00052 k += pwr[p / 2];
00053 m = l + a[j];
00054 t = array[k];
00055 array[k] = array[m];
00056 array[m] = t;
00057 }
00058 }
00059 for (i = 0; i < size - 1; i += 2)
00060 {
00061 t = array[i] + array[i + 1];
00062 u = array[i] - array[i + 1];
00063 array[i] = t;
00064 array[i + 1] = u;
00065 }
00066 for (i = 0; i < size - 3; i += 4)
00067 {
00068 t = array[i] + array[i + 2];
00069 u = array[i + 1] + array[i + 3];
00070 v = array[i] - array[i + 2];
00071 w = array[i + 1] - array[i + 3];
00072 array[i] = t;
00073 array[i + 1] = u;
00074 array[i + 2] = v;
00075 array[i + 3] = w;
00076 }
00077 for (i = 2; i < p; i++)
00078 {
00079 for (j = 0; j < size; j += pwr[i + 1])
00080 {
00081 t = array[j] + array[j + pwr[i]];
00082 u = array[j] - array[j + pwr[i]];
00083 array[j] = t;
00084 array[j + pwr[i]] = u;
00085 l = j + 1;
00086 m = j + pwr[i] - 1;
00087 n = l + pwr[i];
00088 o = m + pwr[i];
00089 for (k = pwr[p - i - 1]; k <= size / 4; k += pwr[p - i - 1])
00090 {
00091 t = array[n] + tang[k] * array[o];
00092 x = array[o] - sine[k] * t;
00093 y = t + tang[k] * x;
00094 t = array[l] + y;
00095 u = array[m] - x;
00096 v = array[l] - y;
00097 w = array[m] + x;
00098 array[l++] = t;
00099 array[m--] = u;
00100 array[n++] = v;
00101 array[o--] = w;
00102 }
00103 }
00104 }
00105
00106 if (option == FORWARD)
00107 for (i = 0; i < size; array[i++] /= (double) size);
00108
00109 the_end:
00110 if (TIMES)
00111 TIMING (T_EXIT);
00112 }
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130 void
00131 FHT2D (IMAGE * in,
00132 short option,
00133
00134
00135 IMAGE ** out
00136
00137 )
00138 {
00139 register short i, j, k, p, *a = 0, *pwr = 0;
00140 char msg[SLEN];
00141 double *sine = 0, *tang = 0, *array = 0, pinc, psum;
00142 PIXEL t, u, v, w;
00143 IMAGE *img;
00144
00145 if (TIMES)
00146 TIMING (T_FHT2D);
00147 if (NAMES)
00148 {
00149 MESSAGE ('I', "");
00150 MESSAGE ('I', "FHT2D");
00151 MESSAGE ('I', "");
00152 sprintf (msg, " Input image: %s", in->text);
00153 MESSAGE ('I', msg);
00154 if (option == FORWARD)
00155 {
00156 MESSAGE ('I', " Forward FHT");
00157 }
00158 else if (option == INVERSE)
00159 {
00160 MESSAGE ('I', " Inverse FHT");
00161 }
00162 MESSAGE ('I', " ...............");
00163 }
00164
00165
00166 if (!CHECKIMG (in))
00167 {
00168 MESSAGE ('E', " Can't identify image.");
00169 goto the_end;
00170 }
00171 else if (in->nlin != 1L << rnd (log10 ((double) in->nlin) / log10 (2.0)))
00172 {
00173 MESSAGE ('E', " Number of lines must be a power of two.");
00174 goto the_end;
00175 }
00176 else if (in->npix != 1L << rnd (log10 ((double) in->npix) / log10 (2.0)))
00177 {
00178 MESSAGE ('E', " Number of pixels/line must be a power of two.");
00179 goto the_end;
00180 }
00181 else if (option != FORWARD && option != INVERSE)
00182 {
00183 MESSAGE ('E',
00184 " Transformation option must be either FORWARD or INVERSE.");
00185 goto the_end;
00186 }
00187 else if (option == FORWARD && in->nlin != in->npix || option == INVERSE
00188 && in->nlin != in->npix / 2)
00189 {
00190 MESSAGE ('E',
00191 " Number of lines must be equal to number of pixels/line.");
00192 goto the_end;
00193 }
00194
00195
00196 if (!CHECKIMG (*out))
00197 GETMEM (in->nbnd, in->nlin,
00198 (short) ((option == FORWARD ? 2.0 : 0.5) * in->npix), out);
00199 if (!*out)
00200 goto the_end;
00201 img = option == FORWARD ? in : *out;
00202
00203
00204 if (img->nlin == 1)
00205 {
00206 (*out)->data[0][0][0] = in->data[0][0][0];
00207 if (option == FORWARD)
00208 {
00209 (*out)->data[0][0][1] = 0.0;
00210 }
00211 goto the_end;
00212 }
00213
00214
00215 i = max (img->nlin, img->npix);
00216 p = rnd (log10 ((double) i) / log10 (2.0));
00217 if (!(a = (short *) malloc ((1 << p / 2 + p % 2) * sizeof (short))))
00218 {
00219 MESSAGE ('E', " Memory request failed.");
00220 goto the_end;
00221 }
00222
00223
00224 if (!(pwr = (short *) malloc ((p + 1) * sizeof (short))))
00225 {
00226 MESSAGE ('E', " Memory request failed.");
00227 goto the_end;
00228 }
00229 for (j = 0; j <= p; pwr[j] = 1 << j, j++);
00230
00231
00232 if (!(sine = (double *) malloc ((i / 4 + 1) * sizeof (double))))
00233 {
00234 MESSAGE ('E', " Memory request failed.");
00235 goto the_end;
00236 }
00237 if (!(tang = (double *) malloc ((i / 4 + 1) * sizeof (double))))
00238 {
00239 MESSAGE ('E', " Memory request failed.");
00240 goto the_end;
00241 }
00242 for (j = 0; j <= i / 4; j++)
00243 {
00244 sine[j] = sin ((2.0 * PI * (double) j) / (double) i);
00245 tang[j] = tan ((PI * (double) j) / (double) i);
00246 }
00247
00248
00249 if (!(array = (double *) malloc (i * sizeof (double))))
00250 {
00251 MESSAGE ('E', " Memory request failed.");
00252 goto the_end;
00253 }
00254
00255
00256 if (LINES && !PROGRESS (psum = 0.0))
00257 goto the_end;
00258 pinc = 1.0 / (double) img->nbnd / (double) (img->nlin + img->npix);
00259
00260
00261 for (i = 0; i < (*out)->nbnd; i++)
00262 {
00263
00264 if (option == FORWARD)
00265 {
00266
00267 for (j = 0; j < (*out)->nlin / 2; j++)
00268 {
00269 for (k = 0; k < (*out)->npix / 4; k++)
00270 {
00271 (*out)->data[i][j][k] =
00272 in->data[i][in->nlin / 2 + j][in->npix / 2 + k];
00273 (*out)->data[i][(*out)->nlin / 2 + j][k] =
00274 in->data[i][j][in->npix / 2 + k];
00275 (*out)->data[i][j][(*out)->npix / 4 + k] =
00276 in->data[i][in->nlin / 2 + j][k];
00277 (*out)->data[i][(*out)->nlin / 2 + j][(*out)->npix / 4 +
00278 k] =
00279 in->data[i][j][k];
00280 }
00281 }
00282 }
00283 else
00284 {
00285
00286 for (j = 0; j < (*out)->nlin / 2; j++)
00287 {
00288 for (k = 0; k < (*out)->npix / 2; k++)
00289 {
00290 (*out)->data[i][j][k] =
00291 in->data[i][in->nlin / 2 + j][in->npix / 2 + 2 * k] -
00292 in->data[i][in->nlin / 2 + j][in->npix / 2 + 2 * k + 1];
00293 (*out)->data[i][(*out)->nlin / 2 + j][k] =
00294 in->data[i][j][in->npix / 2 + 2 * k] -
00295 in->data[i][j][in->npix / 2 + 2 * k + 1];
00296 (*out)->data[i][j][(*out)->npix / 2 + k] =
00297 in->data[i][in->nlin / 2 + j][2 * k] -
00298 in->data[i][in->nlin / 2 + j][2 * k + 1];
00299 (*out)->data[i][(*out)->nlin / 2 + j][(*out)->npix / 2 +
00300 k] =
00301 in->data[i][j][2 * k] - in->data[i][j][2 * k + 1];
00302 }
00303 }
00304 }
00305
00306
00307 for (j = 0; j < img->nlin; j++)
00308 {
00309 for (k = 0; k < img->npix; k++)
00310 {
00311 array[k] = (double) (*out)->data[i][j][k];
00312 }
00313 FHT1D (option, img->npix, p, a, pwr, sine, tang, array);
00314 for (k = 0; k < img->npix; k++)
00315 {
00316 (*out)->data[i][j][k] = (PIXEL) array[k];
00317 }
00318 if (LINES && !PROGRESS (psum += pinc))
00319 goto the_end;
00320 }
00321
00322
00323 for (k = img->npix - 1; k >= 0; k--)
00324 {
00325 for (j = 0; j < img->nlin; j++)
00326 {
00327 array[j] = (double) (*out)->data[i][j][k];
00328 }
00329 FHT1D (option, img->nlin, p, a, pwr, sine, tang, array);
00330 if (option == FORWARD)
00331 {
00332 for (j = 0; j < img->nlin; j++)
00333 {
00334 (*out)->data[i][j][2 * k] = (PIXEL) array[j];
00335 }
00336 }
00337 else
00338 {
00339 for (j = 0; j < img->nlin; j++)
00340 {
00341 (*out)->data[i][j][k] = (PIXEL) array[j];
00342 }
00343 }
00344 if (LINES && !PROGRESS (psum += pinc))
00345 goto the_end;
00346 }
00347
00348 if (option == FORWARD)
00349 {
00350
00351 (*out)->data[i][0][1] = (*out)->data[i][(*out)->nlin / 2][1] =
00352 (*out)->data[i][0][(*out)->npix / 2 + 1] =
00353 (*out)->data[i][(*out)->nlin / 2][(*out)->npix / 2 + 1] = 0.0;
00354 for (k = 2; k < (*out)->npix / 2; k += 2)
00355 {
00356 (*out)->data[i][0][(*out)->npix - k + 1] =
00357 -((*out)->data[i][0][k + 1] =
00358 ((*out)->data[i][0][(*out)->npix - k] -
00359 (*out)->data[i][0][k]) / 2.0);
00360 (*out)->data[i][0][(*out)->npix - k] = ((*out)->data[i][0][k] =
00361 ((*out)->
00362 data[i][0][(*out)->
00363 npix - k] +
00364 (*out)->
00365 data[i][0][k]) / 2.0);
00366 (*out)->data[i][(*out)->nlin / 2][(*out)->npix - k + 1] =
00367 -((*out)->data[i][(*out)->nlin / 2][k + 1] =
00368 ((*out)->data[i][(*out)->nlin / 2][(*out)->npix - k] -
00369 (*out)->data[i][(*out)->nlin / 2][k]) / 2.0);
00370 (*out)->data[i][(*out)->nlin / 2][(*out)->npix - k] =
00371 ((*out)->data[i][(*out)->nlin / 2][k] =
00372 ((*out)->data[i][(*out)->nlin / 2][(*out)->npix - k] +
00373 (*out)->data[i][(*out)->nlin / 2][k]) / 2.0);
00374 }
00375 for (j = 1; j < (*out)->nlin / 2; j++)
00376 {
00377 (*out)->data[i][(*out)->nlin - j][1] = -((*out)->data[i][j][1] =
00378 ((*out)->
00379 data[i][(*out)->nlin -
00380 j][0] -
00381 (*out)->
00382 data[i][j][0]) / 2.0);
00383 (*out)->data[i][(*out)->nlin - j][0] = ((*out)->data[i][j][0] =
00384 ((*out)->
00385 data[i][(*out)->nlin -
00386 j][0] +
00387 (*out)->
00388 data[i][j][0]) / 2.0);
00389 (*out)->data[i][(*out)->nlin - j][(*out)->npix / 2 + 1] =
00390 -((*out)->data[i][j][(*out)->npix / 2 + 1] =
00391 ((*out)->data[i][(*out)->nlin - j][(*out)->npix / 2] -
00392 (*out)->data[i][j][(*out)->npix / 2]) / 2.0);
00393 (*out)->data[i][(*out)->nlin - j][(*out)->npix / 2] =
00394 ((*out)->data[i][j][(*out)->npix / 2] =
00395 ((*out)->data[i][(*out)->nlin - j][(*out)->npix / 2] +
00396 (*out)->data[i][j][(*out)->npix / 2]) / 2.0);
00397 for (k = 2; k < (*out)->npix / 2; k += 2)
00398 {
00399 t = (*out)->data[i][(*out)->nlin - j][(*out)->npix - k];
00400 u = (*out)->data[i][j][k];
00401 v = (*out)->data[i][j][(*out)->npix - k];
00402 w = (*out)->data[i][(*out)->nlin - j][k];
00403 (*out)->data[i][(*out)->nlin - j][(*out)->npix - k + 1] =
00404 -((*out)->data[i][j][k + 1] = (t - u) / 2.0);
00405 (*out)->data[i][j][(*out)->npix - k] =
00406 ((*out)->data[i][(*out)->nlin - j][k] = (t + u) / 2.0);
00407 (*out)->data[i][j][(*out)->npix - k + 1] =
00408 -((*out)->data[i][(*out)->nlin - j][k + 1] =
00409 (v - w) / 2.0);
00410 (*out)->data[i][(*out)->nlin - j][(*out)->npix - k] =
00411 ((*out)->data[i][j][k] = (v + w) / 2.0);
00412 }
00413 }
00414 }
00415 else
00416 {
00417
00418 for (j = 1; j < (*out)->nlin / 2; j++)
00419 {
00420 for (k = 1; k < (*out)->npix / 2; k++)
00421 {
00422 t =
00423 ((*out)->data[i][j][k] -
00424 (*out)->data[i][(*out)->nlin - j][k] -
00425 (*out)->data[i][j][(*out)->npix - k] +
00426 (*out)->data[i][(*out)->nlin - j][(*out)->npix -
00427 k]) / 2.0;
00428 (*out)->data[i][j][k] -= t;
00429 (*out)->data[i][(*out)->nlin - j][k] += t;
00430 (*out)->data[i][j][(*out)->npix - k] += t;
00431 (*out)->data[i][(*out)->nlin - j][(*out)->npix - k] -= t;
00432 }
00433 }
00434 }
00435
00436
00437 for (j = 0; j < (*out)->nlin / 2; j++)
00438 {
00439 for (k = 0; k < (*out)->npix / 2; k++)
00440 {
00441 t = (*out)->data[i][j][k];
00442 (*out)->data[i][j][k] =
00443 (*out)->data[i][(*out)->nlin / 2 + j][(*out)->npix / 2 + k];
00444 (*out)->data[i][(*out)->nlin / 2 + j][(*out)->npix / 2 + k] = t;
00445 t = (*out)->data[i][(*out)->nlin / 2 + j][k];
00446 (*out)->data[i][(*out)->nlin / 2 + j][k] =
00447 (*out)->data[i][j][(*out)->npix / 2 + k];
00448 (*out)->data[i][j][(*out)->npix / 2 + k] = t;
00449 }
00450 }
00451 }
00452
00453 the_end:
00454 if (a)
00455 free (a);
00456 if (pwr)
00457 free (pwr);
00458 if (sine)
00459 free (sine);
00460 if (tang)
00461 free (tang);
00462 if (array)
00463 free (array);
00464 if (TIMES)
00465 TIMING (T_EXIT);
00466 }