зеркало из https://github.com/mozilla/mozjpeg.git
230 строки
7.1 KiB
C
230 строки
7.1 KiB
C
/*
|
|
* jrevdct.c
|
|
*
|
|
* Copyright (C) 1991, Thomas G. Lane.
|
|
* This file is part of the Independent JPEG Group's software.
|
|
* For conditions of distribution and use, see the accompanying README file.
|
|
*
|
|
* This file contains the basic inverse-DCT transformation subroutine.
|
|
*
|
|
* This implementation is based on Appendix A.2 of the book
|
|
* "Discrete Cosine Transform---Algorithms, Advantages, Applications"
|
|
* by K.R. Rao and P. Yip (Academic Press, Inc, London, 1990).
|
|
* It uses scaled fixed-point arithmetic instead of floating point.
|
|
*/
|
|
|
|
#include "jinclude.h"
|
|
|
|
|
|
/* We assume that right shift corresponds to signed division by 2 with
|
|
* rounding towards minus infinity. This is correct for typical "arithmetic
|
|
* shift" instructions that shift in copies of the sign bit. But some
|
|
* C compilers implement >> with an unsigned shift. For these machines you
|
|
* must define RIGHT_SHIFT_IS_UNSIGNED.
|
|
* RIGHT_SHIFT provides a signed right shift of an INT32 quantity.
|
|
* It is only applied with constant shift counts.
|
|
*/
|
|
|
|
#ifdef RIGHT_SHIFT_IS_UNSIGNED
|
|
#define SHIFT_TEMPS INT32 shift_temp;
|
|
#define RIGHT_SHIFT(x,shft) \
|
|
((shift_temp = (x)) < 0 ? \
|
|
(shift_temp >> (shft)) | ((~0) << (32-(shft))) : \
|
|
(shift_temp >> (shft)))
|
|
#else
|
|
#define SHIFT_TEMPS
|
|
#define RIGHT_SHIFT(x,shft) ((x) >> (shft))
|
|
#endif
|
|
|
|
|
|
/* The poop on this scaling stuff is as follows:
|
|
*
|
|
* We have to do addition and subtraction of the integer inputs, which
|
|
* is no problem, and multiplication by fractional constants, which is
|
|
* a problem to do in integer arithmetic. We multiply all the constants
|
|
* by DCT_SCALE and convert them to integer constants (thus retaining
|
|
* LG2_DCT_SCALE bits of precision in the constants). After doing a
|
|
* multiplication we have to divide the product by DCT_SCALE, with proper
|
|
* rounding, to produce the correct output. The division can be implemented
|
|
* cheaply as a right shift of LG2_DCT_SCALE bits. The DCT equations also
|
|
* specify an additional division by 2 on the final outputs; this can be
|
|
* folded into the right-shift by shifting one more bit (see UNFIXH).
|
|
*
|
|
* If you are planning to recode this in assembler, you might want to set
|
|
* LG2_DCT_SCALE to 15. This loses a bit of precision, but then all the
|
|
* multiplications are between 16-bit quantities (given 8-bit JSAMPLEs!)
|
|
* so you could use a signed 16x16=>32 bit multiply instruction instead of
|
|
* full 32x32 multiply. Unfortunately there's no way to describe such a
|
|
* multiply portably in C, so we've gone for the extra bit of accuracy here.
|
|
*/
|
|
|
|
#ifdef EIGHT_BIT_SAMPLES
|
|
#define LG2_DCT_SCALE 16
|
|
#else
|
|
#define LG2_DCT_SCALE 15 /* lose a little precision to avoid overflow */
|
|
#endif
|
|
|
|
#define ONE ((INT32) 1)
|
|
|
|
#define DCT_SCALE (ONE << LG2_DCT_SCALE)
|
|
|
|
/* In some places we shift the inputs left by a couple more bits, */
|
|
/* so that they can be added to fractional results without too much */
|
|
/* loss of precision. */
|
|
#define LG2_OVERSCALE 2
|
|
#define OVERSCALE (ONE << LG2_OVERSCALE)
|
|
#define OVERSHIFT(x) ((x) <<= LG2_OVERSCALE)
|
|
|
|
/* Scale a fractional constant by DCT_SCALE */
|
|
#define FIX(x) ((INT32) ((x) * DCT_SCALE + 0.5))
|
|
|
|
/* Scale a fractional constant by DCT_SCALE/OVERSCALE */
|
|
/* Such a constant can be multiplied with an overscaled input */
|
|
/* to produce something that's scaled by DCT_SCALE */
|
|
#define FIXO(x) ((INT32) ((x) * DCT_SCALE / OVERSCALE + 0.5))
|
|
|
|
/* Descale and correctly round a value that's scaled by DCT_SCALE */
|
|
#define UNFIX(x) RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1)), LG2_DCT_SCALE)
|
|
|
|
/* Same with an additional division by 2, ie, correctly rounded UNFIX(x/2) */
|
|
#define UNFIXH(x) RIGHT_SHIFT((x) + (ONE << LG2_DCT_SCALE), LG2_DCT_SCALE+1)
|
|
|
|
/* Take a value scaled by DCT_SCALE and round to integer scaled by OVERSCALE */
|
|
#define UNFIXO(x) RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1-LG2_OVERSCALE)),\
|
|
LG2_DCT_SCALE-LG2_OVERSCALE)
|
|
|
|
/* Here are the constants we need */
|
|
/* SIN_i_j is sine of i*pi/j, scaled by DCT_SCALE */
|
|
/* COS_i_j is cosine of i*pi/j, scaled by DCT_SCALE */
|
|
|
|
#define SIN_1_4 FIX(0.707106781)
|
|
#define COS_1_4 SIN_1_4
|
|
|
|
#define SIN_1_8 FIX(0.382683432)
|
|
#define COS_1_8 FIX(0.923879533)
|
|
#define SIN_3_8 COS_1_8
|
|
#define COS_3_8 SIN_1_8
|
|
|
|
#define SIN_1_16 FIX(0.195090322)
|
|
#define COS_1_16 FIX(0.980785280)
|
|
#define SIN_7_16 COS_1_16
|
|
#define COS_7_16 SIN_1_16
|
|
|
|
#define SIN_3_16 FIX(0.555570233)
|
|
#define COS_3_16 FIX(0.831469612)
|
|
#define SIN_5_16 COS_3_16
|
|
#define COS_5_16 SIN_3_16
|
|
|
|
/* OSIN_i_j is sine of i*pi/j, scaled by DCT_SCALE/OVERSCALE */
|
|
/* OCOS_i_j is cosine of i*pi/j, scaled by DCT_SCALE/OVERSCALE */
|
|
|
|
#define OSIN_1_4 FIXO(0.707106781)
|
|
#define OCOS_1_4 OSIN_1_4
|
|
|
|
#define OSIN_1_8 FIXO(0.382683432)
|
|
#define OCOS_1_8 FIXO(0.923879533)
|
|
#define OSIN_3_8 OCOS_1_8
|
|
#define OCOS_3_8 OSIN_1_8
|
|
|
|
#define OSIN_1_16 FIXO(0.195090322)
|
|
#define OCOS_1_16 FIXO(0.980785280)
|
|
#define OSIN_7_16 OCOS_1_16
|
|
#define OCOS_7_16 OSIN_1_16
|
|
|
|
#define OSIN_3_16 FIXO(0.555570233)
|
|
#define OCOS_3_16 FIXO(0.831469612)
|
|
#define OSIN_5_16 OCOS_3_16
|
|
#define OCOS_5_16 OSIN_3_16
|
|
|
|
|
|
/*
|
|
* Perform a 1-dimensional inverse DCT.
|
|
* Note that this code is specialized to the case DCTSIZE = 8.
|
|
*/
|
|
|
|
INLINE
|
|
LOCAL void
|
|
fast_idct_8 (DCTELEM *in, int stride)
|
|
{
|
|
/* many tmps have nonoverlapping lifetime -- flashy register colourers
|
|
* should be able to do this lot very well
|
|
*/
|
|
INT32 in0, in1, in2, in3, in4, in5, in6, in7;
|
|
INT32 tmp10, tmp11, tmp12, tmp13;
|
|
INT32 tmp20, tmp21, tmp22, tmp23;
|
|
INT32 tmp30, tmp31;
|
|
INT32 tmp40, tmp41, tmp42, tmp43;
|
|
INT32 tmp50, tmp51, tmp52, tmp53;
|
|
SHIFT_TEMPS
|
|
|
|
in0 = in[ 0];
|
|
in1 = in[stride ];
|
|
in2 = in[stride*2];
|
|
in3 = in[stride*3];
|
|
in4 = in[stride*4];
|
|
in5 = in[stride*5];
|
|
in6 = in[stride*6];
|
|
in7 = in[stride*7];
|
|
|
|
/* These values are scaled by DCT_SCALE */
|
|
|
|
tmp10 = (in0 + in4) * COS_1_4;
|
|
tmp11 = (in0 - in4) * COS_1_4;
|
|
tmp12 = in2 * SIN_1_8 - in6 * COS_1_8;
|
|
tmp13 = in6 * SIN_1_8 + in2 * COS_1_8;
|
|
|
|
tmp20 = tmp10 + tmp13;
|
|
tmp21 = tmp11 + tmp12;
|
|
tmp22 = tmp11 - tmp12;
|
|
tmp23 = tmp10 - tmp13;
|
|
|
|
/* These values are scaled by OVERSCALE */
|
|
|
|
tmp30 = UNFIXO((in3 + in5) * COS_1_4);
|
|
tmp31 = UNFIXO((in3 - in5) * COS_1_4);
|
|
|
|
OVERSHIFT(in1);
|
|
OVERSHIFT(in7);
|
|
|
|
tmp40 = in1 + tmp30;
|
|
tmp41 = in7 + tmp31;
|
|
tmp42 = in1 - tmp30;
|
|
tmp43 = in7 - tmp31;
|
|
|
|
/* And these are scaled by DCT_SCALE */
|
|
|
|
tmp50 = tmp40 * OCOS_1_16 + tmp41 * OSIN_1_16;
|
|
tmp51 = tmp40 * OSIN_1_16 - tmp41 * OCOS_1_16;
|
|
tmp52 = tmp42 * OCOS_5_16 + tmp43 * OSIN_5_16;
|
|
tmp53 = tmp42 * OSIN_5_16 - tmp43 * OCOS_5_16;
|
|
|
|
in[ 0] = (DCTELEM) UNFIXH(tmp20 + tmp50);
|
|
in[stride ] = (DCTELEM) UNFIXH(tmp21 + tmp53);
|
|
in[stride*2] = (DCTELEM) UNFIXH(tmp22 + tmp52);
|
|
in[stride*3] = (DCTELEM) UNFIXH(tmp23 + tmp51);
|
|
in[stride*4] = (DCTELEM) UNFIXH(tmp23 - tmp51);
|
|
in[stride*5] = (DCTELEM) UNFIXH(tmp22 - tmp52);
|
|
in[stride*6] = (DCTELEM) UNFIXH(tmp21 - tmp53);
|
|
in[stride*7] = (DCTELEM) UNFIXH(tmp20 - tmp50);
|
|
}
|
|
|
|
|
|
/*
|
|
* Perform the inverse DCT on one block of coefficients.
|
|
*
|
|
* A 2-D IDCT can be done by 1-D IDCT on each row
|
|
* followed by 1-D IDCT on each column.
|
|
*/
|
|
|
|
GLOBAL void
|
|
j_rev_dct (DCTBLOCK data)
|
|
{
|
|
int i;
|
|
|
|
for (i = 0; i < DCTSIZE; i++)
|
|
fast_idct_8(data+i*DCTSIZE, 1);
|
|
|
|
for (i = 0; i < DCTSIZE; i++)
|
|
fast_idct_8(data+i, DCTSIZE);
|
|
}
|