dct.cpp

h.263一款不错的网络视频监控系统软件

////////////////////////////////////////////////////////////////////////////
//
//
//    Project     : VideoNet version 1.1.
//    Description : Peer to Peer Video Conferencing over the LAN.
//	  Author      :	Nagareshwar Y Talekar ( [e-mail])
//    Date        : 15-6-2004.
//
//    I have converted origional fast h.263 encoder library from C to C++ 
//	  so that it can be integrated into any windows application easily.
//	  I have removed some of unnecessary codes/files from the
//	  fast h263 library.Also moved definitions and declarations
//	  in their proper .h and .cpp files.
//
//    File description : 
//    Name    : dct.cpp
//
//
/////////////////////////////////////////////////////////////////////////////

/*************************************************
 * libr263: fast H.263 encoder library
 *
 * Copyright (C) 1996, Roalt Aalmoes, Twente University
 * SPA multimedia group
 *
 * Based on Telenor TMN 1.6 encoder (Copyright (C) 1995, Telenor R&D)
 * created by Karl Lillevold 
 *
 * Author encoder: Roalt Aalmoes, <[e-mail]>
 * 
 * Date: 31-07-96
 **************************************************/

/*****************************************************************
 * Some routines are translated from Gisle Bj鴑tegaards's FORTRAN
 * routines by [e-mail]
 *
 *****************************************************************/

/*****************************************************************
 * This source includes sources from Berkeley's MPEG-1 encoder
 * which are copyright of Berkeley University, California, USA
 *****************************************************************/


#include "stdafx.h"
#include "dct.h"


#include <math.h>

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

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

#ifndef FASTDCT

/**********************************************************************
 *
 *	Name:		Dct
 *	Description:	Does dct on an 8x8 block, does zigzag-scanning of
 *			coefficients
 *
 *	Input:		64 pixels in a 1D array
 *	Returns:	64 coefficients in a 1D array
 *	Side effects:	
 *
 *	Date: 930128	Author: [e-mail]
 *
 **********************************************************************/

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

  for (i = 0, k = 0; i < 8; i++, k += 8) {
    for (j = 0; j < 8; j++) {
      b[j] = block[k+j];
    }
    /* Horizontal transform */
    for (j = 0; j < 4; j++) {
      j1 = 7 - j;
      b1[j] = b[j] + b[j1];
      b1[j1] = b[j] - b[j1];
    }
    b[0] = b1[0] + b1[3];
    b[1] = b1[1] + b1[2];
    b[2] = b1[1] - b1[2];
    b[3] = b1[0] - b1[3];
    b[4] = b1[4];
    b[5] = (b1[6] - b1[5]) * f0;
    b[6] = (b1[6] + b1[5]) * f0;
    b[7] = b1[7];
    d[i][0] = (b[0] + b[1]) * f4;
    d[i][4] = (b[0] - b[1]) * f4;
    d[i][2] = b[2] * f6 + b[3] * f2;
    d[i][6] = b[3] * f6 - b[2] * f2;
    b1[4] = b[4] + b[5];
    b1[7] = b[7] + b[6];
    b1[5] = b[4] - b[5];
    b1[6] = b[7] - b[6];
    d[i][1] = b1[4] * f7 + b1[7] * f1;
    d[i][5] = b1[5] * f3 + b1[6] * f5;
    d[i][7] = b1[7] * f7 - b1[4] * f1;
    d[i][3] = b1[6] * f3 - b1[5] * f5;
  }
  /* Vertical transform */
  for (i = 0; i < 8; i++) {
    for (j = 0; j < 4; j++) {
      j1 = 7 - j;
      b1[j] = d[j][i] + d[j1][i];
      b1[j1] = d[j][i] - d[j1][i];
    }
    b[0] = b1[0] + b1[3];
    b[1] = b1[1] + b1[2];
    b[2] = b1[1] - b1[2];
    b[3] = b1[0] - b1[3];
    b[4] = b1[4];
    b[5] = (b1[6] - b1[5]) * f0;
    b[6] = (b1[6] + b1[5]) * f0;
    b[7] = b1[7];
    d[0][i] = (b[0] + b[1]) * f4;
    d[4][i] = (b[0] - b[1]) * f4;
    d[2][i] = b[2] * f6 + b[3] * f2;
    d[6][i] = b[3] * f6 - b[2] * f2;
    b1[4] = b[4] + b[5];
    b1[7] = b[7] + b[6];
    b1[5] = b[4] - b[5];
    b1[6] = b[7] - b[6];
    d[1][i] = b1[4] * f7 + b1[7] * f1;
    d[5][i] = b1[5] * f3 + b1[6] * f5;
    d[7][i] = b1[7] * f7 - b1[4] * f1;
    d[3][i] = b1[6] * f3 - b1[5] * f5;
  }
  /* Zigzag - scanning */
  for (i = 0; i < 8; i++) {
    for (j = 0; j < 8; j++) {
      *(coeff + zigzag[i][j]) = (int)(d[i][j]);
    }
  }
  return 0;
}

#else
				/* Start of MPEG DCT operation */

typedef unsigned char uint8;
typedef char int8;
typedef unsigned short uint16;
typedef short int16;

#ifdef LONG_32		
typedef unsigned long uint32;
typedef long int32;
#else
typedef unsigned int uint32;
typedef int int32;
#endif

/* this is ansi.h */


#undef _ANSI_ARGS_
#undef const

#ifdef NON_ANSI_COMPILER
#define _ANSI_ARGS_(x)       ()
#define CONST

#else

#define _ANSI_ARGS_(x)   x
#define CONST const

#ifdef __cplusplus
#define VARARGS (...)
#else
#define VARARGS ()
#endif

#endif


#define DCTSIZE     8   /* you really don't want to change this */
#define DCTSIZE_SQ 64   /* you really don't want to change this */

#define DCTSIZE2    DCTSIZE*DCTSIZE
typedef short DCTELEM;
typedef DCTELEM DCTBLOCK[DCTSIZE2];
typedef DCTELEM DCTBLOCK_2D[DCTSIZE][DCTSIZE];

/* 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 proper signed right shift of an int32 quantity.
 * It is only applied with constant shift counts.  SHIFT_TEMPS must be
 * included in the variables of any routine using RIGHT_SHIFT.
 */

#ifdef RIGHT_SHIFT_IS_UNSIGNED
#define SHIFT_TEMPS     int32 shift_temp

/*
#define RIGHT_SHIFT(x,shft)  ((shift_temp = (x)) < 0 ? (shift_temp >> (shft)) |((~((int32) 0)) << (32-(shft))) : (shift_temp >> (shft)))
     
	 */   

#else
#define SHIFT_TEMPS
#define RIGHT_SHIFT(x,shft)     ((x) >> (shft))
#endif
          
#define LG2_DCT_SCALE 16

#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
     
/*==================*
 * TYPE DEFINITIONS *
 *==================*/

/*  
 *  your basic Block type
 */
typedef int32 Block[DCTSIZE][DCTSIZE];
typedef int32 FlatBlock[DCTSIZE_SQ];
typedef	    int32   LumBlock[2*DCTSIZE][2*DCTSIZE];
typedef	    int32   ChromBlock[DCTSIZE][DCTSIZE];



/* Prototypes */
void reference_fwd_dct(Block block, Block dest);

void mp_fwd_dct_fast(Block data2d, Block dest2d);

/*void mp_fwd_dct_fast _ANSI_ARGS_((int16 *data2d, int16 *dest2d));
 */
void init_fdct(void);

/*
 * --------------------------------------------------------------
 *
 * 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
 *
 * --------------------------------------------------------------
 */
/*
__inline__ voidmp_fwd_dct_fast(data2d, dest2d)
    Block data2d, dest2d;
    */
__inline__ void mp_fwd_dct_fast(Block data2d, Block dest2d)
{
    int32 *data = (int32 *) data2d;	/* this algorithm wants
					 * a 1-d array */
    int32 *dest = (int32 *) dest2d;
    int rowctr, columncounter;
    register int32 *inptr, *outptr;
    int32 workspace[DCTSIZE_SQ];
    int32 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
    int32 tmp10, tmp11, tmp12, tmp13;
    int32 tmp14, tmp15, tmp16, tmp17;
    int32 tmp25, tmp26;

    SHIFT_TEMPS;

    /*
     * 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;
   
    /* PASS ONE */
	
    for (rowctr = DCTSIZE - 1; rowctr >= 0; rowctr--) {
      /*
       * many tmps have nonoverlapping lifetime -- flashy
       * register colourers should be able to do this lot
       * very well
       */
      /* 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] = (int32) UNFIXH((tmp10 + tmp11) * SIN_1_4);
      outptr[DCTSIZE * 4] = (int32) UNFIXH((tmp10 - tmp11) * COS_1_4);
      
      outptr[DCTSIZE * 2] = (int32) UNFIXH(tmp13 * COS_1_8 + tmp12 * SIN_1_8);
      outptr[DCTSIZE * 6] = (int32) 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] = (int32) UNFIXH(tmp17 * OCOS_1_16 + tmp14 * OSIN_1_16);
      outptr[DCTSIZE * 7] = (int32) UNFIXH(tmp17 * OCOS_7_16 - tmp14 * OSIN_7_16);
      outptr[DCTSIZE * 5] = (int32) UNFIXH(tmp26 * OCOS_5_16 + tmp25 * OSIN_5_16);
      outptr[DCTSIZE * 3] = (int32) 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;

    columncounter = 0;
    /* PASS TWO */

    for (rowctr = DCTSIZE - 1; rowctr >= 0; rowctr--) {
      /*
       * many tmps have nonoverlapping lifetime -- flashy
       * register colourers should be able to do this lot
       * very well
       */
      /* SHIFT_TEMPS */