annotate Fotopoulos/common.c @ 18:3bdb67e76858

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

Repositories maintained by Peter Meerwald, pmeerw@pmeerw.net.