|
0
|
1 #include "wm.h"
|
|
|
2 #include "dct.h"
|
|
|
3
|
|
|
4 #define INVROOT2 0.7071067814
|
|
|
5 #define SWAP(A, B) {double t = A; A = B; B = t;}
|
|
|
6
|
|
|
7 int N;
|
|
|
8 int M;
|
|
|
9
|
|
|
10 double *dct_NxN_tmp = NULL;
|
|
|
11 double *dct_NxN_costable = NULL;
|
|
|
12 int dct_NxN_log2N = 0;
|
|
|
13
|
|
|
14 static const unsigned int JPEG_lumin_quant_table[NJPEG][NJPEG] = {
|
|
|
15 {16, 11, 10, 16, 24, 40, 51, 61},
|
|
|
16 {12, 12, 14, 19, 26, 58, 60, 55},
|
|
|
17 {14, 13, 16, 24, 40, 57, 69, 56},
|
|
|
18 {14, 17, 22, 29, 51, 87, 80, 62},
|
|
|
19 {18, 22, 37, 56, 68, 109, 103, 77},
|
|
|
20 {24, 35, 55, 64, 81, 104, 113, 92},
|
|
|
21 {49, 64, 78, 87, 103, 121, 120, 101},
|
|
|
22 {72, 92, 95, 98, 112, 100, 103, 99}};
|
|
|
23
|
|
|
24 static void initcosarray()
|
|
|
25 {
|
|
|
26 int i,group,base,item,nitems,halfN;
|
|
|
27 double factor;
|
|
|
28
|
|
|
29 dct_NxN_log2N = -1;
|
|
|
30 do{
|
|
|
31 dct_NxN_log2N++;
|
|
|
32 if ((1<<dct_NxN_log2N)>N){
|
|
|
33 fprintf(stderr, "dct_NxN: %d not a power of 2\n", N);
|
|
|
34 exit(1);
|
|
|
35 }
|
|
|
36 }while((1<<dct_NxN_log2N)<N);
|
|
|
37 if (dct_NxN_costable) free(dct_NxN_costable);
|
|
3
|
38 dct_NxN_costable = malloc(N * sizeof(double));
|
|
0
|
39 #ifdef DEBUG
|
|
|
40 if(!dct_NxN_costable){
|
|
|
41 fprintf(stderr, "Unable to allocate C array\n");
|
|
|
42 exit(1);
|
|
|
43 }
|
|
|
44 #endif
|
|
|
45 halfN=N/2;
|
|
|
46 for(i=0;i<=halfN-1;i++) dct_NxN_costable[halfN+i]=4*i+1;
|
|
|
47 for(group=1;group<=dct_NxN_log2N-1;group++){
|
|
|
48 base= 1<<(group-1);
|
|
|
49 nitems=base;
|
|
|
50 factor = 1.0*(1<<(dct_NxN_log2N-group));
|
|
|
51 for(item=1; item<=nitems;item++) dct_NxN_costable[base+item-1]=factor*dct_NxN_costable[halfN+item-1];
|
|
|
52 }
|
|
|
53
|
|
|
54 for(i=1;i<=N-1;i++) dct_NxN_costable[i] = 1.0/(2.0*cos(dct_NxN_costable[i]*M_PI/(2.0*N)));
|
|
|
55 }
|
|
|
56
|
|
|
57 void init_dct_NxN(int width, int height) {
|
|
|
58 #ifdef DEBUG
|
|
|
59 if (width != height || width <= 0) {
|
|
|
60 fprintf(stderr, "init_dct_NxN(): dimensions out of range\n");
|
|
|
61 exit(1);
|
|
|
62 }
|
|
|
63 #endif
|
|
|
64
|
|
|
65 if (dct_NxN_tmp && M != height)
|
|
|
66 free(dct_NxN_tmp);
|
|
|
67
|
|
|
68 N = width;
|
|
|
69 M = height;
|
|
|
70
|
|
3
|
71 dct_NxN_tmp = malloc(height * sizeof(double));
|
|
0
|
72 #ifdef DEBUG
|
|
|
73 if (!dct_NxN_tmp) {
|
|
|
74 fprintf(stderr, "init_dct_NxN(): failed to allocate memory\n");
|
|
|
75 exit(1);
|
|
|
76 }
|
|
|
77 #endif
|
|
|
78
|
|
|
79 initcosarray();
|
|
|
80 }
|
|
|
81
|
|
|
82 static void bitrev(double *f, int len)
|
|
|
83 {
|
|
|
84 int i,j,m;
|
|
|
85
|
|
|
86 if (len<=2) return; /* No action necessary if n=1 or n=2 */
|
|
|
87 j=1;
|
|
|
88 for(i=1; i<=len; i++){
|
|
|
89 if(i<j)
|
|
|
90 SWAP(f[j-1], f[i-1]);
|
|
|
91 m = len>>1;
|
|
|
92 while(j>m){
|
|
|
93 j=j-m;
|
|
|
94 m=(m+1)>>1;
|
|
|
95 }
|
|
|
96 j=j+m;
|
|
|
97 }
|
|
|
98 }
|
|
|
99
|
|
|
100 static void inv_sums(double *f)
|
|
|
101 {
|
|
|
102 int stepsize,stage,curptr,nthreads,thread,step,nsteps;
|
|
|
103
|
|
|
104 for(stage=1; stage <=dct_NxN_log2N-1; stage++){
|
|
|
105 nthreads = 1<<(stage-1);
|
|
|
106 stepsize = nthreads<<1;
|
|
|
107 nsteps = (1<<(dct_NxN_log2N-stage)) - 1;
|
|
|
108 for(thread=1; thread<=nthreads; thread++){
|
|
|
109 curptr=N-thread;
|
|
|
110 for(step=1; step<=nsteps; step++){
|
|
|
111 f[curptr] += f[curptr-stepsize];
|
|
|
112 curptr -= stepsize;
|
|
|
113 }
|
|
|
114 }
|
|
|
115 }
|
|
|
116 }
|
|
|
117
|
|
|
118 static void fwd_sums(double *f)
|
|
|
119 {
|
|
|
120 int stepsize,stage,curptr,nthreads,thread,step,nsteps;
|
|
|
121
|
|
|
122 for(stage=dct_NxN_log2N-1; stage >=1; stage--){
|
|
|
123 nthreads = 1<<(stage-1);
|
|
|
124 stepsize = nthreads<<1;
|
|
|
125 nsteps = (1<<(dct_NxN_log2N-stage)) - 1;
|
|
|
126 for(thread=1; thread<=nthreads; thread++){
|
|
|
127 curptr=nthreads +thread-1;
|
|
|
128 for(step=1; step<=nsteps; step++){
|
|
|
129 f[curptr] += f[curptr+stepsize];
|
|
|
130 curptr += stepsize;
|
|
|
131 }
|
|
|
132 }
|
|
|
133 }
|
|
|
134 }
|
|
|
135
|
|
|
136 static void scramble(double *f,int len){
|
|
|
137 int i,ii1,ii2;
|
|
|
138
|
|
|
139 bitrev(f,len);
|
|
|
140 bitrev(&f[0], len>>1);
|
|
|
141 bitrev(&f[len>>1], len>>1);
|
|
|
142 ii1=len-1;
|
|
|
143 ii2=len>>1;
|
|
|
144 for(i=0; i<(len>>2); i++){
|
|
|
145 SWAP(f[ii1], f[ii2]);
|
|
|
146 ii1--;
|
|
|
147 ii2++;
|
|
|
148 }
|
|
|
149 }
|
|
|
150
|
|
|
151 static void unscramble(double *f,int len)
|
|
|
152 {
|
|
|
153 int i,ii1,ii2;
|
|
|
154
|
|
|
155 ii1 = len-1;
|
|
|
156 ii2 = len>>1;
|
|
|
157 for(i=0; i<(len>>2); i++){
|
|
|
158 SWAP(f[ii1], f[ii2]);
|
|
|
159 ii1--;
|
|
|
160 ii2++;
|
|
|
161 }
|
|
|
162 bitrev(&f[0], len>>1);
|
|
|
163 bitrev(&f[len>>1], len>>1);
|
|
|
164 bitrev(f,len);
|
|
|
165 }
|
|
|
166
|
|
|
167 static void inv_butterflies(double *f)
|
|
|
168 {
|
|
|
169 int stage,ii1,ii2,butterfly,ngroups,group,wingspan,increment,baseptr;
|
|
|
170 double Cfac,T;
|
|
|
171
|
|
|
172 for(stage=1; stage<=dct_NxN_log2N;stage++){
|
|
|
173 ngroups=1<<(dct_NxN_log2N-stage);
|
|
|
174 wingspan=1<<(stage-1);
|
|
|
175 increment=wingspan<<1;
|
|
|
176 for(butterfly=1; butterfly<=wingspan; butterfly++){
|
|
|
177 Cfac = dct_NxN_costable[wingspan+butterfly-1];
|
|
|
178 baseptr=0;
|
|
|
179 for(group=1; group<=ngroups; group++){
|
|
|
180 ii1=baseptr+butterfly-1;
|
|
|
181 ii2=ii1+wingspan;
|
|
|
182 T=Cfac * f[ii2];
|
|
|
183 f[ii2]=f[ii1]-T;
|
|
|
184 f[ii1]=f[ii1]+T;
|
|
|
185 baseptr += increment;
|
|
|
186 }
|
|
|
187 }
|
|
|
188 }
|
|
|
189 }
|
|
|
190
|
|
|
191 static void fwd_butterflies(double *f)
|
|
|
192 {
|
|
|
193 int stage,ii1,ii2,butterfly,ngroups,group,wingspan,increment,baseptr;
|
|
|
194 double Cfac,T;
|
|
|
195
|
|
|
196 for(stage=dct_NxN_log2N; stage>=1;stage--){
|
|
|
197 ngroups=1<<(dct_NxN_log2N-stage);
|
|
|
198 wingspan=1<<(stage-1);
|
|
|
199 increment=wingspan<<1;
|
|
|
200 for(butterfly=1; butterfly<=wingspan; butterfly++){
|
|
|
201 Cfac = dct_NxN_costable[wingspan+butterfly-1];
|
|
|
202 baseptr=0;
|
|
|
203 for(group=1; group<=ngroups; group++){
|
|
|
204 ii1=baseptr+butterfly-1;
|
|
|
205 ii2=ii1+wingspan;
|
|
|
206 T= f[ii2];
|
|
|
207 f[ii2]=Cfac *(f[ii1]-T);
|
|
|
208 f[ii1]=f[ii1]+T;
|
|
|
209 baseptr += increment;
|
|
|
210 }
|
|
|
211 }
|
|
|
212 }
|
|
|
213 }
|
|
|
214
|
|
|
215 static void ifct_noscale(double *f)
|
|
|
216 {
|
|
|
217 f[0] *= INVROOT2;
|
|
|
218 inv_sums(f);
|
|
|
219 bitrev(f,N);
|
|
|
220 inv_butterflies(f);
|
|
|
221 unscramble(f,N);
|
|
|
222 }
|
|
|
223
|
|
|
224 static void fct_noscale(double *f)
|
|
|
225 {
|
|
|
226 scramble(f,N);
|
|
|
227 fwd_butterflies(f);
|
|
|
228 bitrev(f,N);
|
|
|
229 fwd_sums(f);
|
|
|
230 f[0] *= INVROOT2;
|
|
|
231 }
|
|
|
232
|
|
|
233 void fdct_NxN(gray **pixels, double **dcts) {
|
|
|
234 int u,v;
|
|
|
235 double two_over_sqrtncolsnrows = 2.0/sqrt((double) N*M);
|
|
|
236
|
|
|
237 for (u=0; u < N; u++)
|
|
|
238 for (v=0; v < M; v++)
|
|
3
|
239 dcts[u][v] = ((int) pixels[u][v] - 128);
|
|
0
|
240
|
|
|
241 for (u=0; u<=M-1; u++){
|
|
|
242 fct_noscale(dcts[u]);
|
|
|
243 }
|
|
|
244 for (v=0; v<=N-1; v++){
|
|
|
245 for (u=0; u<=M-1; u++){
|
|
|
246 dct_NxN_tmp[u] = dcts[u][v];
|
|
|
247 }
|
|
|
248 fct_noscale(dct_NxN_tmp);
|
|
|
249 for (u=0; u<=M-1; u++){
|
|
|
250 dcts[u][v] = dct_NxN_tmp[u]*two_over_sqrtncolsnrows;
|
|
|
251 }
|
|
|
252 }
|
|
|
253 }
|
|
|
254
|
|
|
255 void idct_NxN(double **dcts, gray **pixels) {
|
|
|
256 int u,v;
|
|
|
257 double two_over_sqrtncolsnrows = 2.0/sqrt((double) N*M);
|
|
|
258
|
|
|
259 double **tmp;
|
|
|
260
|
|
|
261 tmp = alloc_coeffs(N, N);
|
|
|
262 for (u=0;u<N;u++)
|
|
|
263 for (v=0;v<M;v++)
|
|
|
264 tmp[u][v] = dcts[u][v];
|
|
|
265
|
|
|
266 for (u=0; u<=M-1; u++){
|
|
|
267 ifct_noscale(tmp[u]);
|
|
|
268 }
|
|
|
269 for (v=0; v<=N-1; v++){
|
|
|
270 for (u=0; u<=M-1; u++){
|
|
|
271 dct_NxN_tmp[u] = tmp[u][v];
|
|
|
272 }
|
|
|
273 ifct_noscale(dct_NxN_tmp);
|
|
|
274 for (u=0; u<=M-1; u++){
|
|
|
275 tmp[u][v] = dct_NxN_tmp[u]*two_over_sqrtncolsnrows;
|
|
|
276 }
|
|
|
277 }
|
|
|
278
|
|
|
279 for (u=0;u<N;u++)
|
|
|
280 for (v=0;v<M;v++)
|
|
|
281 pixels[u][v] = PIXELRANGE(tmp[u][v] + 128.5);
|
|
|
282 free(tmp);
|
|
|
283 }
|
|
|
284
|
|
|
285 void fdct_inplace_NxN(double **coeffs) {
|
|
|
286 int u,v;
|
|
|
287 double two_over_sqrtncolsnrows = 2.0/sqrt((double) N*M);
|
|
|
288
|
|
|
289 for (u=0; u<=M-1; u++)
|
|
|
290 fct_noscale(coeffs[u]);
|
|
|
291
|
|
|
292 for (v=0; v<=N-1; v++){
|
|
|
293 for (u=0; u<=M-1; u++)
|
|
|
294 dct_NxN_tmp[u] = coeffs[u][v];
|
|
|
295
|
|
|
296 fct_noscale(dct_NxN_tmp);
|
|
|
297 for (u=0; u<=M-1; u++)
|
|
|
298 coeffs[u][v] = dct_NxN_tmp[u]*two_over_sqrtncolsnrows;
|
|
|
299 }
|
|
|
300 }
|
|
|
301
|
|
|
302 void idct_inplace_NxN(double **coeffs) {
|
|
|
303 int u,v;
|
|
|
304 double two_over_sqrtncolsnrows = 2.0/sqrt((double) N*M);
|
|
|
305
|
|
|
306 for (u=0; u<=M-1; u++)
|
|
|
307 ifct_noscale(coeffs[u]);
|
|
|
308
|
|
|
309 for (v=0; v<=N-1; v++) {
|
|
|
310 for (u=0; u<=M-1; u++)
|
|
|
311 dct_NxN_tmp[u] = coeffs[u][v];
|
|
|
312
|
|
|
313 ifct_noscale(dct_NxN_tmp);
|
|
|
314 for (u=0; u<=M-1; u++)
|
|
|
315 coeffs[u][v] = dct_NxN_tmp[u]*two_over_sqrtncolsnrows;
|
|
|
316 }
|
|
|
317
|
|
|
318 }
|
|
|
319
|
|
|
320 double **dct_NxM_costable_x = NULL;
|
|
|
321 double **dct_NxM_costable_y = NULL;
|
|
|
322
|
|
|
323 void init_dct_NxM(int cols, int rows) {
|
|
|
324 int i, j;
|
|
|
325 double cx = sqrt(2.0 / cols);
|
|
|
326 double cy = sqrt(2.0 / rows);
|
|
|
327
|
|
|
328 #ifdef DEBUG
|
|
|
329 if (cols <= 0 || rows <= 0) {
|
|
|
330 fprintf(stderr, "init_dct_NxM(): dimensions out of range\n");
|
|
|
331 exit(1);
|
|
|
332 }
|
|
|
333 #endif
|
|
|
334
|
|
|
335 if (dct_NxM_costable_x && N != cols) {
|
|
|
336 free_coeffs(dct_NxM_costable_x);
|
|
|
337 dct_NxM_costable_x = NULL;
|
|
|
338 }
|
|
|
339
|
|
|
340 if (dct_NxM_costable_y && M != rows) {
|
|
|
341 free_coeffs(dct_NxM_costable_y);
|
|
|
342 dct_NxM_costable_y = NULL;
|
|
|
343 }
|
|
|
344
|
|
|
345 if (!dct_NxM_costable_x)
|
|
|
346 dct_NxM_costable_x = alloc_coeffs(cols, cols);
|
|
|
347 if (!dct_NxM_costable_y)
|
|
|
348 dct_NxM_costable_y = alloc_coeffs(rows, rows);
|
|
|
349
|
|
|
350 N = cols;
|
|
|
351 M = rows;
|
|
|
352
|
|
|
353 for (i = 0; i < cols; i++) {
|
|
|
354 for (j = 0; j < cols; j++) {
|
|
|
355 dct_NxM_costable_x[i][j] = cx * cos((M_PI * ((2*i + 1) * j)) / (double) (2 * N));
|
|
|
356 }
|
|
|
357 }
|
|
|
358
|
|
|
359 for (i = 0; i < rows; i++) {
|
|
|
360 for (j = 0; j < rows; j++) {
|
|
|
361 dct_NxM_costable_y[i][j] = cy * cos((M_PI * ((2*i + 1) * j)) / (double) (2 * M));
|
|
|
362 }
|
|
|
363 }
|
|
|
364 }
|
|
|
365
|
|
|
366 void fdct_NxM(gray **pixels, double **dcts) {
|
|
|
367 int x, y;
|
|
|
368 int i, j;
|
|
|
369 double t;
|
|
|
370 double cx0 = sqrt(1.0 / N);
|
|
|
371 double cy0 = sqrt(1.0 / M);
|
|
|
372
|
|
|
373 t = 0.0;
|
|
|
374 for (x = 0; x < N; x++)
|
|
|
375 for (y = 0; y < M; y++)
|
|
|
376 t += ((int) pixels[y][x] - 128);
|
|
|
377 dcts[0][0] = cx0 * cy0 * t;
|
|
|
378
|
|
|
379 for (i = 1; i < N; i++) {
|
|
|
380 t = 0.0;
|
|
12
|
381 for (x = 0; x < N; x++) {
|
|
|
382 double s = 0.0;
|
|
|
383 for (y = 0; y < M; y++) {
|
|
|
384 s += ((int) pixels[y][x] - 128);
|
|
|
385 }
|
|
|
386 t += s * dct_NxM_costable_x[x][i];
|
|
|
387 }
|
|
0
|
388 dcts[0][i] = cy0 * t;
|
|
|
389 }
|
|
|
390
|
|
|
391 for (j = 1; j < M; j++) {
|
|
|
392 t = 0.0;
|
|
12
|
393 for (y = 0; y < M; y++) {
|
|
|
394 double s = 0.0;
|
|
|
395 for (x = 0; x < N; x++) {
|
|
|
396 s += ((int) pixels[y][x] - 128);
|
|
|
397 }
|
|
|
398 t += s * dct_NxM_costable_y[y][j];
|
|
|
399 }
|
|
0
|
400 dcts[j][0] = cx0 * t;
|
|
|
401 }
|
|
|
402
|
|
12
|
403 for (i = 1; i < N; i++) {
|
|
|
404 for (j = 1; j < M; j++) {
|
|
0
|
405 t = 0.0;
|
|
12
|
406 for (x = 0; x < N; x++) {
|
|
|
407 double s = 0;
|
|
|
408 for (y = 0; y < M; y++) {
|
|
|
409 s += ((int) pixels[y][x] - 128) * dct_NxM_costable_y[y][j];
|
|
|
410 }
|
|
|
411 t += s * dct_NxM_costable_x[x][i];
|
|
|
412 }
|
|
0
|
413 dcts[j][i] = t;
|
|
|
414 }
|
|
12
|
415 }
|
|
0
|
416 }
|
|
|
417
|
|
|
418 void idct_NxM(double **dcts, gray **pixels) {
|
|
|
419 int x, y;
|
|
|
420 int i, j;
|
|
|
421 double cx0 = sqrt(1.0 / N);
|
|
|
422 double cy0 = sqrt(1.0 / M);
|
|
|
423 double t;
|
|
|
424
|
|
|
425 for (x = 0; x < N; x++) {
|
|
|
426 for (y = 0; y < M; y++) {
|
|
|
427
|
|
|
428 t = cx0 * cy0 * dcts[0][0];
|
|
|
429
|
|
|
430 for (i = 1; i < N; i++)
|
|
|
431 t += cy0 * dcts[0][i] * dct_NxM_costable_x[x][i];
|
|
|
432
|
|
|
433 for (j = 1; j < M; j++)
|
|
|
434 t += cx0 * dcts[j][0] * dct_NxM_costable_y[y][j];
|
|
|
435
|
|
12
|
436 for (i = 1; i < N; i++) {
|
|
|
437 double s = 0.0;
|
|
|
438 for (j = 1; j < M; j++) {
|
|
|
439 s += dcts[j][i] * dct_NxM_costable_y[y][j];
|
|
|
440 }
|
|
|
441 t += s * dct_NxM_costable_x[x][i];
|
|
|
442 }
|
|
0
|
443
|
|
|
444 pixels[y][x] = PIXELRANGE((int) (t + 128.5));
|
|
|
445 }
|
|
|
446 }
|
|
|
447 }
|
|
|
448
|
|
|
449 double C[NJPEG][NJPEG];
|
|
|
450 double Ct[NJPEG][NJPEG];
|
|
|
451 int Quantum[NJPEG][NJPEG];
|
|
|
452
|
|
|
453 void init_quantum_8x8(int quality) {
|
|
|
454 int i;
|
|
|
455 int j;
|
|
|
456
|
|
|
457 for (i = 0; i < NJPEG; i++)
|
|
|
458 for ( j = 0 ; j < NJPEG ; j++ )
|
|
|
459 Quantum[ i ][ j ] = 1 + ( ( 1 + i + j ) * quality );
|
|
|
460 }
|
|
|
461
|
|
|
462 void init_quantum_JPEG_lumin(int quality) {
|
|
|
463 int i;
|
|
|
464 int j;
|
|
|
465
|
|
|
466 if (quality < 50)
|
|
|
467 quality = 5000 / quality;
|
|
|
468 else
|
|
|
469 quality = 200 - quality * 2;
|
|
|
470
|
|
|
471 for (i = 0; i < NJPEG; i++)
|
|
|
472 for (j = 0 ; j < NJPEG ; j++)
|
|
|
473 if (quality)
|
|
|
474 Quantum[i][j] = (JPEG_lumin_quant_table[i][j] * quality + 50) / 100;
|
|
|
475 else
|
|
|
476 Quantum[i][j] = JPEG_lumin_quant_table[i][j];
|
|
|
477 }
|
|
|
478
|
|
|
479 void init_quantum_JPEG_chromin(int quality) {
|
|
|
480 int i;
|
|
|
481 int j;
|
|
|
482
|
|
|
483 if (quality < 50)
|
|
|
484 quality = 5000 / quality;
|
|
|
485 else
|
|
|
486 quality = 200 - quality * 2;
|
|
|
487
|
|
|
488 for (i = 0; i < NJPEG; i++)
|
|
|
489 for (j = 0 ; j < NJPEG ; j++)
|
|
|
490 if (quality)
|
|
|
491 Quantum[i][j] = (JPEG_lumin_quant_table[i][j] * quality + 50) / 100;
|
|
|
492 else
|
|
|
493 Quantum[i][j] = JPEG_lumin_quant_table[i][j];
|
|
|
494 }
|
|
|
495
|
|
|
496 void quantize_8x8(double **transform) {
|
|
|
497 int i;
|
|
|
498 int j;
|
|
|
499
|
|
|
500 for (i = 0; i < NJPEG; i++)
|
|
|
501 for (j = 0; j < NJPEG; j++)
|
|
|
502 transform[i][j] = ROUND(transform[i][j] / Quantum[i][j]);
|
|
|
503 }
|
|
|
504
|
|
|
505 void dequantize_8x8(double **transform) {
|
|
|
506 int i;
|
|
|
507 int j;
|
|
|
508
|
|
|
509 for (i = 0; i < NJPEG; i++)
|
|
|
510 for (j = 0; j < NJPEG; j++)
|
|
|
511 transform[i][j] = ROUND(transform[i][j] * Quantum[i][j]);
|
|
|
512 }
|
|
|
513
|
|
|
514 void init_dct_8x8() {
|
|
|
515 int i;
|
|
|
516 int j;
|
|
|
517 double pi = atan( 1.0 ) * 4.0;
|
|
|
518
|
|
|
519 for ( j = 0 ; j < NJPEG ; j++ ) {
|
|
|
520 C[ 0 ][ j ] = 1.0 / sqrt( (double) NJPEG );
|
|
|
521 Ct[ j ][ 0 ] = C[ 0 ][ j ];
|
|
|
522 }
|
|
|
523
|
|
|
524 for ( i = 1 ; i < NJPEG ; i++ )
|
|
|
525 for ( j = 0 ; j < NJPEG ; j++ ) {
|
|
|
526 C[ i ][ j ] = sqrt( 2.0 / NJPEG ) * cos( pi * ( 2 * j + 1 ) * i / ( 2.0 * NJPEG ) );
|
|
|
527 Ct[ j ][ i ] = C[ i ][ j ];
|
|
|
528 }
|
|
|
529 }
|
|
|
530
|
|
|
531 /*
|
|
|
532 * The Forward DCT routine implements the matrix function:
|
|
|
533 *
|
|
|
534 * DCT = C * pixels * Ct
|
|
|
535 */
|
|
|
536
|
|
|
537 void fdct_8x8(gray **input, double **output) {
|
|
|
538 double temp[NJPEG][NJPEG];
|
|
|
539 double temp1;
|
|
|
540 int i;
|
|
|
541 int j;
|
|
|
542 int k;
|
|
|
543
|
|
|
544 /* MatrixMultiply( temp, input, Ct ); */
|
|
|
545 for ( i = 0 ; i < NJPEG ; i++ ) {
|
|
|
546 for ( j = 0 ; j < NJPEG ; j++ ) {
|
|
|
547 temp[ i ][ j ] = 0.0;
|
|
|
548 for ( k = 0 ; k < NJPEG ; k++ )
|
|
3
|
549 temp[ i ][ j ] += ( (int) input[ i ][ k ] - 128 ) *
|
|
0
|
550 Ct[ k ][ j ];
|
|
|
551 }
|
|
|
552 }
|
|
|
553
|
|
|
554 /* MatrixMultiply( output, C, temp ); */
|
|
|
555 for ( i = 0 ; i < NJPEG ; i++ ) {
|
|
|
556 for ( j = 0 ; j < NJPEG ; j++ ) {
|
|
|
557 temp1 = 0.0;
|
|
|
558 for ( k = 0 ; k < NJPEG ; k++ )
|
|
|
559 temp1 += C[ i ][ k ] * temp[ k ][ j ];
|
|
|
560 output[ i ][ j ] = temp1;
|
|
|
561 }
|
|
|
562 }
|
|
|
563 }
|
|
|
564
|
|
|
565 void fdct_block_8x8(gray **input, int col, int row, double **output) {
|
|
8
|
566 int i;
|
|
0
|
567 gray *input_array[NJPEG];
|
|
|
568
|
|
|
569 for (i = 0; i < NJPEG; i++)
|
|
|
570 input_array[i] = &input[row + i][col];
|
|
|
571
|
|
|
572 fdct_8x8(input_array, output);
|
|
|
573 }
|
|
|
574
|
|
|
575 /*
|
|
|
576 * The Inverse DCT routine implements the matrix function:
|
|
|
577 *
|
|
|
578 * pixels = C * DCT * Ct
|
|
|
579 */
|
|
|
580
|
|
|
581 void idct_8x8(double **input, gray **output) {
|
|
|
582 double temp[ NJPEG ][ NJPEG ];
|
|
|
583 double temp1;
|
|
|
584 int i;
|
|
|
585 int j;
|
|
|
586 int k;
|
|
|
587
|
|
|
588 /* MatrixMultiply( temp, input, C ); */
|
|
|
589 for ( i = 0 ; i < NJPEG ; i++ ) {
|
|
|
590 for ( j = 0 ; j < NJPEG ; j++ ) {
|
|
|
591 temp[ i ][ j ] = 0.0;
|
|
|
592 for ( k = 0 ; k < NJPEG ; k++ )
|
|
|
593 temp[ i ][ j ] += input[ i ][ k ] * C[ k ][ j ];
|
|
|
594 }
|
|
|
595 }
|
|
|
596
|
|
|
597 /* MatrixMultiply( output, Ct, temp ); */
|
|
|
598 for ( i = 0 ; i < NJPEG ; i++ ) {
|
|
|
599 for ( j = 0 ; j < NJPEG ; j++ ) {
|
|
|
600 temp1 = 0.0;
|
|
|
601 for ( k = 0 ; k < NJPEG ; k++ )
|
|
|
602 temp1 += Ct[ i ][ k ] * temp[ k ][ j ];
|
|
3
|
603 temp1 += 128.0;
|
|
0
|
604 output[i][j] = PIXELRANGE(ROUND(temp1));
|
|
|
605 }
|
|
|
606 }
|
|
|
607 }
|
|
|
608
|
|
|
609 void idct_block_8x8(double **input, gray **output, int col, int row) {
|
|
8
|
610 int i;
|
|
0
|
611 gray *output_array[NJPEG];
|
|
|
612
|
|
|
613 for (i = 0; i < NJPEG; i++)
|
|
|
614 output_array[i] = &output[row + i][col];
|
|
|
615
|
|
|
616 idct_8x8(input, output_array);
|
|
|
617 }
|
|
|
618
|
|
|
619 int is_middle_frequency_coeff_8x8(int coeff) {
|
|
|
620 switch (coeff) {
|
|
|
621 case 3:
|
|
|
622 case 10:
|
|
|
623 case 17:
|
|
|
624 case 24:
|
|
|
625 return 1;
|
|
|
626 case 4:
|
|
|
627 case 11:
|
|
|
628 case 18:
|
|
|
629 case 25:
|
|
|
630 case 32:
|
|
|
631 return 2;
|
|
|
632 case 5:
|
|
|
633 case 12:
|
|
|
634 case 19:
|
|
|
635 case 26:
|
|
|
636 case 33:
|
|
|
637 case 40:
|
|
|
638 return 3;
|
|
|
639 case 13:
|
|
|
640 case 20:
|
|
|
641 case 27:
|
|
|
642 case 34:
|
|
|
643 case 41:
|
|
|
644 return 4;
|
|
|
645 case 28:
|
|
|
646 case 35:
|
|
|
647 return 5;
|
|
|
648 default:
|
|
|
649 return 0;
|
|
|
650 }
|
|
|
651 }
|
