0
|
1 #include <stdio.h>
|
|
2 #include <stdlib.h>
|
|
3 #include <string.h>
|
|
4 #include <math.h>
|
|
5 #include <float.h>
|
|
6 #include "pgm.h"
|
|
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 {
|
|
153 int ii, stepsize, stage, curptr, nthreads, thread, step, nsteps;
|
|
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 {
|
|
185 int ii, stepsize, stage, curptr, nthreads, thread, step, nsteps;
|
|
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
|