home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Interactive Guide / c-cplusplus-interactive-guide.iso / c_ref / csource1 / mpegplay / jrevdct.c < prev    next >
C/C++ Source or Header  |  1993-02-02  |  44KB  |  1,462 lines
  1. /*
  2.  * jrevdct.c
  3.  *
  4.  * Copyright (C) 1991, 1992, Thomas G. Lane.
  5.  * This file is part of the Independent JPEG Group's software.
  6.  * For conditions of distribution and use, see the accompanying README file.
  7.  *
  8.  * This file contains the basic inverse-DCT transformation subroutine.
  9.  *
  10.  * This implementation is based on an algorithm described in
  11.  *   C. Loeffler, A. Ligtenberg and G. Moschytz, "Practical Fast 1-D DCT
  12.  *   Algorithms with 11 Multiplications", Proc. Int'l. Conf. on Acoustics,
  13.  *   Speech, and Signal Processing 1989 (ICASSP '89), pp. 988-991.
  14.  * The primary algorithm described there uses 11 multiplies and 29 adds.
  15.  * We use their alternate method with 12 multiplies and 32 adds.
  16.  * The advantage of this method is that no data path contains more than one
  17.  * multiplication; this allows a very simple and accurate implementation in
  18.  * scaled fixed-point arithmetic, with a minimal number of shifts.
  19.  * 
  20.  * I've made lots of modifications to attempt to take advantage of the
  21.  * sparse nature of the DCT matrices we're getting.  Although the logic
  22.  * is cumbersome, it's straightforward and the resulting code is much
  23.  * faster.
  24.  *
  25.  * A better way to do this would be to pass in the DCT block as a sparse
  26.  * matrix, perhaps with the difference cases encoded.
  27.  */
  28.  
  29. #include <string.h>
  30. #include "video.h"
  31. #include "proto.h"
  32.  
  33. #define GLOBAL            /* a function referenced thru EXTERNs */
  34.   
  35. /* We assume that right shift corresponds to signed division by 2 with
  36.  * rounding towards minus infinity.  This is correct for typical "arithmetic
  37.  * shift" instructions that shift in copies of the sign bit.  But some
  38.  * C compilers implement >> with an unsigned shift.  For these machines you
  39.  * must define RIGHT_SHIFT_IS_UNSIGNED.
  40.  * RIGHT_SHIFT provides a proper signed right shift of an INT32 quantity.
  41.  * It is only applied with constant shift counts.  SHIFT_TEMPS must be
  42.  * included in the variables of any routine using RIGHT_SHIFT.
  43.  */
  44.   
  45. #ifdef RIGHT_SHIFT_IS_UNSIGNED
  46. #define SHIFT_TEMPS    INT32 shift_temp;
  47. #define RIGHT_SHIFT(x,shft)  \
  48.     ((shift_temp = (x)) < 0 ? \
  49.      (shift_temp >> (shft)) | ((~((INT32) 0)) << (32-(shft))) : \
  50.      (shift_temp >> (shft)))
  51. #else
  52. #define SHIFT_TEMPS
  53. #define RIGHT_SHIFT(x,shft)    ((x) >> (shft))
  54. #endif
  55.  
  56. /*
  57.  * This routine is specialized to the case DCTSIZE = 8.
  58.  */
  59.  
  60. #if DCTSIZE != 8
  61.   Sorry, this code only copes with 8x8 DCTs. /* deliberate syntax err */
  62. #endif
  63.  
  64.  
  65. /*
  66.  * A 2-D IDCT can be done by 1-D IDCT on each row followed by 1-D IDCT
  67.  * on each column.  Direct algorithms are also available, but they are
  68.  * much more complex and seem not to be any faster when reduced to code.
  69.  *
  70.  * The poop on this scaling stuff is as follows:
  71.  *
  72.  * Each 1-D IDCT step produces outputs which are a factor of sqrt(N)
  73.  * larger than the true IDCT outputs.  The final outputs are therefore
  74.  * a factor of N larger than desired; since N=8 this can be cured by
  75.  * a simple right shift at the end of the algorithm.  The advantage of
  76.  * this arrangement is that we save two multiplications per 1-D IDCT,
  77.  * because the y0 and y4 inputs need not be divided by sqrt(N).
  78.  *
  79.  * We have to do addition and subtraction of the integer inputs, which
  80.  * is no problem, and multiplication by fractional constants, which is
  81.  * a problem to do in integer arithmetic.  We multiply all the constants
  82.  * by CONST_SCALE and convert them to integer constants (thus retaining
  83.  * CONST_BITS bits of precision in the constants).  After doing a
  84.  * multiplication we have to divide the product by CONST_SCALE, with proper
  85.  * rounding, to produce the correct output.  This division can be done
  86.  * cheaply as a right shift of CONST_BITS bits.  We postpone shifting
  87.  * as long as possible so that partial sums can be added together with
  88.  * full fractional precision.
  89.  *
  90.  * The outputs of the first pass are scaled up by PASS1_BITS bits so that
  91.  * they are represented to better-than-integral precision.  These outputs
  92.  * require BITS_IN_JSAMPLE + PASS1_BITS + 3 bits; this fits in a 16-bit word
  93.  * with the recommended scaling.  (To scale up 12-bit sample data further, an
  94.  * intermediate INT32 array would be needed.)
  95.  *
  96.  * To avoid overflow of the 32-bit intermediate results in pass 2, we must
  97.  * have BITS_IN_JSAMPLE + CONST_BITS + PASS1_BITS <= 26.  Error analysis
  98.  * shows that the values given below are the most effective.
  99.  */
  100.  
  101. #ifdef EIGHT_BIT_SAMPLES
  102. #define PASS1_BITS  2
  103. #else
  104. #define PASS1_BITS  1        /* lose a little precision to avoid overflow */
  105. #endif
  106.  
  107. #define ONE    ((INT32) 1)
  108.  
  109. #define CONST_SCALE (ONE << CONST_BITS)
  110.  
  111. /* Convert a positive real constant to an integer scaled by CONST_SCALE.
  112.  * IMPORTANT: if your compiler doesn't do this arithmetic at compile time,
  113.  * you will pay a significant penalty in run time.  In that case, figure
  114.  * the correct integer constant values and insert them by hand.
  115.  */
  116.  
  117. #define FIX(x)    ((INT32) ((x) * CONST_SCALE + 0.5))
  118.  
  119. /* Descale and correctly round an INT32 value that's scaled by N bits.
  120.  * We assume RIGHT_SHIFT rounds towards minus infinity, so adding
  121.  * the fudge factor is correct for either sign of X.
  122.  */
  123.  
  124. #define DESCALE(x,n)  RIGHT_SHIFT((x) + (ONE << ((n)-1)), n)
  125.  
  126. /* Multiply an INT32 variable by an INT32 constant to yield an INT32 result.
  127.  * For 8-bit samples with the recommended scaling, all the variable
  128.  * and constant values involved are no more than 16 bits wide, so a
  129.  * 16x16->32 bit multiply can be used instead of a full 32x32 multiply;
  130.  * this provides a useful speedup on many machines.
  131.  * There is no way to specify a 16x16->32 multiply in portable C, but
  132.  * some C compilers will do the right thing if you provide the correct
  133.  * combination of casts.
  134.  * NB: for 12-bit samples, a full 32-bit multiplication will be needed.
  135.  */
  136.  
  137. #ifdef EIGHT_BIT_SAMPLES
  138. #ifdef SHORTxSHORT_32        /* may work if 'int' is 32 bits */
  139. #define MULTIPLY(var,const)  (((INT16) (var)) * ((INT16) (const)))
  140. #endif
  141. #ifdef SHORTxLCONST_32        /* known to work with Microsoft C 6.0 */
  142. #define MULTIPLY(var,const)  (((INT16) (var)) * ((INT32) (const)))
  143. #endif
  144. #endif
  145.  
  146. #ifndef MULTIPLY        /* default definition */
  147. #define MULTIPLY(var,const)  ((var) * (const))
  148. #endif
  149.  
  150. /* Precomputed idct value arrays. */
  151.  
  152. static DCTELEM PreIDCT[64][64];
  153.  
  154. /* Pre compute singleton coefficient IDCT values. */
  155. void
  156. init_pre_idct() {
  157.   int i;
  158.   void j_rev_dct();
  159.  
  160.   for (i=0; i<64; i++) {
  161.     memset((char *) PreIDCT[i], 0, 64*sizeof(DCTELEM));
  162.     PreIDCT[i][i] = 2048;
  163.     j_rev_dct(PreIDCT[i]);
  164.   }
  165. }
  166.  
  167. #ifndef ORIG_DCT
  168.   
  169.  
  170. /*
  171.  * Perform the inverse DCT on one block of coefficients.
  172.  */
  173.  
  174. void
  175. j_rev_dct_sparse (data, pos)
  176.      DCTBLOCK data;
  177.      int pos;
  178. {
  179.   register DCTELEM *dataptr;
  180.   short int val;
  181.   DCTELEM *ndataptr;
  182.   int scale, coeff, rr;
  183.   register int *dp;
  184.   register int v;
  185.  
  186.   /* If DC Coefficient. */
  187.   
  188.   if (pos == 0) {
  189.     dp = (int *)data;
  190.     v = *data;
  191.     /* Compute 32 bit value to assign.  This speeds things up a bit */
  192.     if (v < 0) val = (v-3)>>3;
  193.     else val = (v+4)>>3;
  194.     v = val | (val << 16);
  195.     dp[0] = v;      dp[1] = v;      dp[2] = v;      dp[3] = v;
  196.     dp[4] = v;      dp[5] = v;      dp[6] = v;      dp[7] = v;
  197.     dp[8] = v;      dp[9] = v;      dp[10] = v;      dp[11] = v;
  198.     dp[12] = v;      dp[13] = v;      dp[14] = v;      dp[15] = v;
  199.     dp[16] = v;      dp[17] = v;      dp[18] = v;      dp[19] = v;
  200.     dp[20] = v;      dp[21] = v;      dp[22] = v;      dp[23] = v;
  201.     dp[24] = v;      dp[25] = v;      dp[26] = v;      dp[27] = v;
  202.     dp[28] = v;      dp[29] = v;      dp[30] = v;      dp[31] = v;
  203.     return;
  204.   }
  205.   
  206.   /* Some other coefficient. */
  207.   dataptr = (DCTELEM *)data;
  208.   coeff = dataptr[pos];
  209.   ndataptr = PreIDCT[pos];
  210.  
  211.   for (rr=0; rr<4; rr++) {
  212.     dataptr[0] = (ndataptr[0] * coeff) >> (CONST_BITS-2);
  213.     dataptr[1] = (ndataptr[1] * coeff) >> (CONST_BITS-2);
  214.     dataptr[2] = (ndataptr[2] * coeff) >> (CONST_BITS-2);
  215.     dataptr[3] = (ndataptr[3] * coeff) >> (CONST_BITS-2);
  216.     dataptr[4] = (ndataptr[4] * coeff) >> (CONST_BITS-2);
  217.     dataptr[5] = (ndataptr[5] * coeff) >> (CONST_BITS-2);
  218.     dataptr[6] = (ndataptr[6] * coeff) >> (CONST_BITS-2);
  219.     dataptr[7] = (ndataptr[7] * coeff) >> (CONST_BITS-2);
  220.     dataptr[8] = (ndataptr[8] * coeff) >> (CONST_BITS-2);
  221.     dataptr[9] = (ndataptr[