annotate Meerwald/dct.c @ 18:3bdb67e76858

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

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