diff options
Diffstat (limited to 'converter/ppm/ppmtompeg/mfwddct.c')
| -rw-r--r-- | converter/ppm/ppmtompeg/mfwddct.c | 393 |
1 files changed, 393 insertions, 0 deletions
diff --git a/converter/ppm/ppmtompeg/mfwddct.c b/converter/ppm/ppmtompeg/mfwddct.c new file mode 100644 index 00000000..4643bf25 --- /dev/null +++ b/converter/ppm/ppmtompeg/mfwddct.c @@ -0,0 +1,393 @@ + +/* + * mfwddct.c (derived from jfwddct.c, which carries the following info) + * + * Copyright (C) 1991, 1992, 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 DCT (Discrete Cosine Transform) 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 "all.h" + +#include "dct.h" +#include "mtypes.h" +#include "opts.h" + +/* + * 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. + */ + +#define EIGHT_BIT_SAMPLES +#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 + +/* Prototypes */ +void reference_fwd_dct _ANSI_ARGS_((Block block, Block dest)); +void mp_fwd_dct_fast _ANSI_ARGS_((Block data2d, Block dest2d)); +void init_fdct _ANSI_ARGS_((void)); + +/* + * -------------------------------------------------------------- + * + * mp_fwd_dct_block2 -- + * + * Select the appropriate mp_fwd_dct routine + * + * Results: None + * + * Side effects: None + * + * -------------------------------------------------------------- + */ +extern boolean pureDCT; +void +mp_fwd_dct_block2(data, dest) + Block data, dest; +{ + if (pureDCT) reference_fwd_dct(data, dest); + else mp_fwd_dct_fast(data, dest); +} + +/* + * -------------------------------------------------------------- + * + * mp_fwd_dct_fast -- + * + * Perform the forward DCT on one block of samples. + * + * A 2-D DCT can be done by 1-D DCT on each row followed by 1-D DCT on each + * column. + * + * Results: None + * + * Side effects: Overwrites the input data + * + * -------------------------------------------------------------- + */ + +void +mp_fwd_dct_fast(data2d, dest2d) + Block data2d, dest2d; +{ + int16 *data = (int16 *) data2d; /* this algorithm wants + * a 1-d array */ + int16 *dest = (int16 *) dest2d; + int pass, rowctr; + register int16 *inptr, *outptr; + int16 workspace[DCTSIZE_SQ]; + SHIFT_TEMPS + +#ifdef ndef + { + int y; + + printf("fwd_dct (beforehand):\n"); + for (y = 0; y < 8; y++) + printf("%4d %4d %4d %4d %4d %4d %4d %4d\n", + data2d[y][0], data2d[y][1], + data2d[y][2], data2d[y][3], + data2d[y][4], data2d[y][5], + data2d[y][6], data2d[y][7]); + } +#endif + + /* + * Each iteration of the inner loop performs one 8-point 1-D DCT. It + * reads from a *row* of the input matrix and stores into a *column* + * of the output matrix. In the first pass, we read from the data[] + * array and store into the local workspace[]. In the second pass, + * we read from the workspace[] array and store into data[], thus + * performing the equivalent of a columnar DCT pass with no variable + * array indexing. + */ + + inptr = data; /* initialize pointers for first pass */ + outptr = workspace; + for (pass = 1; pass >= 0; pass--) { + for (rowctr = DCTSIZE - 1; rowctr >= 0; rowctr--) { + /* + * many tmps have nonoverlapping lifetime -- flashy + * register colorers should be able to do this lot + * very well + */ + int32 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7; + int32 tmp10, tmp11, tmp12, tmp13; + int32 tmp14, tmp15, tmp16, tmp17; + int32 tmp25, tmp26; + /* SHIFT_TEMPS */ + + /* temp0 through tmp7: -512 to +512 */ + /* if I-block, then -256 to +256 */ + tmp0 = inptr[7] + inptr[0]; + tmp1 = inptr[6] + inptr[1]; + tmp2 = inptr[5] + inptr[2]; + tmp3 = inptr[4] + inptr[3]; + tmp4 = inptr[3] - inptr[4]; + tmp5 = inptr[2] - inptr[5]; + tmp6 = inptr[1] - inptr[6]; + tmp7 = inptr[0] - inptr[7]; + + /* tmp10 through tmp13: -1024 to +1024 */ + /* if I-block, then -512 to +512 */ + tmp10 = tmp3 + tmp0; + tmp11 = tmp2 + tmp1; + tmp12 = tmp1 - tmp2; + tmp13 = tmp0 - tmp3; + + outptr[0] = (int16) UNFIXH((tmp10 + tmp11) * SIN_1_4); + outptr[DCTSIZE * 4] = (int16) UNFIXH((tmp10 - tmp11) * COS_1_4); + + outptr[DCTSIZE * 2] = (int16) UNFIXH(tmp13 * COS_1_8 + tmp12 * SIN_1_8); + outptr[DCTSIZE * 6] = (int16) UNFIXH(tmp13 * SIN_1_8 - tmp12 * COS_1_8); + + tmp16 = UNFIXO((tmp6 + tmp5) * SIN_1_4); + tmp15 = UNFIXO((tmp6 - tmp5) * COS_1_4); + + OVERSHIFT(tmp4); + OVERSHIFT(tmp7); + + /* + * tmp4, tmp7, tmp15, tmp16 are overscaled by + * OVERSCALE + */ + + tmp14 = tmp4 + tmp15; + tmp25 = tmp4 - tmp15; + tmp26 = tmp7 - tmp16; + tmp17 = tmp7 + tmp16; + + outptr[DCTSIZE] = (int16) UNFIXH(tmp17 * OCOS_1_16 + tmp14 * OSIN_1_16); + outptr[DCTSIZE * 7] = (int16) UNFIXH(tmp17 * OCOS_7_16 - tmp14 * OSIN_7_16); + outptr[DCTSIZE * 5] = (int16) UNFIXH(tmp26 * OCOS_5_16 + tmp25 * OSIN_5_16); + outptr[DCTSIZE * 3] = (int16) UNFIXH(tmp26 * OCOS_3_16 - tmp25 * OSIN_3_16); + + inptr += DCTSIZE; /* advance inptr to next row */ + outptr++; /* advance outptr to next column */ + } + /* end of pass; in case it was pass 1, set up for pass 2 */ + inptr = workspace; + outptr = dest; + } +#ifdef ndef + { + int y; + + printf("fwd_dct (afterward):\n"); + for (y = 0; y < 8; y++) + printf("%4d %4d %4d %4d %4d %4d %4d %4d\n", + dest2d[y][0], dest2d[y][1], + dest2d[y][2], dest2d[y][3], + dest2d[y][4], dest2d[y][5], + dest2d[y][6], dest2d[y][7]); + } +#endif +} + + +/* Modifies from the MPEG2 verification coder */ +/* fdctref.c, forward discrete cosine transform, double precision */ + +/* Copyright (C) 1994, MPEG Software Simulation Group. All Rights Reserved. */ + +/* + * Disclaimer of Warranty + * + * These software programs are available to the user without any license fee or + * royalty on an "as is" basis. The MPEG Software Simulation Group disclaims + * any and all warranties, whether express, implied, or statuary, including any + * implied warranties or merchantability or of fitness for a particular + * purpose. In no event shall the copyright-holder be liable for any + * incidental, punitive, or consequential damages of any kind whatsoever + * arising from the use of these programs. + * + * This disclaimer of warranty extends to the user of these programs and user's + * customers, employees, agents, transferees, successors, and assigns. + * + * The MPEG Software Simulation Group does not represent or warrant that the + * programs furnished hereunder are free of infringement of any third-party + * patents. + * + * Commercial implementations of MPEG-1 and MPEG-2 video, including shareware, + * are subject to royalty fees to patent holders. Many of these patents are + * general enough such that they are unavoidable regardless of implementation + * design. + * + */ + +#ifndef PI +#ifdef M_PI +#define PI M_PI +#else +#define PI 3.14159265358979323846 +#endif +#endif + +/* private data */ +static double trans_coef[8][8]; /* transform coefficients */ + +void init_fdct() +{ + int i, j; + double s; + + for (i=0; i<8; i++) + { + s = (i==0) ? sqrt(0.125) : 0.5; + + for (j=0; j<8; j++) + trans_coef[i][j] = s * cos((PI/8.0)*i*(j+0.5)); + } +} + +void reference_fwd_dct(block, dest) +Block block, dest; +{ + int i, j, k; + double s; + double tmp[64]; + + if (DoLaplace) { + LaplaceNum++; + } + + for (i=0; i<8; i++) + for (j=0; j<8; j++) + { + s = 0.0; + + for (k=0; k<8; k++) + s += trans_coef[j][k] * block[i][k]; + + tmp[8*i+j] = s; + } + + for (i=0; i<8; i++) + for (j=0; j<8; j++) + { + s = 0.0; + + for (k=0; k<8; k++) + s += trans_coef[i][k] * tmp[8*k+j]; + + if (collect_quant) { + fprintf(collect_quant_fp, "%d %f\n", 8*i+j, s); + } + if (DoLaplace) { + L1[LaplaceCnum][i*8+j] += s*s; + L2[LaplaceCnum][i*8+j] += s; + } + + + dest[i][j] = (int)floor(s+0.499999); + /* + * reason for adding 0.499999 instead of 0.5: + * s is quite often x.5 (at least for i and/or j = 0 or 4) + * and setting the rounding threshold exactly to 0.5 leads to an + * extremely high arithmetic implementation dependency of the result; + * s being between x.5 and x.500001 (which is now incorrectly rounded + * downwards instead of upwards) is assumed to occur less often + * (if at all) + */ + } +} |