Mercurial > hg > wm

comparison Meerwald/dct.c @ 0:be303a3f5ea8

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

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