dct.c

quicktime linux播放器v1

/**************************************************************************
 *                                                                        *
 * This code is developed by Adam Li.  This software is an                *
 * implementation of a part of one or more MPEG-4 Video tools as          *
 * specified in ISO/IEC 14496-2 standard.  Those intending to use this    *
 * software module in hardware or software products are advised that its  *
 * use may infringe existing patents or copyrights, and any such use      *
 * would be at such party's own risk.  The original developer of this     *
 * software module and his/her company, and subsequent editors and their  *
 * companies (including Project Mayo), will have no liability for use of  *
 * this software or modifications or derivatives thereof.                 *
 *                                                                        *
 * Project Mayo gives users of the Codec a license to this software       *
 * module or modifications thereof for use in hardware or software        *
 * products claiming conformance to the MPEG-4 Video Standard as          *
 * described in the Open DivX license.                                    *
 *                                                                        *
 * The complete Open DivX license can be found at                         *
 * http://www.projectmayo.com/opendivx/license.php .                      *
 *                                                                        *
 **************************************************************************/

/**************************************************************************
 *
 *  text_dct.c
 *
 *  Copyright (C) 2001  Project Mayo
 *
 *  Adam Li
 *
 *  DivX Advance Research Center <[e-mail]>
 *
 **************************************************************************/

/* This file contains some functions for fDCT/iDCT transformation.        */

/* Some codes of this project come from SSG MPEG-2 implementation.        */

/* Please see seperate acknowledgement file for a list of contributors.   */

#include <math.h>
#include "enc_portab.h"

/* The first part of it is for the forward DCT */

#ifndef PI
# ifdef M_PI
#  define PI M_PI
# else
#  define PI 3.14159265358979323846
# endif
#endif

/* private data */
static double c[8][8];				 /* 

						    transform coefficients */

void fdct_enc(int16_t *block)
{
    int i, j, k;
    double s;
    double tmp[64];

    for (i = 0; i < 8; i++)
	for (j = 0; j < 8; j++)
	{
	    s = 0.0;

	    for (k = 0; k < 8; k++)
		s += c[j][k] * block[8 * i + k];

	    tmp[8 * i + j] = s;
	}

    for (j = 0; j < 8; j++)
	for (i = 0; i < 8; i++)
	{
	    s = 0.0;

	    for (k = 0; k < 8; k++)
		s += c[i][k] * tmp[8 * k + j];

	    block[8 * 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)
	     */
	}
}

void init_fdct_enc()
{
    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++)
	    c[i][j] = s * cos((PI / 8.0) * i * (j + 0.5));
    }
}




/* the second part of it is for the inverse DCT */

#define W1 2841					 /* 
						    2048*sqrt(2)*cos(1*pi/16) */
#define W2 2676					 /* 
						    2048*sqrt(2)*cos(2*pi/16) */
#define W3 2408					 /* 
						    2048*sqrt(2)*cos(3*pi/16) */
#define W5 1609					 /* 
						    2048*sqrt(2)*cos(5*pi/16) */
#define W6 1108					 /* 
						    2048*sqrt(2)*cos(6*pi/16) */
#define W7 565					 /* 
						    2048*sqrt(2)*cos(7*pi/16) */

/* private data */
static int16_t iclip[1024];			 /* 

						    clipping table */
static int16_t *iclp;

/* private prototypes */
static void idctrow_enc(int16_t *blk);
static void idctcol_enc(int16_t *blk);

/* two dimensional inverse discrete cosine transform */
void idct_enc(int16_t *block)
{
    int i;

    for (i = 0; i < 8; i++)
	idctrow_enc(block + 8 * i);

    for (i = 0; i < 8; i++)
	idctcol_enc(block + i);
}

void init_idct_enc()
{
    int i;

    iclp = iclip + 512;
    for (i = -512; i < 512; i++)
	iclp[i] = (i < -256) ? -256 : ((i > 255) ? 255 : i);
}

/* row (horizontal) IDCT
 *
 *           7                       pi         1
 * dst[k] = sum c[l] * src[l] * cos( -- * ( k + - ) * l )
 *          l=0                      8          2
 *
 * where: c[0]    = 128
 *        c[1..7] = 128*sqrt(2)
 */

static void idctrow_enc(blk)
     int16_t *blk;
{
    int x0, x1, x2, x3, x4, x5, x6, x7, x8;

    /* 
       shortcut */
    if (!((x1 = blk[4] << 11) | (x2 = blk[6]) | (x3 = blk[2]) |
	  (x4 = blk[1]) | (x5 = blk[7]) | (x6 = blk[5]) | (x7 = blk[3])))
    {
	blk[0] = blk[1] = blk[2] = blk[3] = blk[4] = blk[5] = blk[6] =
	    blk[7] = blk[0] << 3;
	return;
    }

    x0 = (blk[0] << 11) + 128;			 /* 
						    for proper rounding in the fourth stage */

    /* 
       first stage */
    x8 = W7 * (x4 + x5);
    x4 = x8 + (W1 - W7) * x4;
    x5 = x8 - (W1 + W7) * x5;
    x8 = W3 * (x6 + x7);
    x6 = x8 - (W3 - W5) * x6;
    x7 = x8 - (W3 + W5) * x7;

    /* 
       second stage */
    x8 = x0 + x1;
    x0 -= x1;
    x1 = W6 * (x3 + x2);
    x2 = x1 - (W2 + W6) * x2;
    x3 = x1 + (W2 - W6) * x3;
    x1 = x4 + x6;
    x4 -= x6;
    x6 = x5 + x7;
    x5 -= x7;

    /* 
       third stage */
    x7 = x8 + x3;
    x8 -= x3;
    x3 = x0 + x2;
    x0 -= x2;
    x2 = (181 * (x4 + x5) + 128) >> 8;
    x4 = (181 * (x4 - x5) + 128) >> 8;

    /* 
       fourth stage */
    blk[0] = (x7 + x1) >> 8;
    blk[1] = (x3 + x2) >> 8;
    blk[2] = (x0 + x4) >> 8;
    blk[3] = (x8 + x6) >> 8;
    blk[4] = (x8 - x6) >> 8;
    blk[5] = (x0 - x4) >> 8;
    blk[6] = (x3 - x2) >> 8;
    blk[7] = (x7 - x1) >> 8;
}

/* column (vertical) IDCT
 *
 *             7                         pi         1
 * dst[8*k] = sum c[l] * src[8*l] * cos( -- * ( k + - ) * l )
 *            l=0                        8          2
 *
 * where: c[0]    = 1/1024
 *        c[1..7] = (1/1024)*sqrt(2)
 */
static void idctcol_enc(blk)
     int16_t *blk;
{
    int x0, x1, x2, x3, x4, x5, x6, x7, x8;

    /* 
       shortcut */
    if (!((x1 = (blk[8 * 4] << 8)) | (x2 = blk[8 * 6]) | (x3 = blk[8 * 2]) |
	  (x4 = blk[8 * 1]) | (x5 = blk[8 * 7]) | (x6 = blk[8 * 5]) | (x7 =
								       blk[8 *
									   3])))
    {
	blk[8 * 0] = blk[8 * 1] = blk[8 * 2] = blk[8 * 3] = blk[8 * 4] =
	    blk[8 * 5] = blk[8 * 6] = blk[8 * 7] =
	    iclp[(blk[8 * 0] + 32) >> 6];
	return;
    }

    x0 = (blk[8 * 0] << 8) + 8192;

    /* 
       first stage */
    x8 = W7 * (x4 + x5) + 4;
    x4 = (x8 + (W1 - W7) * x4) >> 3;
    x5 = (x8 - (W1 + W7) * x5) >> 3;
    x8 = W3 * (x6 + x7) + 4;
    x6 = (x8 - (W3 - W5) * x6) >> 3;
    x7 = (x8 - (W3 + W5) * x7) >> 3;

    /* 
       second stage */
    x8 = x0 + x1;
    x0 -= x1;
    x1 = W6 * (x3 + x2) + 4;
    x2 = (x1 - (W2 + W6) * x2) >> 3;
    x3 = (x1 + (W2 - W6) * x3) >> 3;
    x1 = x4 + x6;
    x4 -= x6;
    x6 = x5 + x7;
    x5 -= x7;

    /* 
       third stage */
    x7 = x8 + x3;
    x8 -= x3;
    x3 = x0 + x2;
    x0 -= x2;
    x2 = (181 * (x4 + x5) + 128) >> 8;
    x4 = (181 * (x4 - x5) + 128) >> 8;

    /* 
       fourth stage */
    blk[8 * 0] = iclp[(x7 + x1) >> 14];
    blk[8 * 1] = iclp[(x3 + x2) >> 14];
    blk[8 * 2] = iclp[(x0 + x4) >> 14];
    blk[8 * 3] = iclp[(x8 + x6) >> 14];
    blk[8 * 4] = iclp[(x8 - x6) >> 14];
    blk[8 * 5] = iclp[(x0 - x4) >> 14];
    blk[8 * 6] = iclp[(x3 - x2) >> 14];
    blk[8 * 7] = iclp[(x7 - x1) >> 14];
}