Mercurial > hg > wm
annotate Fotopoulos-dir/common.c @ 24:9f20bce6184e v0.7
move directories, support netpbm 11
| author | Peter Meerwald-Stadler <pmeerw@pmeerw.net> |
|---|---|
| date | Fri, 20 Dec 2024 13:08:59 +0100 |
| parents | Fotopoulos/common.c@cbecc570129d |
| children |
| rev | line source |
|---|---|
| 0 | 1 #include <stdio.h> |
| 2 #include <stdlib.h> | |
| 3 #include <string.h> | |
| 4 #include <math.h> | |
| 5 #include <float.h> | |
|
24
9f20bce6184e
move directories, support netpbm 11
Peter Meerwald-Stadler <pmeerw@pmeerw.net>
parents:
13
diff
changeset
|
6 #include "netpbm/pgm.h" |
| 0 | 7 #include "common.h" |
| 8 | |
| 9 #define IA 16807 | |
| 10 #define IM 2147483647 | |
| 11 #define AM (1.0/IM) | |
| 12 #define IQ 127773 | |
| 13 #define IR 2836 | |
| 14 #define MASK 123459876 | |
| 15 | |
| 16 static int NN = 0; | |
| 17 static int m = 0; | |
| 18 static double two_over_N = 0; | |
| 19 static double root2_over_rootN = 0; | |
| 20 static double *C = NULL; | |
| 21 | |
| 22 static gray maxval; | |
| 23 static int format; | |
| 24 | |
| 25 void open_image(FILE *in, int *width, int *height) | |
| 26 { | |
| 27 | |
| 28 pgm_readpgminit(in, width, height, &maxval, &format); | |
| 29 } | |
| 30 | |
| 31 void load_image(int **im, FILE *in, int width, int height) | |
| 32 { | |
| 33 int col, row; | |
| 34 gray *rowbuf; | |
| 35 | |
| 36 rowbuf = malloc(sizeof(gray) * width); | |
| 37 | |
| 38 for (row = 0; row < height; row++) { | |
| 39 pgm_readpgmrow(in, rowbuf, width, maxval, format); | |
| 40 for (col = 0; col < width; col++) | |
| 41 im[row][col] = rowbuf[col]; | |
| 42 } | |
| 43 | |
| 44 free(rowbuf); | |
| 45 } | |
| 46 | |
| 47 void save_image(int **im, FILE *out, int width, int height) | |
| 48 { | |
| 49 int col, row; | |
| 50 gray *rowbuf; | |
| 51 | |
| 52 pgm_writepgminit(out, width, height, 255, 0); | |
| 53 | |
| 54 rowbuf = malloc(sizeof(gray) * width); | |
| 55 | |
| 56 for (row = 0; row < height; row++) { | |
| 57 for (col = 0; col < width; col++) | |
| 58 rowbuf[col] = im[row][col]; | |
| 59 pgm_writepgmrow(out, rowbuf, width, 255, 0); | |
| 60 } | |
| 61 | |
| 62 free(rowbuf); | |
| 63 } | |
| 64 | |
| 65 int ** imatrix(int nrows, int ncols) | |
| 66 { | |
| 67 int **m; | |
| 68 int i, j; | |
| 69 m = (int **) malloc (nrows * sizeof(int *)); | |
| 70 for (i = 0; i < nrows; i++) { | |
| 71 m[i] = (int *) malloc (ncols * sizeof(int)); | |
| 72 if (!m[i]) fprintf(stderr, "\nIt's not working"); | |
| 73 } | |
| 74 for (i = 0; i < nrows; i++) | |
| 75 for (j = 0; j < ncols; j++) | |
| 76 m[i][j] = 0; | |
| 77 return m; | |
| 78 } | |
| 79 | |
| 80 void freematrix(int **I, int rows) | |
| 81 { | |
| 82 int k; | |
| 83 for (k = 0; k < rows; k++) | |
| 84 free (I[k]); | |
| 85 } | |
| 86 | |
| 87 float ran0(long int *idum) | |
| 88 { | |
| 89 long int k; | |
| 90 float ans; | |
| 91 *idum ^= MASK; | |
| 92 k = (*idum) / IQ; | |
| 93 *idum = IA * (*idum - k * IQ) - IR * k; | |
| 94 if (*idum < 0) *idum += IM; | |
| 95 ans = AM * (*idum); | |
| 96 *idum ^= MASK; | |
| 97 return ans; | |
| 98 } | |
| 99 | |
| 100 float gasdev(long int *idum) | |
| 101 { | |
| 102 | |
| 103 float v1; | |
| 104 v1 = (float) sqrt( -2.0 * log(ran0(idum))) * cos(2 * PI * ran0(idum)); | |
| 105 return v1; | |
| 106 } | |
| 107 | |
| 108 | |
| 109 void put_image_from_int_2_double(int **i, double *f, int N) | |
| 110 { | |
| 111 int l, j, k; | |
| 112 k = 0; | |
| 113 for (l = 0; l < N; l++) | |
| 114 for (j = 0; j < N; j++) | |
| 115 f[k++] = (double) i[l][j]; | |
| 116 } | |
| 117 | |
| 118 void put_image_from_double_2_int(double *f, int **i, int N) | |
| 119 { | |
| 120 int l, j, k; | |
| 121 k = 0; | |
| 122 for (l = 0; l < N; l++) | |
| 123 for (j = 0; j < N; j++) | |
| 124 i[l][j] = (int) f[k++]; | |
| 125 } | |
| 126 | |
| 127 | |
| 128 void bitrev(double *f, int len) | |
| 129 { | |
| 130 int i, j, m, halflen; | |
| 131 double temp; | |
| 132 | |
| 133 if (len <= 2) return ; /* No action necessary if n=1 or n=2 */ | |
| 134 halflen = len >> 1; | |
| 135 j = 1; | |
| 136 for (i = 1; i <= len; i++) { | |
| 137 if (i < j) { | |
| 138 temp = f[j - 1]; | |
| 139 f[j - 1] = f[i - 1]; | |
| 140 f[i - 1] = temp; | |
| 141 } | |
| 142 m = halflen; | |
| 143 while (j > m) { | |
| 144 j = j - m; | |
| 145 m = (m + 1) >> 1; | |
| 146 } | |
| 147 j = j + m; | |
| 148 } | |
| 149 } | |
| 150 | |
| 151 void inv_sums(double *f) | |
| 152 { | |
| 13 | 153 int stepsize, stage, curptr, nthreads, thread, step, nsteps; |
| 0 | 154 |
| 155 for (stage = 1; stage <= m - 1; stage++) { | |
| 156 nthreads = 1 << (stage - 1); | |
| 157 stepsize = nthreads << 1; | |
| 158 nsteps = (1 << (m - stage)) - 1; | |
| 159 for (thread = 1; thread <= nthreads; thread++) { | |
| 160 curptr = NN - thread; | |
| 161 for (step = 1; step <= nsteps; step++) { | |
| 162 f[curptr] += f[curptr - stepsize]; | |
| 163 curptr -= stepsize; | |
| 164 } | |
| 165 } | |
| 166 } | |
| 167 } | |
| 168 | |
| 169 //-------------------------------------------------------- | |
| 170 //Foreign code - FCT from Bath Univerity | |
| 171 //-------------------------------------------------------- | |
| 172 void rarrwrt(double f[], int n) | |
| 173 { | |
| 174 int i; | |
| 175 | |
| 176 for (i = 0; i <= n - 1; i++) { | |
| 177 fprintf(stderr, "%4d : %f\n", i, f[i]); | |
| 178 } | |
| 179 } | |
| 180 | |
| 181 /* fast DCT based on IEEE signal proc, 1992 #8, yugoslavian authors. */ | |
| 182 | |
| 183 void fwd_sums(double *f) | |
| 184 { | |
| 13 | 185 int stepsize, stage, curptr, nthreads, thread, step, nsteps; |
| 0 | 186 |
| 187 for (stage = m - 1; stage >= 1; stage--) { | |
| 188 nthreads = 1 << (stage - 1); | |
| 189 stepsize = nthreads << 1; | |
| 190 nsteps = (1 << (m - stage)) - 1; | |
| 191 for (thread = 1; thread <= nthreads; thread++) { | |
| 192 curptr = nthreads + thread - 1; | |
| 193 for (step = 1; step <= nsteps; step++) { | |
| 194 f[curptr] += f[curptr + stepsize]; | |
| 195 curptr += stepsize; | |
| 196 } | |
| 197 } | |
| 198 } | |
| 199 } | |
| 200 | |
| 201 void scramble(double *f, int len) | |
| 202 { | |
| 203 double temp; | |
| 204 int i, ii1, ii2, halflen, qtrlen; | |
| 205 | |
| 206 halflen = len >> 1; | |
| 207 qtrlen = halflen >> 1; | |
| 208 bitrev(f, len); | |
| 209 bitrev(&f[0], halflen); | |
| 210 bitrev(&f[halflen], halflen); | |
| 211 ii1 = len - 1; | |
| 212 ii2 = halflen; | |
| 213 for (i = 0; i <= qtrlen - 1; i++) { | |
| 214 temp = f[ii1]; | |
| 215 f[ii1] = f[ii2]; | |
| 216 f[ii2] = temp; | |
| 217 ii1--; | |
| 218 ii2++; | |
| 219 } | |
| 220 } | |
| 221 | |
| 222 void unscramble(double *f, int len) | |
| 223 { | |
| 224 double temp; | |
| 225 int i, ii1, ii2, halflen, qtrlen; | |
| 226 | |
| 227 halflen = len >> 1; | |
| 228 qtrlen = halflen >> 1; | |
| 229 ii1 = len - 1; | |
| 230 ii2 = halflen; | |
| 231 for (i = 0; i <= qtrlen - 1; i++) { | |
| 232 temp = f[ii1]; | |
| 233 f[ii1] = f[ii2]; | |
| 234 f[ii2] = temp; | |
| 235 ii1--; | |
| 236 ii2++; | |
| 237 } | |
| 238 bitrev(&f[0], halflen); | |
| 239 bitrev(&f[halflen], halflen); | |
| 240 bitrev(f, len); | |
| 241 } | |
| 242 | |
| 243 void initcosarray(int length) | |
| 244 { | |
| 245 int i, group, base, item, nitems, halfN; | |
| 246 double factor; | |
| 247 | |
| 248 m = -1; | |
| 249 do { | |
| 250 m++; | |
| 251 NN = 1 << m; | |
| 252 if (NN > length) { | |
| 253 fprintf(stderr, "ERROR in FCT-- length %d not a power of 2\n", length); | |
| 254 exit(1); | |
| 255 } | |
| 256 } while (NN < length); | |
| 257 if (C != NULL) free(C); | |
| 258 C = (double *)calloc(NN, sizeof(double)); | |
| 259 if (C == NULL) { | |
| 260 fprintf(stderr, "Unable to allocate C array\n"); | |
| 261 exit(1); | |
| 262 } | |
| 263 halfN = NN / 2; | |
| 264 two_over_N = 2.0 / (double)NN; | |
| 265 root2_over_rootN = sqrt(2.0 / (double)NN); | |
| 266 for (i = 0; i <= halfN - 1; i++) C[halfN + i] = 4 * i + 1; | |
| 267 for (group = 1; group <= m - 1; group++) { | |
| 268 base = 1 << (group - 1); | |
| 269 nitems = base; | |
| 270 factor = 1.0 * (1 << (m - group)); | |
| 271 for (item = 1; item <= nitems; item++) C[base + item - 1] = factor * C[halfN + item - 1]; | |
| 272 } | |
| 273 | |
| 274 //printf("before taking cos, C array =\n"); rarrwrt(C,N); | |
| 275 for (i = 1; i <= NN - 1; i++) C[i] = 1.0 / (2.0 * cos(C[i] * PI / (2.0 * NN))); | |
| 276 //printf("After taking cos, Carray = \n"); rarrwrt(C,N); | |
| 277 } | |
| 278 | |
| 279 | |
| 280 void inv_butterflies(double *f) | |
| 281 { | |
| 282 int stage, ii1, ii2, butterfly, ngroups, group, wingspan, increment, baseptr; | |
| 283 double Cfac, T; | |
| 284 | |
| 285 for (stage = 1; stage <= m; stage++) { | |
| 286 ngroups = 1 << (m - stage); | |
| 287 wingspan = 1 << (stage - 1); | |
| 288 increment = wingspan << 1; | |
| 289 for (butterfly = 1; butterfly <= wingspan; butterfly++) { | |
| 290 Cfac = C[wingspan + butterfly - 1]; | |
| 291 baseptr = 0; | |
| 292 for (group = 1; group <= ngroups; group++) { | |
| 293 ii1 = baseptr + butterfly - 1; | |
| 294 ii2 = ii1 + wingspan; | |
| 295 T = Cfac * f[ii2]; | |
| 296 f[ii2] = f[ii1]-T; | |
| 297 f[ii1] = f[ii1] + T; | |
| 298 baseptr += increment; | |
| 299 } | |
| 300 } | |
| 301 } | |
| 302 } | |
| 303 | |
| 304 void fwd_butterflies(double *f) | |
| 305 { | |
| 306 int stage, ii1, ii2, butterfly, ngroups, group, wingspan, increment, baseptr; | |
| 307 double Cfac, T; | |
| 308 | |
| 309 for (stage = m; stage >= 1; stage--) { | |
| 310 ngroups = 1 << (m - stage); | |
| 311 wingspan = 1 << (stage - 1); | |
| 312 increment = wingspan << 1; | |
| 313 for (butterfly = 1; butterfly <= wingspan; butterfly++) { | |
| 314 Cfac = C[wingspan + butterfly - 1]; | |
| 315 baseptr = 0; | |
| 316 for (group = 1; group <= ngroups; group++) { | |
| 317 ii1 = baseptr + butterfly - 1; | |
| 318 ii2 = ii1 + wingspan; | |
| 319 T = f[ii2]; | |
| 320 f[ii2] = Cfac * (f[ii1]-T); | |
| 321 f[ii1] = f[ii1] + T; | |
| 322 baseptr += increment; | |
| 323 } | |
| 324 } | |
| 325 } | |
| 326 } | |
| 327 | |
| 328 void ifct_noscale(double *f, int length) | |
| 329 { | |
| 330 if (length != NN) initcosarray(length); | |
| 331 f[0] *= INVROOT2; | |
| 332 inv_sums(f); | |
| 333 bitrev(f, NN); | |
| 334 inv_butterflies(f); | |
| 335 unscramble(f, NN); | |
| 336 } | |
| 337 | |
| 338 void fct_noscale(double *f, int length) | |
| 339 { | |
| 340 if (length != NN) initcosarray(length); | |
| 341 scramble(f, NN); | |
| 342 fwd_butterflies(f); | |
| 343 bitrev(f, NN); | |
| 344 fwd_sums(f); | |
| 345 f[0] *= INVROOT2; | |
| 346 } | |
| 347 | |
| 348 void ifct_defn_scaling(double *f, int length) | |
| 349 { | |
| 350 ifct_noscale(f, length); | |
| 351 } | |
| 352 | |
| 353 void fct_defn_scaling(double *f, int length) | |
| 354 { | |
| 355 int i; | |
| 356 | |
| 357 fct_noscale(f, length); | |
| 358 for (i = 0; i <= NN - 1; i++) f[i] *= two_over_N; | |
| 359 } | |
| 360 | |
| 361 void ifct(double *f, int length) | |
| 362 { | |
| 363 /* CALL THIS FOR INVERSE 1D DCT DON-MONRO PREFERRED SCALING */ | |
| 364 int i; | |
| 365 | |
| 366 if (length != NN) initcosarray(length); /* BGS patch June 1997 */ | |
| 367 for (i = 0; i <= NN - 1; i++) f[i] *= root2_over_rootN; | |
| 368 ifct_noscale(f, length); | |
| 369 } | |
| 370 | |
| 371 void fct(double *f, int length) | |
| 372 { | |
| 373 /* CALL THIS FOR FORWARD 1D DCT DON-MONRO PREFERRED SCALING */ | |
| 374 int i; | |
| 375 | |
| 376 fct_noscale(f, length); | |
| 377 for (i = 0; i <= NN - 1; i++) f[i] *= root2_over_rootN; | |
| 378 } | |
| 379 | |
| 380 /**************************************************************** | |
| 381 2D FAST DCT SECTION | |
| 382 ****************************************************************/ | |
| 383 | |
| 384 static double *g = NULL; | |
| 385 static double two_over_sqrtncolsnrows = 0.0; | |
| 386 static int ncolsvalue = 0; | |
| 387 static int nrowsvalue = 0; | |
| 388 | |
| 389 void initfct2d(int nrows, int ncols) | |
| 390 { | |
| 391 if ((nrows <= 0) || (ncols < 0)) { | |
| 392 fprintf(stderr, "FCT2D -- ncols=%d or nrows=%d is <=0\n", nrows, ncols); | |
| 393 exit(1); | |
| 394 } | |
| 395 if (g != NULL) free(g); | |
| 396 g = (double *)calloc(nrows, sizeof(double)); | |
| 397 if (g == NULL) { | |
| 398 fprintf(stderr, "FCT2D -- Unable to allocate g array\n"); | |
| 399 exit(1); | |
| 400 } | |
| 401 ncolsvalue = ncols; | |
| 402 nrowsvalue = nrows; | |
| 403 two_over_sqrtncolsnrows = 2.0 / sqrt(ncols * 1.0 * nrows); | |
| 404 } | |
| 405 | |
| 406 void fct2d(double f[], int nrows, int ncols) | |
| 407 /* CALL THIS FOR FORWARD 2d DCT DON-MONRO PREFERRED SCALING */ | |
| 408 { | |
| 409 int u, v; | |
| 410 | |
| 411 if ((ncols != ncolsvalue) || (nrows != nrowsvalue)){ | |
| 412 initfct2d(nrows, ncols); | |
| 413 } | |
| 414 for (u = 0; u <= nrows - 1; u++){ | |
| 415 fct_noscale(&f[u*ncols], ncols); | |
| 416 } | |
| 417 for (v = 0; v <= ncols - 1; v++){ | |
| 418 for (u = 0; u <= nrows - 1; u++) { | |
| 419 g[u] = f[u * ncols + v]; | |
| 420 } | |
| 421 fct_noscale(g, nrows); | |
| 422 for (u = 0; u <= nrows - 1; u++) { | |
| 423 f[u*ncols + v] = g[u] * two_over_sqrtncolsnrows; | |
| 424 } | |
| 425 } | |
| 426 } | |
| 427 | |
| 428 void ifct2d(double f[], int nrows, int ncols) | |
| 429 /* CALL THIS FOR INVERSE 2d DCT DON-MONRO PREFERRED SCALING */ | |
| 430 { | |
| 431 int u, v; | |
| 432 | |
| 433 if ((ncols != ncolsvalue) || (nrows != nrowsvalue)){ | |
| 434 initfct2d(nrows, ncols); | |
| 435 } | |
| 436 for (u = 0; u <= nrows - 1; u++){ | |
| 437 ifct_noscale(&f[u*ncols], ncols); | |
| 438 } | |
| 439 for (v = 0; v <= ncols - 1; v++){ | |
| 440 for (u = 0; u <= nrows - 1; u++) { | |
| 441 g[u] = f[u * ncols + v]; | |
| 442 } | |
| 443 ifct_noscale(g, nrows); | |
| 444 for (u = 0; u <= nrows - 1; u++) { | |
| 445 f[u*ncols + v] = g[u] * two_over_sqrtncolsnrows; | |
| 446 } | |
| 447 } | |
| 448 } | |
| 449 | |
| 450 void matmul(double **a, double **b, double **r, int N) | |
| 451 { | |
| 452 int i, j, k; | |
| 453 | |
| 454 for (i = 0; i < N; i++) | |
| 455 for (j = 0; j < N; j++) { | |
| 456 r[i][j] = 0.0; | |
| 457 for (k = 0; k < N; k++) | |
| 458 r[i][j] += a[i][k] * b[k][j]; | |
| 459 } | |
| 460 } | |
| 461 | |
| 462 void hartley(double **in, double **out, int N) | |
| 463 { | |
| 464 int k, n; | |
| 465 double **h; | |
| 466 | |
| 467 h = dmatrix(N, N); | |
| 468 //Building up the transformation matrix | |
| 469 for (k = 0; k < N; k++) | |
| 470 for (n = 0; n < N; n++) | |
| 471 h[k][n] = (cos(2 * PI * k * n / N) + sin(2 * PI * k * n / N)) / sqrt(N); | |
| 472 | |
| 473 // Now we have to multiply the input with the transformation matrix | |
| 474 matmul(h, in, out, N); | |
| 475 freematrix_d(h, N); | |
| 476 return ; | |
| 477 } | |
| 478 | |
| 479 double ** dmatrix(int nrows, int ncols) | |
| 480 { | |
| 481 double **m; | |
| 482 int i, j; | |
| 483 m = (double **) malloc (nrows * sizeof(double *)); | |
| 484 for (i = 0; i < nrows; i++) { | |
| 485 m[i] = (double *) malloc (ncols * sizeof(double)); | |
| 486 if (!m[i]) printf("\nIt's not working"); | |
| 487 } | |
| 488 for (i = 0; i < nrows; i++) | |
| 489 for (j = 0; j < ncols; j++) | |
| 490 m[i][j] = 0.0; | |
| 491 return m; | |
| 492 } | |
| 493 | |
| 494 void freematrix_d(double **I, int rows) | |
| 495 { | |
| 496 int k; | |
| 497 for (k = 0; k < rows; k++) | |
| 498 free (I[k]); | |
| 499 } | |
| 500 | |
| 501 void matrix_i2d(int **i, double **d, int N) | |
| 502 { | |
| 503 int x, y; | |
| 504 for (x = 0; x < N; x++) | |
| 505 for (y = 0; y < N; y++) | |
| 506 d[y][x] = i[y][x]; | |
| 507 } | |
| 508 | |
| 509 void matrix_d2i(double **d, int **i, int N) | |
| 510 { | |
| 511 int x, y; | |
| 512 for (x = 0; x < N; x++) | |
| 513 for (y = 0; y < N; y++) | |
| 514 i[y][x] = d[y][x]; | |
| 515 } | |
| 516 | |
| 517 double * dvector(long int N) | |
| 518 { | |
| 519 double *m; | |
| 520 m = (double *) malloc (N * sizeof(double)); | |
| 521 if (!m) printf("\nIt's not working"); | |
| 522 return m; | |
| 523 } | |
| 524 | |
| 525 void put_matrix_2_vector(double **i, double *f, int N) | |
| 526 { | |
| 527 int l, j, k; | |
| 528 k = 0; | |
| 529 for (l = 0; l < N; l++) | |
| 530 for (j = 0; j < N; j++) | |
| 531 f[k++] = i[l][j]; | |
| 532 } | |
| 533 | |
| 534 void put_vector_2_matrix(double *f, double **i, int N) | |
| 535 { | |
| 536 int l, j, k; | |
| 537 k = 0; | |
| 538 for (l = 0; l < N; l++) | |
| 539 for (j = 0; j < N; j++) | |
| 540 i[l][j] = f[k++]; | |
| 541 } | |
| 542 | |
| 543 void set_in_binary() { | |
| 544 #if defined(EMX) | |
| 545 _fsetmode(in, "b"); | |
| 546 #elif defined(MINGW) | |
| 547 setmode(STDIN_FILENO, O_BINARY); | |
| 548 #endif | |
| 549 } | |
| 550 | |
| 551 void set_out_binary() { | |
| 552 #if defined(EMX) | |
| 553 _fsetmode(out, "b"); | |
| 554 #elif defined(MINGW) | |
| 555 setmode(STDOUT_FILENO, O_BINARY); | |
| 556 #endif | |
| 557 } | |
| 558 | |
| 559 void wm_init2() { | |
| 560 set_out_binary(); | |
| 561 } | |
| 562 | |
| 563 void wm_init1() { | |
| 564 set_in_binary(); | |
| 565 } | |
| 566 | |
| 567 void wm_init() { | |
| 568 set_in_binary(); | |
| 569 set_out_binary(); | |
| 570 } | |
| 571 |