1 /*************************************************************
2 Copyright (C) 1990, 1991, 1993 Andy C. Hung, all rights reserved.
3 PUBLIC DOMAIN LICENSE: Stanford University Portable Video Research
4 Group. If you use this software, you agree to the following: This
5 program package is purely experimental, and is licensed "as is".
6 Permission is granted to use, modify, and distribute this program
7 without charge for any purpose, provided this license/ disclaimer
8 notice appears in the copies. No warranty or maintenance is given,
9 either expressed or implied. In no event shall the author(s) be
10 liable to you or a third party for any special, incidental,
11 consequential, or other damages, arising out of the use or inability
12 to use the program for any purpose (or the loss of data), even if we
13 have been advised of such possibilities. Any public reference or
14 advertisement of this source code should refer to it as the Portable
15 Video Research Group (PVRG) code, and not by any author(s) (or
16 Stanford University) name.
17 *************************************************************/
18
19 /*
20 ************************************************************
21 transform.c
22
23 This file contains the reference DCT, the zig-zag and quantization
24 algorithms.
25
26 ************************************************************
27 */
28
29 /*LABEL transform.c */
30
31 /* Include files */
32
33 #include "globals.h"
34 #include "dct.h"
35 #include <math.h>
36 #include <stdlib.h> /* exit */
37
38 /*PUBLIC*/
39
40 extern void ReferenceDct();
41 extern void ReferenceIDct();
42 extern void TransposeMatrix();
43 extern void Quantize();
44 extern void IQuantize();
45 extern void PreshiftDctMatrix();
46 extern void PostshiftIDctMatrix();
47 extern void BoundDctMatrix();
48 extern void BoundIDctMatrix();
49 extern void ZigzagMatrix();
50 extern void IZigzagMatrix();
51 extern int *ScaleMatrix();
52 extern void PrintMatrix();
53 extern void ClearMatrix();
54
55 static void DoubleReferenceDct1D();
56 static void DoubleReferenceIDct1D();
57 static void DoubleTransposeMatrix();
58
59 /*PRIVATE*/
60
61 /* The transposition indices */
62
63 int transpose_index[] = /* Is a transpose map for matrix transp. */
64 {0, 8, 16, 24, 32, 40, 48, 56,
65 1, 9, 17, 25, 33, 41, 49, 57,
66 2, 10, 18, 26, 34, 42, 50, 58,
67 3, 11, 19, 27, 35, 43, 51, 59,
68 4, 12, 20, 28, 36, 44, 52, 60,
69 5, 13, 21, 29, 37, 45, 53, 61,
70 6, 14, 22, 30, 38, 46, 54, 62,
71 7, 15, 23, 31, 39, 47, 55, 63};
72
73 int zigzag_index[] = /* Is zig-zag map for matrix -> scan array */
74 {0, 1, 5, 6, 14, 15, 27, 28,
75 2, 4, 7, 13, 16, 26, 29, 42,
76 3, 8, 12, 17, 25, 30, 41, 43,
77 9, 11, 18, 24, 31, 40, 44, 53,
78 10, 19, 23, 32, 39, 45, 52, 54,
79 20, 22, 33, 38, 46, 51, 55, 60,
80 21, 34, 37, 47, 50, 56, 59, 61,
81 35, 36, 48, 49, 57, 58, 62, 63};
82
83 int izigzag_index[] =
84 {0, 1, 8, 16, 9, 2, 3, 10,
85 17, 24, 32, 25, 18, 11, 4, 5,
86 12, 19, 26, 33, 40, 48, 41, 34,
87 27, 20, 13, 6, 7, 14, 21, 28,
88 35, 42, 49, 56, 57, 50, 43, 36,
89 29, 22, 15, 23, 30, 37, 44, 51,
90 58, 59, 52, 45, 38, 31, 39, 46,
91 53, 60, 61, 54, 47, 55, 62, 63};
92
93 /*Some definitions */
94
95 #define MakeMatrix() (int *) calloc(BLOCKSIZE,sizeof(int))
96 #define FixedMultiply(s,x,y) x = ((x * y) >> s);
97 #define DCT_OFFSET 128
98
99 /*START*/
100
101 /*BFUNC
102
103 ReferenceDct() does a reference DCT on the input (matrix) and output
104 (new matrix).
105
106 EFUNC*/
107
ReferenceDct(matrix,newmatrix)108 void ReferenceDct(matrix,newmatrix)
109 int *matrix;
110 int *newmatrix;
111 {
112 BEGIN("ReferenceDct")
113 int *mptr;
114 double *sptr,*dptr;
115 double sourcematrix[BLOCKSIZE],destmatrix[BLOCKSIZE];
116
117 for(sptr=sourcematrix,mptr=matrix;mptr<matrix+BLOCKSIZE;mptr++)
118 { /* Convert to doubles */
119 *(sptr++) = (double) *mptr;
120 }
121 for(dptr = destmatrix,sptr=sourcematrix;
122 sptr<sourcematrix+BLOCKSIZE;sptr+=BLOCKWIDTH)
123 { /* Do DCT on rows */
124 DoubleReferenceDct1D(sptr,dptr);
125 dptr+=BLOCKWIDTH;
126 }
127 DoubleTransposeMatrix(destmatrix,sourcematrix); /* Transpose */
128 for(dptr = destmatrix,sptr=sourcematrix;
129 sptr<sourcematrix+BLOCKSIZE;sptr+=BLOCKWIDTH)
130 { /* Do DCT on columns */
131 DoubleReferenceDct1D(sptr,dptr);
132 dptr+=BLOCKWIDTH;
133 }
134 DoubleTransposeMatrix(destmatrix,sourcematrix); /* Transpose */
135 for(sptr = sourcematrix,mptr=newmatrix;
136 mptr<newmatrix+BLOCKSIZE;sptr++)
137 { /* NB: Inversion on counter */
138 *(mptr++) = (int) (*sptr > 0 ? (*(sptr)+0.5):(*(sptr)-0.5));
139 }
140 }
141
142 /*BFUNC
143
144 DoubleReferenceDCT1D() does a 8 point dct on an array of double
145 input and places the result in a double output.
146
147 EFUNC*/
148
DoubleReferenceDct1D(ivect,ovect)149 static void DoubleReferenceDct1D(ivect,ovect)
150 double *ivect;
151 double *ovect;
152 {
153 BEGIN("DoubleReferenceDct1D")
154 double *mptr,*iptr,*optr;
155
156 for(mptr=DctMatrix,optr=ovect;optr<ovect+BLOCKWIDTH;optr++)
157 { /* 1d dct is just matrix multiply */
158 for(*optr=0,iptr=ivect;iptr<ivect+BLOCKWIDTH;iptr++)
159 {
160 *optr += *iptr*(*(mptr++));
161 }
162 }
163 }
164
165 /*BFUNC
166
167 ReferenceIDct() is used to perform a reference 8x8 inverse dct. It is
168 a balanced IDCT. It takes the input (matrix) and puts it into the
169 output (newmatrix).
170
171 EFUNC*/
172
ReferenceIDct(matrix,newmatrix)173 void ReferenceIDct(matrix,newmatrix)
174 int *matrix;
175 int *newmatrix;
176 {
177 BEGIN("ReferenceIDct")
178 int *mptr;
179 double *sptr,*dptr;
180 double sourcematrix[BLOCKSIZE],destmatrix[BLOCKSIZE];
181
182 for(sptr = sourcematrix,mptr=matrix;mptr<matrix+BLOCKSIZE;mptr++)
183 {
184 *(sptr++) = (double) *mptr;
185 }
186 for(dptr = destmatrix,sptr=sourcematrix;
187 sptr<sourcematrix+BLOCKSIZE;sptr+=BLOCKWIDTH)
188 {
189 DoubleReferenceIDct1D(sptr,dptr);
190 dptr+=BLOCKWIDTH;
191 }
192 DoubleTransposeMatrix(destmatrix,sourcematrix);
193 for(dptr = destmatrix,sptr=sourcematrix;
194 sptr<sourcematrix+BLOCKSIZE;sptr+=BLOCKWIDTH)
195 {
196 DoubleReferenceIDct1D(sptr,dptr);
197 dptr+=BLOCKWIDTH;
198 }
199 DoubleTransposeMatrix(destmatrix,sourcematrix);
200 for(sptr = sourcematrix,mptr=newmatrix;mptr<newmatrix+BLOCKSIZE;sptr++)
201 { /* NB: Inversion on counter */
202 *(mptr++) = (int) (*sptr > 0 ? (*(sptr)+0.5):(*(sptr)-0.5));
203 }
204 }
205
206 /*BFUNC
207
208 DoubleReferenceIDct1D() does an 8 point inverse dct on ivect and
209 puts the output in ovect.
210
211 EFUNC*/
212
DoubleReferenceIDct1D(ivect,ovect)213 static void DoubleReferenceIDct1D(ivect,ovect)
214 double *ivect;
215 double *ovect;
216 {
217 BEGIN("DoubleReferenceIDct1D")
218 double *mptr,*iptr,*optr;
219
220 for(mptr = IDctMatrix,optr=ovect;optr<ovect+BLOCKWIDTH;optr++)
221 {
222 for(*optr=0,iptr=ivect;iptr<ivect+BLOCKWIDTH;iptr++)
223 {
224 *optr += *iptr*(*(mptr++));
225 }
226 }
227 }
228
229 /*BFUNC
230
231 TransposeMatrix transposes an input matrix and puts the output in
232 newmatrix.
233
234 EFUNC*/
235
TransposeMatrix(matrix,newmatrix)236 void TransposeMatrix(matrix,newmatrix)
237 int *matrix;
238 int *newmatrix;
239 {
240 BEGIN("TransposeMatrix")
241 int *tptr;
242
243 for(tptr=transpose_index;tptr<transpose_index+BLOCKSIZE;tptr++)
244 {
245 *(newmatrix++) = matrix[*tptr];
246 }
247 }
248
249 /*BFUNC
250
251 DoubleTransposeMatrix transposes a double input matrix and puts the
252 double output in newmatrix.
253
254 EFUNC*/
255
DoubleTransposeMatrix(matrix,newmatrix)256 static void DoubleTransposeMatrix(matrix,newmatrix)
257 double *matrix;
258 double *newmatrix;
259 {
260 BEGIN("DoubleTransposeMatrix")
261 int *tptr;
262
263 for(tptr=transpose_index;tptr<transpose_index+BLOCKSIZE;tptr++)
264 {
265 *(newmatrix++) = matrix[*tptr];
266 }
267 }
268
269 /*BFUNC
270
271 Quantize() quantizes an input matrix and puts the output in qmatrix.
272
273 EFUNC*/
274
Quantize(matrix,qmatrix)275 void Quantize(matrix,qmatrix)
276 int *matrix;
277 int *qmatrix;
278 {
279 BEGIN("Quantize")
280 int *mptr;
281
282 if (!qmatrix)
283 {
284 WHEREAMI();
285 printf("No quantization matrix specified!\n");
286 exit(ERROR_BOUNDS);
287 }
288 for(mptr=matrix;mptr<matrix+BLOCKSIZE;mptr++)
289 {
290 if (*mptr > 0) /* Rounding is different for +/- coeffs */
291 {
292 *mptr = (*mptr + *qmatrix/2)/ (*qmatrix);
293 qmatrix++;
294 }
295 else
296 {
297 *mptr = (*mptr - *qmatrix/2)/ (*qmatrix);
298 qmatrix++;
299 }
300 }
301 }
302
303 /*BFUNC
304
305 IQuantize() takes an input matrix and does an inverse quantization
306 and puts the output int qmatrix.
307
308 EFUNC*/
309
IQuantize(matrix,qmatrix)310 void IQuantize(matrix,qmatrix)
311 int *matrix;
312 int *qmatrix;
313 {
314 BEGIN("IQuantize")
315 int *mptr;
316
317 if (!qmatrix)
318 {
319 WHEREAMI();
320 printf("No quantization matrix specified!\n");
321 exit(ERROR_BOUNDS);
322 }
323 for(mptr=matrix;mptr<matrix+BLOCKSIZE;mptr++)
324 {
325 *mptr = *mptr*(*qmatrix);
326 qmatrix++;
327 }
328 }
329
330 /*BFUNC
331
332 PreshiftDctMatrix() subtracts 128 (2048) from all 64 elements of an 8x8
333 matrix. This results in a balanced DCT without any DC bias.
334
335 EFUNC*/
336
337
PreshiftDctMatrix(matrix,shift)338 void PreshiftDctMatrix(matrix,shift)
339 int *matrix;
340 int shift;
341 {
342 BEGIN("PreshiftDctMatrix")
343 int *mptr;
344
345 for(mptr=matrix;mptr<matrix+BLOCKSIZE;mptr++) {*mptr -= shift;}
346 }
347
348 /*BFUNC
349
350 PostshiftIDctMatrix() adds 128 (2048) to all 64 elements of an 8x8 matrix.
351 This results in strictly positive values for all pixel coefficients.
352
353 EFUNC*/
354
PostshiftIDctMatrix(matrix,shift)355 void PostshiftIDctMatrix(matrix,shift)
356 int *matrix;
357 int shift;
358 {
359 BEGIN("PostshiftIDctMatrix")
360 int *mptr;
361
362 for(mptr=matrix;mptr<matrix+BLOCKSIZE;mptr++) {*mptr += shift;}
363 }
364
365 /*BFUNC
366
367 BoundDctMatrix() clips the Dct matrix such that it is no larger than
368 a 10 (1023) bit word or 14 bit word (4095).
369
370 EFUNC*/
371
BoundDctMatrix(matrix,Bound)372 void BoundDctMatrix(matrix,Bound)
373 int *matrix;
374 int Bound;
375 {
376 BEGIN("BoundDctMatrix")
377 int *mptr;
378
379 for(mptr=matrix;mptr<matrix+BLOCKSIZE;mptr++)
380 {
381 if (*mptr+Bound < 0)
382 *mptr = -Bound;
383 else if (*mptr-Bound > 0)
384 *mptr = Bound;
385 }
386 }
387
388 /*BFUNC
389
390 BoundIDctMatrix bounds the inverse dct matrix so that no pixel has a
391 value greater than 255 (4095) or less than 0.
392
393 EFUNC*/
394
BoundIDctMatrix(matrix,Bound)395 void BoundIDctMatrix(matrix,Bound)
396 int *matrix;
397 int Bound;
398 {
399 BEGIN("BoundIDctMatrix")
400 int *mptr;
401
402 for(mptr=matrix;mptr<matrix+BLOCKSIZE;mptr++)
403 {
404 if (*mptr < 0) {*mptr = 0;}
405 else if (*mptr > Bound) {*mptr = Bound;}
406 }
407 }
408
409 /*BFUNC
410
411 IZigzagMatrix() performs an inverse zig-zag translation on the
412 input imatrix and places the output in omatrix.
413
414 EFUNC*/
415
IZigzagMatrix(imatrix,omatrix)416 void IZigzagMatrix(imatrix,omatrix)
417 int *imatrix;
418 int *omatrix;
419 {
420 BEGIN("IZigzagMatrix")
421 int *tptr;
422
423 for(tptr=zigzag_index;tptr<zigzag_index+BLOCKSIZE;tptr++)
424 {
425 *(omatrix++) = imatrix[*tptr];
426 }
427 }
428
429
430 /*BFUNC
431
432 ZigzagMatrix() performs a zig-zag translation on the input imatrix
433 and puts the output in omatrix.
434
435 EFUNC*/
436
ZigzagMatrix(imatrix,omatrix)437 void ZigzagMatrix(imatrix,omatrix)
438 int *imatrix;
439 int *omatrix;
440 {
441 BEGIN("ZigzagMatrix")
442 int *tptr;
443
444 for(tptr=zigzag_index;tptr<zigzag_index+BLOCKSIZE;tptr++)
445 {
446 omatrix[*tptr] = *(imatrix++);
447 }
448 }
449
450
451 /*BFUNC
452
453 ScaleMatrix() does a matrix scale appropriate to the old Q-factor. It
454 returns the matrix created.
455
456 EFUNC*/
457
ScaleMatrix(Numerator,Denominator,LongFlag,Matrix)458 int *ScaleMatrix(Numerator,Denominator,LongFlag,Matrix)
459 int Numerator;
460 int Denominator;
461 int LongFlag;
462 int *Matrix;
463 {
464 BEGIN("ScaleMatrix")
465 int *Temp,*tptr;
466 int Limit;
467
468 Limit=((LongFlag)?65535:255);
469 if (!(Temp = MakeMatrix()))
470 {
471 WHEREAMI();
472 printf("Cannot allocate space for new matrix.\n");
473 exit(ERROR_MEMORY);
474 }
475 for(tptr=Temp;tptr<Temp+BLOCKSIZE;tptr++)
476 {
477 *tptr = (*(Matrix++) * Numerator)/Denominator;
478 if (*tptr > Limit)
479 *tptr = Limit;
480 else if (*tptr < 1)
481 *tptr = 1;
482 }
483 return(Temp);
484 }
485
486 /*BFUNC
487
488 PrintMatrix() prints an 8x8 matrix in row/column form.
489
490 EFUNC*/
491
PrintMatrix(matrix)492 void PrintMatrix(matrix)
493 int *matrix;
494 {
495 BEGIN("PrintMatrix")
496 int i,j;
497
498 if (matrix)
499 {
500 for(i=0;i<BLOCKHEIGHT;i++)
501 {
502 for(j=0;j<BLOCKWIDTH;j++) {printf("%6d ",*(matrix++));}
503 printf("\n");
504 }
505 }
506 else {printf("Null\n");}
507 }
508
509
510 /*BFUNC
511
512 ClearMatrix() sets all the elements of a matrix to be zero.
513
514 EFUNC*/
515
ClearMatrix(matrix)516 void ClearMatrix(matrix)
517 int *matrix;
518 {
519 BEGIN("ClearMatrix")
520 int *mptr;
521
522 for(mptr=matrix;mptr<matrix+BLOCKSIZE;mptr++) {*mptr = 0;}
523 }
524
525 /*END*/
526