home *** CD-ROM | disk | FTP | other *** search
/ Photo CD Demo 1 / Demo.bin / compresn / jpegv3sr / jrevdct.c < prev    next >
C/C++ Source or Header  |  1992-02-07  |  8KB  |  222 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 Appendix A.2 of the book
  11.  * "Discrete Cosine Transform---Algorithms, Advantages, Applications"
  12.  * by K.R. Rao and P. Yip  (Academic Press, Inc, London, 1990).
  13.  * It uses scaled fixed-point arithmetic instead of floating point.
  14.  */
  15.  
  16. #include "jinclude.h"
  17.  
  18. /*
  19.  * This routine is specialized to the case DCTSIZE = 8.
  20.  */
  21.  
  22. #if DCTSIZE != 8
  23.   Sorry, this code only copes with 8x8 DCTs. /* deliberate syntax err */
  24. #endif
  25.  
  26.  
  27. /* The poop on this scaling stuff is as follows:
  28.  *
  29.  * We have to do addition and subtraction of the integer inputs, which
  30.  * is no problem, and multiplication by fractional constants, which is
  31.  * a problem to do in integer arithmetic.  We multiply all the constants
  32.  * by DCT_SCALE and convert them to integer constants (thus retaining
  33.  * LG2_DCT_SCALE bits of precision in the constants).  After doing a
  34.  * multiplication we have to divide the product by DCT_SCALE, with proper
  35.  * rounding, to produce the correct output.  The division can be implemented
  36.  * cheaply as a right shift of LG2_DCT_SCALE bits.  The DCT equations also
  37.  * specify an additional division by 2 on the final outputs; this can be
  38.  * folded into the right-shift by shifting one more bit (see UNFIXH).
  39.  *
  40.  * If you are planning to recode this in assembler, you might want to set
  41.  * LG2_DCT_SCALE to 15.  This loses a bit of precision, but then all the
  42.  * multiplications are between 16-bit quantities (given 8-bit JSAMPLEs!)
  43.  * so you could use a signed 16x16=>32 bit multiply instruction instead of
  44.  * full 32x32 multiply.  Unfortunately there's no way to describe such a
  45.  * multiply portably in C, so we've gone for the extra bit of accuracy here.
  46.  */
  47.  
  48. #ifdef EIGHT_BIT_SAMPLES
  49. #define LG2_DCT_SCALE 16
  50. #else
  51. #define LG2_DCT_SCALE 15    /* lose a little precision to avoid overflow */
  52. #endif
  53.  
  54. #define ONE    ((INT32) 1)
  55.  
  56. #define DCT_SCALE (ONE << LG2_DCT_SCALE)
  57.  
  58. /* In some places we shift the inputs left by a couple more bits, */
  59. /* so that they can be added to fractional results without too much */
  60. /* loss of precision. */
  61. #define LG2_OVERSCALE 2
  62. #define OVERSCALE  (ONE << LG2_OVERSCALE)
  63. #define OVERSHIFT(x)  ((x) <<= LG2_OVERSCALE)
  64.  
  65. /* Scale a fractional constant by DCT_SCALE */
  66. #define FIX(x)    ((INT32) ((x) * DCT_SCALE + 0.5))
  67.  
  68. /* Scale a fractional constant by DCT_SCALE/OVERSCALE */
  69. /* Such a constant can be multiplied with an overscaled input */
  70. /* to produce something that's scaled by DCT_SCALE */
  71. #define FIXO(x)  ((INT32) ((x) * DCT_SCALE / OVERSCALE + 0.5))
  72.  
  73. /* Descale and correctly round a value that's scaled by DCT_SCALE */
  74. #define UNFIX(x)   RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1)), LG2_DCT_SCALE)
  75.  
  76. /* Same with an additional division by 2, ie, correctly rounded UNFIX(x/2) */
  77. #define UNFIXH(x)  RIGHT_SHIFT((x) + (ONE << LG2_DCT_SCALE), LG2_DCT_SCALE+1)
  78.  
  79. /* Take a value scaled by DCT_SCALE and round to integer scaled by OVERSCALE */
  80. #define UNFIXO(x)  RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1-LG2_OVERSCALE)),\
  81.                    LG2_DCT_SCALE-LG2_OVERSCALE)
  82.  
  83. /* Here are the constants we need */
  84. /* SIN_i_j is sine of i*pi/j, scaled by DCT_SCALE */
  85. /* COS_i_j is cosine of i*pi/j, scaled by DCT_SCALE */
  86.  
  87. #define SIN_1_4 FIX(0.707106781)
  88. #define COS_1_4 SIN_1_4
  89.  
  90. #define SIN_1_8 FIX(0.382683432)
  91. #define COS_1_8 FIX(0.923879533)
  92. #define SIN_3_8 COS_1_8
  93. #define COS_3_8 SIN_1_8
  94.  
  95. #define SIN_1_16 FIX(0.195090322)
  96. #define COS_1_16 FIX(0.980785280)
  97. #define SIN_7_16 COS_1_16
  98. #define COS_7_16 SIN_1_16
  99.  
  100. #define SIN_3_16 FIX(0.555570233)
  101. #define COS_3_16 FIX(0.831469612)
  102. #define SIN_5_16 COS_3_16
  103. #define COS_5_16 SIN_3_16
  104.  
  105. /* OSIN_i_j is sine of i*pi/j, scaled by DCT_SCALE/OVERSCALE */
  106. /* OCOS_i_j is cosine of i*pi/j, scaled by DCT_SCALE/OVERSCALE */
  107.  
  108. #define OSIN_1_4 FIXO(0.707106781)
  109. #define OCOS_1_4 OSIN_1_4
  110.  
  111. #define OSIN_1_8 FIXO(0.382683432)
  112. #define OCOS_1_8 FIXO(0.923879533)
  113. #define OSIN_3_8 OCOS_1_8
  114. #define OCOS_3_8 OSIN_1_8
  115.  
  116. #define OSIN_1_16 FIXO(0.195090322)
  117. #define OCOS_1_16 FIXO(0.980785280)
  118. #define OSIN_7_16 OCOS_1_16
  119. #define OCOS_7_16 OSIN_1_16
  120.  
  121. #define OSIN_3_16 FIXO(0.555570233)
  122. #define OCOS_3_16 FIXO(0.831469612)
  123. #define OSIN_5_16 OCOS_3_16
  124. #define OCOS_5_16 OSIN_3_16
  125.  
  126.  
  127. /*
  128.  * Perform the inverse DCT on one block of coefficients.
  129.  *
  130.  * A 2-D IDCT can be done by 1-D IDCT on each row
  131.  * followed by 1-D IDCT on each column.
  132.  */
  133.  
  134. GLOBAL void
  135. j_rev_dct (DCTBLOCK data)
  136. {
  137.   int pass, rowctr;
  138.   register DCTELEM *inptr, *outptr;
  139.   DCTBLOCK workspace;
  140.  
  141.   /* Each iteration of the inner loop performs one 8-point 1-D IDCT.
  142.    * It reads from a *row* of the input matrix and stores into a *column*
  143.    * of the output matrix.  In the first pass, we read from the data[] array
  144.    * and store into the local workspace[].  In the second pass, we read from
  145.    * the workspace[] array and store into data[], thus performing the
  146.    * equivalent of a columnar IDCT pass with no variable array indexing.
  147.    */
  148.  
  149.   inptr = data;            /* initialize pointers for first pass */
  150.   outptr = workspace;
  151.   for (pass = 1; pass >= 0; pass--) {
  152.     for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  153.       /* many tmps have nonoverlapping lifetime -- flashy register colourers
  154.        * should be able to do this lot very well
  155.        */
  156.       INT32 in0, in1, in2, in3, in4, in5, in6, in7;
  157.       INT32 tmp10, tmp11, tmp12, tmp13;
  158.       INT32 tmp20, tmp21, tmp22, tmp23;
  159.       INT32 tmp30, tmp31;
  160.       INT32 tmp40, tmp41, tmp42, tmp43;
  161.       INT32 tmp50, tmp51, tmp52, tmp53;
  162.       SHIFT_TEMPS
  163.     
  164.       in0 = inptr[0];
  165.       in1 = inptr[1];
  166.       in2 = inptr[2];
  167.       in3 = inptr[3];
  168.       in4 = inptr[4];
  169.       in5 = inptr[5];
  170.       in6 = inptr[6];
  171.       in7 = inptr[7];
  172.       
  173.       /* These values are scaled by DCT_SCALE */
  174.       
  175.       tmp10 = (in0 + in4) * COS_1_4;
  176.       tmp11 = (in0 - in4) * COS_1_4;
  177.       tmp12 = in2 * SIN_1_8 - in6 * COS_1_8;
  178.       tmp13 = in6 * SIN_1_8 + in2 * COS_1_8;
  179.       
  180.       tmp20 = tmp10 + tmp13;
  181.       tmp21 = tmp11 + tmp12;
  182.       tmp22 = tmp11 - tmp12;
  183.       tmp23 = tmp10 - tmp13;
  184.       
  185.       /* These values are scaled by OVERSCALE */
  186.       
  187.       tmp30 = UNFIXO((in3 + in5) * COS_1_4);
  188.       tmp31 = UNFIXO((in3 - in5) * COS_1_4);
  189.       
  190.       OVERSHIFT(in1);
  191.       OVERSHIFT(in7);
  192.       
  193.       tmp40 = in1 + tmp30;
  194.       tmp41 = in7 + tmp31;
  195.       tmp42 = in1 - tmp30;
  196.       tmp43 = in7 - tmp31;
  197.       
  198.       /* And these are scaled by DCT_SCALE */
  199.       
  200.       tmp50 = tmp40 * OCOS_1_16 + tmp41 * OSIN_1_16;
  201.       tmp51 = tmp40 * OSIN_1_16 - tmp41 * OCOS_1_16;
  202.       tmp52 = tmp42 * OCOS_5_16 + tmp43 * OSIN_5_16;
  203.       tmp53 = tmp42 * OSIN_5_16 - tmp43 * OCOS_5_16;
  204.       
  205.       outptr[        0] = (DCTELEM) UNFIXH(tmp20 + tmp50);
  206.       outptr[DCTSIZE  ] = (DCTELEM) UNFIXH(tmp21 + tmp53);
  207.       outptr[DCTSIZE*2] = (DCTELEM) UNFIXH(tmp22 + tmp52);
  208.       outptr[DCTSIZE*3] = (DCTELEM) UNFIXH(tmp23 + tmp51);
  209.       outptr[DCTSIZE*4] = (DCTELEM) UNFIXH(tmp23 - tmp51);
  210.       outptr[DCTSIZE*5] = (DCTELEM) UNFIXH(tmp22 - tmp52);
  211.       outptr[DCTSIZE*6] = (DCTELEM) UNFIXH(tmp21 - tmp53);
  212.       outptr[DCTSIZE*7] = (DCTELEM) UNFIXH(tmp20 - tmp50);
  213.       
  214.       inptr += DCTSIZE;        /* advance inptr to next row */
  215.       outptr++;            /* advance outptr to next column */
  216.     }
  217.     /* end of pass; in case it was pass 1, set up for pass 2 */
  218.     inptr = workspace;
  219.     outptr = data;
  220.   }
  221. }
  222.