Home | History | Annotate | Download | only in jpeg
      1 /*
      2  * jidctint.c
      3  *
      4  * Copyright (C) 1991-1998, 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 a slow-but-accurate integer implementation of the
      9  * inverse DCT (Discrete Cosine Transform).  In the IJG code, this routine
     10  * must also perform dequantization of the input coefficients.
     11  *
     12  * A 2-D IDCT can be done by 1-D IDCT on each column followed by 1-D IDCT
     13  * on each row (or vice versa, but it's more convenient to emit a row at
     14  * a time).  Direct algorithms are also available, but they are much more
     15  * complex and seem not to be any faster when reduced to code.
     16  *
     17  * This implementation is based on an algorithm described in
     18  *   C. Loeffler, A. Ligtenberg and G. Moschytz, "Practical Fast 1-D DCT
     19  *   Algorithms with 11 Multiplications", Proc. Int'l. Conf. on Acoustics,
     20  *   Speech, and Signal Processing 1989 (ICASSP '89), pp. 988-991.
     21  * The primary algorithm described there uses 11 multiplies and 29 adds.
     22  * We use their alternate method with 12 multiplies and 32 adds.
     23  * The advantage of this method is that no data path contains more than one
     24  * multiplication; this allows a very simple and accurate implementation in
     25  * scaled fixed-point arithmetic, with a minimal number of shifts.
     26  */
     27 
     28 #define JPEG_INTERNALS
     29 #include "jinclude.h"
     30 #include "jpeglib.h"
     31 #include "jdct.h"		/* Private declarations for DCT subsystem */
     32 
     33 #ifdef DCT_ISLOW_SUPPORTED
     34 
     35 
     36 /*
     37  * This module is specialized to the case DCTSIZE = 8.
     38  */
     39 
     40 #if DCTSIZE != 8
     41   Sorry, this code only copes with 8x8 DCTs. /* deliberate syntax err */
     42 #endif
     43 
     44 
     45 /*
     46  * The poop on this scaling stuff is as follows:
     47  *
     48  * Each 1-D IDCT step produces outputs which are a factor of sqrt(N)
     49  * larger than the true IDCT outputs.  The final outputs are therefore
     50  * a factor of N larger than desired; since N=8 this can be cured by
     51  * a simple right shift at the end of the algorithm.  The advantage of
     52  * this arrangement is that we save two multiplications per 1-D IDCT,
     53  * because the y0 and y4 inputs need not be divided by sqrt(N).
     54  *
     55  * We have to do addition and subtraction of the integer inputs, which
     56  * is no problem, and multiplication by fractional constants, which is
     57  * a problem to do in integer arithmetic.  We multiply all the constants
     58  * by CONST_SCALE and convert them to integer constants (thus retaining
     59  * CONST_BITS bits of precision in the constants).  After doing a
     60  * multiplication we have to divide the product by CONST_SCALE, with proper
     61  * rounding, to produce the correct output.  This division can be done
     62  * cheaply as a right shift of CONST_BITS bits.  We postpone shifting
     63  * as long as possible so that partial sums can be added together with
     64  * full fractional precision.
     65  *
     66  * The outputs of the first pass are scaled up by PASS1_BITS bits so that
     67  * they are represented to better-than-integral precision.  These outputs
     68  * require BITS_IN_JSAMPLE + PASS1_BITS + 3 bits; this fits in a 16-bit word
     69  * with the recommended scaling.  (To scale up 12-bit sample data further, an
     70  * intermediate INT32 array would be needed.)
     71  *
     72  * To avoid overflow of the 32-bit intermediate results in pass 2, we must
     73  * have BITS_IN_JSAMPLE + CONST_BITS + PASS1_BITS <= 26.  Error analysis
     74  * shows that the values given below are the most effective.
     75  */
     76 
     77 #if BITS_IN_JSAMPLE == 8
     78 #define CONST_BITS  13
     79 #define PASS1_BITS  2
     80 #else
     81 #define CONST_BITS  13
     82 #define PASS1_BITS  1		/* lose a little precision to avoid overflow */
     83 #endif
     84 
     85 /* Some C compilers fail to reduce "FIX(constant)" at compile time, thus
     86  * causing a lot of useless floating-point operations at run time.
     87  * To get around this we use the following pre-calculated constants.
     88  * If you change CONST_BITS you may want to add appropriate values.
     89  * (With a reasonable C compiler, you can just rely on the FIX() macro...)
     90  */
     91 
     92 #if CONST_BITS == 13
     93 #define FIX_0_298631336  ((INT32)  2446)	/* FIX(0.298631336) */
     94 #define FIX_0_390180644  ((INT32)  3196)	/* FIX(0.390180644) */
     95 #define FIX_0_541196100  ((INT32)  4433)	/* FIX(0.541196100) */
     96 #define FIX_0_765366865  ((INT32)  6270)	/* FIX(0.765366865) */
     97 #define FIX_0_899976223  ((INT32)  7373)	/* FIX(0.899976223) */
     98 #define FIX_1_175875602  ((INT32)  9633)	/* FIX(1.175875602) */
     99 #define FIX_1_501321110  ((INT32)  12299)	/* FIX(1.501321110) */
    100 #define FIX_1_847759065  ((INT32)  15137)	/* FIX(1.847759065) */
    101 #define FIX_1_961570560  ((INT32)  16069)	/* FIX(1.961570560) */
    102 #define FIX_2_053119869  ((INT32)  16819)	/* FIX(2.053119869) */
    103 #define FIX_2_562915447  ((INT32)  20995)	/* FIX(2.562915447) */
    104 #define FIX_3_072711026  ((INT32)  25172)	/* FIX(3.072711026) */
    105 #else
    106 #define FIX_0_298631336  FIX(0.298631336)
    107 #define FIX_0_390180644  FIX(0.390180644)
    108 #define FIX_0_541196100  FIX(0.541196100)
    109 #define FIX_0_765366865  FIX(0.765366865)
    110 #define FIX_0_899976223  FIX(0.899976223)
    111 #define FIX_1_175875602  FIX(1.175875602)
    112 #define FIX_1_501321110  FIX(1.501321110)
    113 #define FIX_1_847759065  FIX(1.847759065)
    114 #define FIX_1_961570560  FIX(1.961570560)
    115 #define FIX_2_053119869  FIX(2.053119869)
    116 #define FIX_2_562915447  FIX(2.562915447)
    117 #define FIX_3_072711026  FIX(3.072711026)
    118 #endif
    119 
    120 
    121 /* Multiply an INT32 variable by an INT32 constant to yield an INT32 result.
    122  * For 8-bit samples with the recommended scaling, all the variable
    123  * and constant values involved are no more than 16 bits wide, so a
    124  * 16x16->32 bit multiply can be used instead of a full 32x32 multiply.
    125  * For 12-bit samples, a full 32-bit multiplication will be needed.
    126  */
    127 
    128 #if BITS_IN_JSAMPLE == 8
    129 #define MULTIPLY(var,const)  MULTIPLY16C16(var,const)
    130 #else
    131 #define MULTIPLY(var,const)  ((var) * (const))
    132 #endif
    133 
    134 
    135 /* Dequantize a coefficient by multiplying it by the multiplier-table
    136  * entry; produce an int result.  In this module, both inputs and result
    137  * are 16 bits or less, so either int or short multiply will work.
    138  */
    139 
    140 #define DEQUANTIZE(coef,quantval)  (((ISLOW_MULT_TYPE) (coef)) * (quantval))
    141 
    142 
    143 /*
    144  * Perform dequantization and inverse DCT on one block of coefficients.
    145  */
    146 
    147 GLOBAL(void)
    148 jpeg_idct_islow (j_decompress_ptr cinfo, jpeg_component_info * compptr,
    149 		 JCOEFPTR coef_block,
    150 		 JSAMPARRAY output_buf, JDIMENSION output_col)
    151 {
    152   INT32 tmp0, tmp1, tmp2, tmp3;
    153   INT32 tmp10, tmp11, tmp12, tmp13;
    154   INT32 z1, z2, z3, z4, z5;
    155   JCOEFPTR inptr;
    156   ISLOW_MULT_TYPE * quantptr;
    157   int * wsptr;
    158   JSAMPROW outptr;
    159   JSAMPLE *range_limit = IDCT_range_limit(cinfo);
    160   int ctr;
    161   int workspace[DCTSIZE2];	/* buffers data between passes */
    162   SHIFT_TEMPS
    163 
    164   /* Pass 1: process columns from input, store into work array. */
    165   /* Note results are scaled up by sqrt(8) compared to a true IDCT; */
    166   /* furthermore, we scale the results by 2**PASS1_BITS. */
    167 
    168   inptr = coef_block;
    169   quantptr = (ISLOW_MULT_TYPE *) compptr->dct_table;
    170   wsptr = workspace;
    171   for (ctr = DCTSIZE; ctr > 0; ctr--) {
    172     /* Due to quantization, we will usually find that many of the input
    173      * coefficients are zero, especially the AC terms.  We can exploit this
    174      * by short-circuiting the IDCT calculation for any column in which all
    175      * the AC terms are zero.  In that case each output is equal to the
    176      * DC coefficient (with scale factor as needed).
    177      * With typical images and quantization tables, half or more of the
    178      * column DCT calculations can be simplified this way.
    179      */
    180 
    181     if (inptr[DCTSIZE*1] == 0 && inptr[DCTSIZE*2] == 0 &&
    182 	inptr[DCTSIZE*3] == 0 && inptr[DCTSIZE*4] == 0 &&
    183 	inptr[DCTSIZE*5] == 0 && inptr[DCTSIZE*6] == 0 &&
    184 	inptr[DCTSIZE*7] == 0) {
    185       /* AC terms all zero */
    186       int dcval = DEQUANTIZE(inptr[DCTSIZE*0], quantptr[DCTSIZE*0]) << PASS1_BITS;
    187 
    188       wsptr[DCTSIZE*0] = dcval;
    189       wsptr[DCTSIZE*1] = dcval;
    190       wsptr[DCTSIZE*2] = dcval;
    191       wsptr[DCTSIZE*3] = dcval;
    192       wsptr[DCTSIZE*4] = dcval;
    193       wsptr[DCTSIZE*5] = dcval;
    194       wsptr[DCTSIZE*6] = dcval;
    195       wsptr[DCTSIZE*7] = dcval;
    196 
    197       inptr++;			/* advance pointers to next column */
    198       quantptr++;
    199       wsptr++;
    200       continue;
    201     }
    202 
    203     /* Even part: reverse the even part of the forward DCT. */
    204     /* The rotator is sqrt(2)*c(-6). */
    205 
    206     z2 = DEQUANTIZE(inptr[DCTSIZE*2], quantptr[DCTSIZE*2]);
    207     z3 = DEQUANTIZE(inptr[DCTSIZE*6], quantptr[DCTSIZE*6]);
    208 
    209     z1 = MULTIPLY(z2 + z3, FIX_0_541196100);
    210     tmp2 = z1 + MULTIPLY(z3, - FIX_1_847759065);
    211     tmp3 = z1 + MULTIPLY(z2, FIX_0_765366865);
    212 
    213     z2 = DEQUANTIZE(inptr[DCTSIZE*0], quantptr[DCTSIZE*0]);
    214     z3 = DEQUANTIZE(inptr[DCTSIZE*4], quantptr[DCTSIZE*4]);
    215 
    216     tmp0 = (z2 + z3) << CONST_BITS;
    217     tmp1 = (z2 - z3) << CONST_BITS;
    218 
    219     tmp10 = tmp0 + tmp3;
    220     tmp13 = tmp0 - tmp3;
    221     tmp11 = tmp1 + tmp2;
    222     tmp12 = tmp1 - tmp2;
    223 
    224     /* Odd part per figure 8; the matrix is unitary and hence its
    225      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
    226      */
    227 
    228     tmp0 = DEQUANTIZE(inptr[DCTSIZE*7], quantptr[DCTSIZE*7]);
    229     tmp1 = DEQUANTIZE(inptr[DCTSIZE*5], quantptr[DCTSIZE*5]);
    230     tmp2 = DEQUANTIZE(inptr[DCTSIZE*3], quantptr[DCTSIZE*3]);
    231     tmp3 = DEQUANTIZE(inptr[DCTSIZE*1], quantptr[DCTSIZE*1]);
    232 
    233     z1 = tmp0 + tmp3;
    234     z2 = tmp1 + tmp2;
    235     z3 = tmp0 + tmp2;
    236     z4 = tmp1 + tmp3;
    237     z5 = MULTIPLY(z3 + z4, FIX_1_175875602); /* sqrt(2) * c3 */
    238 
    239     tmp0 = MULTIPLY(tmp0, FIX_0_298631336); /* sqrt(2) * (-c1+c3+c5-c7) */
    240     tmp1 = MULTIPLY(tmp1, FIX_2_053119869); /* sqrt(2) * ( c1+c3-c5+c7) */
    241     tmp2 = MULTIPLY(tmp2, FIX_3_072711026); /* sqrt(2) * ( c1+c3+c5-c7) */
    242     tmp3 = MULTIPLY(tmp3, FIX_1_501321110); /* sqrt(2) * ( c1+c3-c5-c7) */
    243     z1 = MULTIPLY(z1, - FIX_0_899976223); /* sqrt(2) * (c7-c3) */
    244     z2 = MULTIPLY(z2, - FIX_2_562915447); /* sqrt(2) * (-c1-c3) */
    245     z3 = MULTIPLY(z3, - FIX_1_961570560); /* sqrt(2) * (-c3-c5) */
    246     z4 = MULTIPLY(z4, - FIX_0_390180644); /* sqrt(2) * (c5-c3) */
    247 
    248     z3 += z5;
    249     z4 += z5;
    250 
    251     tmp0 += z1 + z3;
    252     tmp1 += z2 + z4;
    253     tmp2 += z2 + z3;
    254     tmp3 += z1 + z4;
    255 
    256     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
    257 
    258     wsptr[DCTSIZE*0] = (int) DESCALE(tmp10 + tmp3, CONST_BITS-PASS1_BITS);
    259     wsptr[DCTSIZE*7] = (int) DESCALE(tmp10 - tmp3, CONST_BITS-PASS1_BITS);
    260     wsptr[DCTSIZE*1] = (int) DESCALE(tmp11 + tmp2, CONST_BITS-PASS1_BITS);
    261     wsptr[DCTSIZE*6] = (int) DESCALE(tmp11 - tmp2, CONST_BITS-PASS1_BITS);
    262     wsptr[DCTSIZE*2] = (int) DESCALE(tmp12 + tmp1, CONST_BITS-PASS1_BITS);
    263     wsptr[DCTSIZE*5] = (int) DESCALE(tmp12 - tmp1, CONST_BITS-PASS1_BITS);
    264     wsptr[DCTSIZE*3] = (int) DESCALE(tmp13 + tmp0, CONST_BITS-PASS1_BITS);
    265     wsptr[DCTSIZE*4] = (int) DESCALE(tmp13 - tmp0, CONST_BITS-PASS1_BITS);
    266 
    267     inptr++;			/* advance pointers to next column */
    268     quantptr++;
    269     wsptr++;
    270   }
    271 
    272   /* Pass 2: process rows from work array, store into output array. */
    273   /* Note that we must descale the results by a factor of 8 == 2**3, */
    274   /* and also undo the PASS1_BITS scaling. */
    275 
    276   wsptr = workspace;
    277   for (ctr = 0; ctr < DCTSIZE; ctr++) {
    278     outptr = output_buf[ctr] + output_col;
    279     /* Rows of zeroes can be exploited in the same way as we did with columns.
    280      * However, the column calculation has created many nonzero AC terms, so
    281      * the simplification applies less often (typically 5% to 10% of the time).
    282      * On machines with very fast multiplication, it's possible that the
    283      * test takes more time than it's worth.  In that case this section
    284      * may be commented out.
    285      */
    286 
    287 #ifndef NO_ZERO_ROW_TEST
    288     if (wsptr[1] == 0 && wsptr[2] == 0 && wsptr[3] == 0 && wsptr[4] == 0 &&
    289 	wsptr[5] == 0 && wsptr[6] == 0 && wsptr[7] == 0) {
    290       /* AC terms all zero */
    291       JSAMPLE dcval = range_limit[(int) DESCALE((INT32) wsptr[0], PASS1_BITS+3)
    292 				  & RANGE_MASK];
    293 
    294       outptr[0] = dcval;
    295       outptr[1] = dcval;
    296       outptr[2] = dcval;
    297       outptr[3] = dcval;
    298       outptr[4] = dcval;
    299       outptr[5] = dcval;
    300       outptr[6] = dcval;
    301       outptr[7] = dcval;
    302 
    303       wsptr += DCTSIZE;		/* advance pointer to next row */
    304       continue;
    305     }
    306 #endif
    307 
    308     /* Even part: reverse the even part of the forward DCT. */
    309     /* The rotator is sqrt(2)*c(-6). */
    310 
    311     z2 = (INT32) wsptr[2];
    312     z3 = (INT32) wsptr[6];
    313 
    314     z1 = MULTIPLY(z2 + z3, FIX_0_541196100);
    315     tmp2 = z1 + MULTIPLY(z3, - FIX_1_847759065);
    316     tmp3 = z1 + MULTIPLY(z2, FIX_0_765366865);
    317 
    318     tmp0 = ((INT32) wsptr[0] + (INT32) wsptr[4]) << CONST_BITS;
    319     tmp1 = ((INT32) wsptr[0] - (INT32) wsptr[4]) << CONST_BITS;
    320 
    321     tmp10 = tmp0 + tmp3;
    322     tmp13 = tmp0 - tmp3;
    323     tmp11 = tmp1 + tmp2;
    324     tmp12 = tmp1 - tmp2;
    325 
    326     /* Odd part per figure 8; the matrix is unitary and hence its
    327      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
    328      */
    329 
    330     tmp0 = (INT32) wsptr[7];
    331     tmp1 = (INT32) wsptr[5];
    332     tmp2 = (INT32) wsptr[3];
    333     tmp3 = (INT32) wsptr[1];
    334 
    335     z1 = tmp0 + tmp3;
    336     z2 = tmp1 + tmp2;
    337     z3 = tmp0 + tmp2;
    338     z4 = tmp1 + tmp3;
    339     z5 = MULTIPLY(z3 + z4, FIX_1_175875602); /* sqrt(2) * c3 */
    340 
    341     tmp0 = MULTIPLY(tmp0, FIX_0_298631336); /* sqrt(2) * (-c1+c3+c5-c7) */
    342     tmp1 = MULTIPLY(tmp1, FIX_2_053119869); /* sqrt(2) * ( c1+c3-c5+c7) */
    343     tmp2 = MULTIPLY(tmp2, FIX_3_072711026); /* sqrt(2) * ( c1+c3+c5-c7) */
    344     tmp3 = MULTIPLY(tmp3, FIX_1_501321110); /* sqrt(2) * ( c1+c3-c5-c7) */
    345     z1 = MULTIPLY(z1, - FIX_0_899976223); /* sqrt(2) * (c7-c3) */
    346     z2 = MULTIPLY(z2, - FIX_2_562915447); /* sqrt(2) * (-c1-c3) */
    347     z3 = MULTIPLY(z3, - FIX_1_961570560); /* sqrt(2) * (-c3-c5) */
    348     z4 = MULTIPLY(z4, - FIX_0_390180644); /* sqrt(2) * (c5-c3) */
    349 
    350     z3 += z5;
    351     z4 += z5;
    352 
    353     tmp0 += z1 + z3;
    354     tmp1 += z2 + z4;
    355     tmp2 += z2 + z3;
    356     tmp3 += z1 + z4;
    357 
    358     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
    359 
    360     outptr[0] = range_limit[(int) DESCALE(tmp10 + tmp3,
    361 					  CONST_BITS+PASS1_BITS+3)
    362 			    & RANGE_MASK];
    363     outptr[7] = range_limit[(int) DESCALE(tmp10 - tmp3,
    364 					  CONST_BITS+PASS1_BITS+3)
    365 			    & RANGE_MASK];
    366     outptr[1] = range_limit[(int) DESCALE(tmp11 + tmp2,
    367 					  CONST_BITS+PASS1_BITS+3)
    368 			    & RANGE_MASK];
    369     outptr[6] = range_limit[(int) DESCALE(tmp11 - tmp2,
    370 					  CONST_BITS+PASS1_BITS+3)
    371 			    & RANGE_MASK];
    372     outptr[2] = range_limit[(int) DESCALE(tmp12 + tmp1,
    373 					  CONST_BITS+PASS1_BITS+3)
    374 			    & RANGE_MASK];
    375     outptr[5] = range_limit[(int) DESCALE(tmp12 - tmp1,
    376 					  CONST_BITS+PASS1_BITS+3)
    377 			    & RANGE_MASK];
    378     outptr[3] = range_limit[(int) DESCALE(tmp13 + tmp0,
    379 					  CONST_BITS+PASS1_BITS+3)
    380 			    & RANGE_MASK];
    381     outptr[4] = range_limit[(int) DESCALE(tmp13 - tmp0,
    382 					  CONST_BITS+PASS1_BITS+3)
    383 			    & RANGE_MASK];
    384 
    385     wsptr += DCTSIZE;		/* advance pointer to next row */
    386   }
    387 }
    388 
    389 #endif /* DCT_ISLOW_SUPPORTED */
    390