dct.c

H.263标准的c语言实现

  return 0;
  #endif
}

/**********************************************************************
 *
 *	Name:           Scan and Descan
 *	Description:	Does scanning of the an 8x8 block quantized dct
 *                      coefficients.
 *                      (zigzag, or alernate horizontal, alternate
 *                      vertical when advanced intra coding is used)
 *
 *	Input:          64 quantized coefficients in a 1D array
 *	Returns:	0
 *	Side effects:	scans the dct qcoeff in the desired order
 *                      (zigzag if advanced intra coding is not used)
 *
 *	Date: 970716
 *      Author: Guy Cote <[e-mail]>
 *
 **********************************************************************/

void Scan (int *qcoeff, int Intra_Mode)
{
  int i, j, k;
  int store_qcoeff[384];
  /* Zigzag - scanning of a Macroblock */
  /* Choose Scanning if Advanced Intra Coding Mode is used */

  for (i = 0; i < 384; i++)
  {
    store_qcoeff[i] = *(qcoeff + i);
  }

  for (k = 0; k < 6; k++)
  {                             /* Do for all blocks */
    for (i = 0; i < 8; i++)
    {
      for (j = 0; j < 8; j++)
      {
        if (advanced_intra_coding)
        {
          switch (Intra_Mode)
          {
            case INTRA_MODE_DC:
              *(qcoeff + zigzag[i][j] + k * 64) = (int) (store_qcoeff[i * 8 + j + k * 64]);
              break;
            case INTRA_MODE_VERT_AC:
              *(qcoeff + alternate_horizontal[i][j] + k * 64) = (int) (store_qcoeff[i * 8 + j + k * 64]);
              break;
            case INTRA_MODE_HORI_AC:
              *(qcoeff + alternate_vertical[i][j] + k * 64) = (int) (store_qcoeff[i * 8 + j + k * 64]);
              break;
            default:
              printf ("Invalid Intra_Mode in Advanced Intra Coding");
              exit (-1);
              break;
          }
        } else
        {
          *(qcoeff + zigzag[i][j] + k * 64) = (int) (store_qcoeff[i * 8 + j + k * 64]);
        }
      }
    }
  }
  return;
}
void DeScan (int *qcoeff, int Intra_Mode)
{
  /* Descan coefficients first */
  /* Choose Scanning if Advanced Intra Coding Mode is used */
  int store_qcoeff[384];
  int i, j, k;

  for (i = 0; i < 384; i++)
  {
    store_qcoeff[i] = *(qcoeff + i);
  }
  for (k = 0; k < 6; k++)
  {                             /* Do for all blocks */
    for (i = 0; i < 8; i++)
    {
      for (j = 0; j < 8; j++)
      {
        if (advanced_intra_coding)
        {                       /* Descan with proper scanning matrix */
          switch (Intra_Mode)
          {
            case INTRA_MODE_DC:
              *(qcoeff + k * 64 + i * 8 + j) = *(store_qcoeff + k * 64 + zigzag[i][j]);
              break;
            case INTRA_MODE_VERT_AC:
              *(qcoeff + k * 64 + i * 8 + j) = *(store_qcoeff + k * 64 + alternate_horizontal[i][j]);
              break;
            case INTRA_MODE_HORI_AC:
              *(qcoeff + k * 64 + i * 8 + j) = *(store_qcoeff + k * 64 + alternate_vertical[i][j]);
              break;
            default:
              printf ("Invalid Intra_Mode in Advanced Intra Coding\n");
              exit (-1);
              break;
          }
        } else
        {
          *(qcoeff + k * 64 + i * 8 + j) = *(store_qcoeff + k * 64 + zigzag[i][j]);
        }
      }
    }
  }
  return;
}

#ifdef FASTIDCT

/**********************************************************************
 *
 *	Name:        idct
 *	Description:    inverse dct on 64 coefficients
 *                      single precision floats
 *
 *	Input:        64 coefficients, block for 64 pixels
 *	Returns:        0
 *	Side effects:
 *
 *	Date: 930128	Author: [e-mail]
 *
 **********************************************************************/

int idct (int *coeff, int *block)
{
  int j1, i, j;
  double b[8], b1[8], d[8][8];
  double f0 = .7071068;
  double f1 = .4903926;
  double f2 = .4619398;
  double f3 = .4157348;
  double f4 = .3535534;
  double f5 = .2777851;
  double f6 = .1913417;
  double f7 = .0975452;
  double e, f, g, h;

  /* Horizontal */

  for (i = 0; i < 8; i++)
  {

    for (j = 0; j < 8; j++)
      b[j] = coeff[j + i * 8];

    e = b[1] * f7 - b[7] * f1;
    h = b[7] * f7 + b[1] * f1;
    f = b[5] * f3 - b[3] * f5;
    g = b[3] * f3 + b[5] * f5;

    b1[0] = (b[0] + b[4]) * f4;
    b1[1] = (b[0] - b[4]) * f4;
    b1[2] = b[2] * f6 - b[6] * f2;
    b1[3] = b[6] * f6 + b[2] * f2;
    b[4] = e + f;
    b1[5] = e - f;
    b1[6] = h - g;
    b[7] = h + g;

    b[5] = (b1[6] - b1[5]) * f0;
    b[6] = (b1[6] + b1[5]) * f0;
    b[0] = b1[0] + b1[3];
    b[1] = b1[1] + b1[2];
    b[2] = b1[1] - b1[2];
    b[3] = b1[0] - b1[3];

    for (j = 0; j < 4; j++)
    {
      j1 = 7 - j;
      d[i][j] = b[j] + b[j1];
      d[i][j1] = b[j] - b[j1];
    }
  }


  /* Vertical */

  for (i = 0; i < 8; i++)
  {
    for (j = 0; j < 8; j++)
    {
      b[j] = d[j][i];
    }
    e = b[1] * f7 - b[7] * f1;
    h = b[7] * f7 + b[1] * f1;
    f = b[5] * f3 - b[3] * f5;
    g = b[3] * f3 + b[5] * f5;

    b1[0] = (b[0] + b[4]) * f4;
    b1[1] = (b[0] - b[4]) * f4;
    b1[2] = b[2] * f6 - b[6] * f2;
    b1[3] = b[6] * f6 + b[2] * f2;
    b[4] = e + f;
    b1[5] = e - f;
    b1[6] = h - g;
    b[7] = h + g;

    b[5] = (b1[6] - b1[5]) * f0;
    b[6] = (b1[6] + b1[5]) * f0;
    b[0] = b1[0] + b1[3];
    b[1] = b1[1] + b1[2];
    b[2] = b1[1] - b1[2];
    b[3] = b1[0] - b1[3];

    for (j = 0; j < 4; j++)
    {
      j1 = 7 - j;
      d[j][i] = b[j] + b[j1];
      d[j1][i] = b[j] - b[j1];
    }
  }

  for (i = 0; i < 8; i++)
  {
    for (j = 0; j < 8; j++)
    {
      *(block + i * 8 + j) = mnint (d[i][j]);
    }
  }
  return 0;
}

#else
/* Perform IEEE 1180 reference (64-bit floating point, separable 8x1
 * direct matrix multiply) Inverse Discrete Cosine Transform */


/* Here we use math.h to generate constants.  Compiler results may vary a
 * little */


/* 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 *coeff, int *block)
{
#if 1
      static int *blk;
	int i;
	static long X0, X1, X2, X3, X4, X5, X6, X7, X8;

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

	for (i = 0; i < 8; i++)		/* idct rows */
	{
		blk = coeff + (i << 3);
		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;
			continue;
		}

		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] = (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);

	}							/* end for ( i = 0; i < 8; ++i ) IDCT-rows */



	for (i = 0; i < 8; i++)		/* idct columns */
	{
		blk = coeff + i;
		/* 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];
			continue;
		}

		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];
	}
    for(i = 0; i<64;i++)
        block[i] = coeff[i];
    return ;
#else
  int i, j, k, v;
  double partial_product;
  double tmp[64];
  int tmp2[64];

  for (i = 0; i < 8; i++)
    for (j = 0; j < 8; j++)
      tmp2[j + i * 8] = coeff[j + i * 8];

  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] * tmp2[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 = (int) floor (partial_product + 0.5);
      block[8 * i + j] = (v < -256) ? -256 : ((v > 255) ? 255 : v);
    }
#endif

}
#endif