dct.c

h.263+ vc++商业源代码

#include "sim.h"

#define C1 2009L /* 2048*cos(1*pi/16) */
#define C2 2676L /* 2048*sqrt(2)*cos(2*pi/16) */
#define C3 1703L /* 2048*cos(3*pi/16) */
#define C4 1448L /* 2048*cos(4*pi/16) */
#define C5 1138L /* 2048*cos(5*pi/16) */
#define C6 1108L /* 2048*sqrt(2)*cos(6*pi/16) */
#define C7 400L /* 2048*cos(7*pi/16) */
#define W1 -1609L /* (C7-C1) */
#define W2 -2408L /*-(C7+C1) */
#define W3 -565L /*(C5-C3) */
#define W4 -2841L /*-(C5+C3) */
#define W5 1567L  /*(C2-C6) */
#define W6 -3784L  /*-(C2+C6) */
#define W7 1.4142143562 

#define Q1 2841L /* 2048*sqrt(2)*cos(1*pi/16) */
#define Q2 2676L /* 2048*sqrt(2)*cos(2*pi/16) */
#define Q3 2408L /* 2048*sqrt(2)*cos(3*pi/16) */
#define Q5 1609L /* 2048*sqrt(2)*cos(5*pi/16) */
#define Q6 1108L /* 2048*sqrt(2)*cos(6*pi/16) */
#define Q7 565L  /* 2048*sqrt(2)*cos(7*pi/16) */

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

unsigned char zig_zag[64] =
{
  0,1,8,16,9,2,3,10,17,24,32,25,18,11,4,5,
  12,19,26,33,40,48,41,34,27,20,13,6,7,14,21,28,
  35,42,49,56,57,50,43,36,29,22,15,23,30,37,44,51,
  58,59,52,45,38,31,39,46,53,60,61,54,47,55,62,63
};

/***********************************************************
 * 
 *    	Name:  		dct
 *      Description: 	fast discrete cosine transform
 * 			based on Loeffler's algorithem 
 * 			(Conf, Int'l, IEEE, Assp, 1989)
 * 			with 11 multiplies, 29 addition
 *   	Input:		pointer to coef
 *	Resturn: 	none
 *      Side effects:   none  
 *
 ***********************************************************/

void dctrow(int *);
void dctcol(int *);

void dctrow(int *blk)
{
 long x0, x1, x2, x3, x4, x5, x6, x7, x8;
 long y0, y1, y2, y3, y4, y5, y6, y7;

 /* stage 1 */
 x0 = blk[0]+blk[7];
 x7 = blk[0]-blk[7];
 x1 = blk[1]+blk[6];
 x6 = blk[1]-blk[6];
 x2 = blk[2]+blk[5];
 x5 = blk[2]-blk[5];
 x3 = blk[3]+blk[4];
 x4 = blk[3]-blk[4];

 /* stage 2 */
 y0 = x0+x3;
 y3 = x0-x3;
 y1 = x1+x2;
 y2 = x1-x2;
 x8 = C3*(x4+x7);
 y4 = W3*x7+x8;
 y7 = W4*x4+x8;
 x8 = C1*(x5+x6);
 y5 = W1*x6+x8;
 y6 = W2*x5+x8;

 /* stage 3 */
 x4 = y4+y6;
 x6 = y4-y6;
 x5 = -y5+y7;
 x7 = y5+y7;

 /* stage 4 */
 blk[0] = (int)(y0+y1);
 blk[4] = (int)(y0-y1);
 x8 = C6*(y3+y2);
 blk[2] = (int)((W5*y3+x8+1024L)>>11);
 blk[6] = (int)((W6*y2+x8+1024L)>>11);
 blk[7] = (int)((-x4+x7+1024L)>>11);
 blk[1] = (int)((x4+x7+1024L)>>11);
 blk[3] = (int)((long)(W7*x5+1024L)>>11);
 blk[5] = (int)((long)(W7*x6+1024L)>>11);
}

void dctcol(int *blk)
{
 long x0, x1, x2, x3, x4, x5, x6, x7, x8;
 long y0, y1, y2, y3, y4, y5, y6, y7;
 
 /* stage 1 */
 x0 = blk[0<<3]+blk[7<<3];
 x7 = blk[0<<3]-blk[7<<3];
 x1 = blk[1<<3]+blk[6<<3];
 x6 = blk[1<<3]-blk[6<<3];
 x2 = blk[2<<3]+blk[5<<3];
 x5 = blk[2<<3]-blk[5<<3];
 x3 = blk[3<<3]+blk[4<<3];
 x4 = blk[3<<3]-blk[4<<3];
 
 /* stage 2 */
 y0 = x0+x3;
 y3 = x0-x3;
 y1 = x1+x2;
 y2 = x1-x2;
 x8 = C3*(x4+x7);
 y4 = W3*x7+x8;
 y7 = W4*x4+x8;
 x8 = C1*(x5+x6);
 y5 = W1*x6+x8;
 y6 = W2*x5+x8;
 
 /* stage 3 */
 x4 = y4+y6;
 x6 = y4-y6;
 x5 = -y5+y7;
 x7 = y5+y7;
 
 /* stage 4 */
 blk[0<<3] = (int)(y0+y1);
 blk[4<<3] = (int)(y0-y1);
 x8 = C6*(y3+y2);
 blk[2<<3] = (int)((W5*y3+x8+1024)>>11);
 blk[6<<3] = (int)((W6*y2+x8+1024)>>11);
 blk[7<<3] = (int)((-x4+x7+1024)>>11);
 blk[1<<3] = (int)((x4+x7+1024)>>11);
 blk[3<<3] = (int)((long)(x5*W7+1024L)>>11);
 blk[5<<3] = (int)((long)(x6*W7+1024L)>>11);
}


void dct(int *tt)
{
  int i;

  for (i=0; i<64; i+=8)
    dctrow(tt+i);

  for (i=0; i<8; i++)
    dctcol(tt+i);

  for(i=0;i<64;i++)
    tt[i] = (tt[i]+4)>>3;
}





/**********************************************************/
/* inverse two dimensional DCT, Chen-Wang algorithm       */
/* (cf. IEEE ASSP-32, pp. 803-816, Aug. 1984)             */
/* 32-bit integer arithmetic (8 bit coefficients)         */
/* 11 mults, 29 adds per DCT                              */
/**********************************************************/
/* coefficients extended to 12 bit for IEEE1180-1990      */
/* compliance                                             */
/**********************************************************/

/* this code assumes >> to be a two's-complement arithmetic */
/* right shift: (-2)>>1 == -1 , (-3)>>1 == -2               */
/* global declarations */

/* private data */
static int iclip[1024]; /* clipping table */
static int *iclp;

/* private prototypes */
static void idctrow (int *);
static void idctcol (int *);

/* 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(int *blk)
{
  long x0, x1, x2, x3, x4, x5, x6, x7, x8;

  /* intcut */
  if (!((x1 = ((long)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]=(int)(((long)blk[0])<<3);
    return;
  }

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

  /* first stage */
  x8 = Q7*(x4+x5);
  x4 = x8 + (Q1-Q7)*x4;
  x5 = x8 - (Q1+Q7)*x5;
  x8 = Q3*(x6+x7);
  x6 = x8 - (Q3-Q5)*x6;
  x7 = x8 - (Q3+Q5)*x7;
  
  /* second stage */
  x8 = x0 + x1;
  x0 -= x1;
  x1 = Q6*(x3+x2);
  x2 = x1 - (Q2+Q6)*x2;
  x3 = x1 + (Q2-Q6)*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] = (int)((x7+x1)>>8);
  blk[1] = (int)((x3+x2)>>8);
  blk[2] = (int)((x0+x4)>>8);
  blk[3] = (int)((x8+x6)>>8);
  blk[4] = (int)((x8-x6)>>8);
  blk[5] = (int)((x0-x4)>>8);
  blk[6] = (int)((x3-x2)>>8);
  blk[7] = (int)((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(int *blk)
{
  long x0, x1, x2, x3, x4, x5, x6, x7, x8;

  /* intcut */
  if (!((x1 = (((long)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[((long)(blk[8*0]+32))>>6];
    return;
  }

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

  /* first stage */
  x8 = Q7*(x4+x5) + 4;
  x4 = (x8+(Q1-Q7)*x4)>>3;
  x5 = (x8-(Q1+Q7)*x5)>>3;
  x8 = Q3*(x6+x7) + 4;
  x6 = (x8-(Q3-Q5)*x6)>>3;
  x7 = (x8-(Q3+Q5)*x7)>>3;
  
  /* second stage */
  x8 = x0 + x1;
  x0 -= x1;
  x1 = Q6*(x3+x2) + 4;
  x2 = (x1-(Q2+Q6)*x2)>>3;
  x3 = (x1+(Q2-Q6)*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];
}

/* two dimensional inverse discrete cosine transform */
void idct(int *tt)
{
  int i;

  for (i=0; i<64; i+=8)
    idctrow(tt+i);

  for (i=0; i<8; i++)
    idctcol(tt+i);
}

void init_idct()
{
  int i;

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



/* private data */
 
/* cosine transform matrix for 8x1 IDCT */
static double c[8][8];
 
/* initialize DCT coefficient matrix */
 
void init_idctref()
{
  int freq, time;
  double scale;
 
  for (freq=0; freq < 8; freq++)
  {
    scale = (freq == 0) ? sqrt(0.125) : 0.5;
    for (time=0; time<8; time++)
      c[freq][time] = scale*cos((PI/8.0)*freq*(time + 0.5));
  }
}
 
/* perform IDCT matrix multiply for 8x8 coefficient block */
 

void idctref(int *block)
{
  int i, j;
  long k, v;
  double partial_product;
  double tmp[64];

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

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

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

  /* Transpose operation is integrated into address mapping by switching
     loop order of i and j */

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

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

      v = (long)floor(partial_product+0.5);
      block[8*i+j] = (int)((v<-256) ? -256 : ((v>255) ? 255 : v));
    }
}



void Quant(int *coeff, int *qcoeff, int QP, int Mode)
{
  int i,j;
  int level;
  int DQP = QP<<1, HQP = QP>>1;
  if (Mode == MODE_INTRA || Mode == MODE_INTRA_Q) 
  { // Intra 
	qcoeff[0] = mmax(1,mmin(254,coeff[0]>>3));

	for (i = 1; i < 64; i++) 
	{
		j = zig_zag[i];
		level = abs(coeff[j]) / DQP;
		qcoeff[i] =  mmin(127,mmax(-127,(coeff[j]>0?level:-level)));
    }
  }
  else 
  { // non Intra 
    for (i = 0; i < 64; i++) 
	{
		j = zig_zag[i];
		level = (abs(coeff[j])-HQP)  / DQP;
		qcoeff[i] = mmin(127,mmax(-127, (coeff[j]>0?level:-level)) );
    }
  }
}

void Dequant(int *qcoeff, int *rcoeff, int QP, int Mode)
{
  int i,j;
  int value;
  for (i = 0; i < 64; i++) 
  {
	j = zig_zag[i];
	value = abs(qcoeff[i]);
	if (value) 
	{
		if (QP % 2)
			rcoeff[j] = QP * ((value<<1) + 1);
		else
			rcoeff[j] = QP * ((value<<1) + 1) - 1;
		rcoeff[j] = qcoeff[i]>0?rcoeff[j]:(-rcoeff[j]);
    }
    else
		rcoeff[j] = 0;
  }
  if (Mode == MODE_INTRA || Mode == MODE_INTRA_Q) /* Intra */
      rcoeff[0] = qcoeff[0]<<3;
}