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