idct.c

linux下的一款播放器

/* ***** BEGIN LICENSE BLOCK *****
 * Source last modified: $Id: idct.c,v 1.3.40.1 2004/07/09 01:56:22 hubbe Exp $
 * 
 * Portions Copyright (c) 1995-2004 RealNetworks, Inc. All Rights Reserved.
 * 
 * The contents of this file, and the files included with this file,
 * are subject to the current version of the RealNetworks Public
 * Source License (the "RPSL") available at
 * http://www.helixcommunity.org/content/rpsl unless you have licensed
 * the file under the current version of the RealNetworks Community
 * Source License (the "RCSL") available at
 * http://www.helixcommunity.org/content/rcsl, in which case the RCSL
 * will apply. You may also obtain the license terms directly from
 * RealNetworks.  You may not use this file except in compliance with
 * the RPSL or, if you have a valid RCSL with RealNetworks applicable
 * to this file, the RCSL.  Please see the applicable RPSL or RCSL for
 * the rights, obligations and limitations governing use of the
 * contents of the file.
 * 
 * Alternatively, the contents of this file may be used under the
 * terms of the GNU General Public License Version 2 or later (the
 * "GPL") in which case the provisions of the GPL are applicable
 * instead of those above. If you wish to allow use of your version of
 * this file only under the terms of the GPL, and not to allow others
 * to use your version of this file under the terms of either the RPSL
 * or RCSL, indicate your decision by deleting the provisions above
 * and replace them with the notice and other provisions required by
 * the GPL. If you do not delete the provisions above, a recipient may
 * use your version of this file under the terms of any one of the
 * RPSL, the RCSL or the GPL.
 * 
 * This file is part of the Helix DNA Technology. RealNetworks is the
 * developer of the Original Code and owns the copyrights in the
 * portions it created.
 * 
 * This file, and the files included with this file, is distributed
 * and made available on an 'AS IS' basis, WITHOUT WARRANTY OF ANY
 * KIND, EITHER EXPRESS OR IMPLIED, AND REALNETWORKS HEREBY DISCLAIMS
 * ALL SUCH WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES
 * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, QUIET
 * ENJOYMENT OR NON-INFRINGEMENT.
 * 
 * Technology Compatibility Kit Test Suite(s) Location:
 *    http://www.helixcommunity.org/content/tck
 * 
 * Contributor(s):
 * 
 * ***** END LICENSE BLOCK ***** */

#include <math.h>
#include "dllindex.h"
#include "h261defs.h"
#include "h261func.h"
#include "clip.h"
#include "h263plus.h"

#ifdef _MACINTOSH
#include <string.h> // for memset
#endif

#ifndef UNREFERENCED_PARAMETER
#define UNREFERENCED_PARAMETER(x) (x) = (x)
#endif

extern S16 Recon [QUANT_MAX - QUANT_MIN + 1] [N_SYM_INDICES];

//  Function prototypes
static void init_inv_dct (void);
static void Fix_PTel_S4000_mismatch( void );
static void truncate_more( S32 * idct_tab_entry );
static void idct2_goertzel( SYMBOL sym[], int nsym, S32 x[8][4],
                            S16 recon[], int intra, int clean, int idct_class );
static int idct2_energy_test( SYMBOL sym[], int nsym,
                            S16 recon[] );                            
static void recon_intra_dc( U8 index, S32 vec[8]);
static void recon_dc( S32 y, S32 vec[8]);
static void recon_hor_ac( S32 y, S32 vec[8]);
static void recon_vert_ac( S32 y, S32 vec[8]);
static void update( S32 x[8], S32 index, S32 table[][8*8*8]);
static double dctfunc (int freq, int index);
static S32 dct_tab_entry (double x, double y);
static S32 combine (S32 iy, S32 ix);

extern void idct2_advanced_intra( SYMBOL sym[], int nsym, S32 x[8][4], S16 recon[],
                            U8 rDCpred, S8 rACpred[8], U8 rDCstore[1], S8 rACstore[8], 
                            U8 cDCpred, S8 cACpred[8], U8 cDCstore[1], S8 cACstore[8],
                            int predtype, int fixedDC, int leftBoundary, int upperBoundary);

#define PI              3.141592654
#define FRACBITS        6   /* Fractional bits in IDCT computation */
#define SCALE_FACTOR    64. /* 2**FRACBITS */
#define MAX_DCT_INDEX   1024
#define IDCT_NTAB1_BITS 5
#define IDCT_NTAB1_SIZE 32  /* 2**NTAB1_BITS entries in table for small values */
#define IDCT_NTAB2_SIZE ((MAX_DCT_INDEX - 1) / IDCT_NTAB1_SIZE)
                            /* Entries to handle values > NTAB1_SIZE */
#define PIXEL_MIN       0
#define PIXEL_MAX       255
#define CLIPMARGIN      300
#define CLIPMIN         (PIXEL_MIN - CLIPMARGIN)
#define CLIPMAX         (PIXEL_MAX + CLIPMARGIN)
#define N_DCT_INDEX     1024


//  Reconstruction levels for QUANT = 1,2,..,31
S16 Recon [QUANT_MAX - QUANT_MIN + 1] [N_SYM_INDICES];

//  Define zigzag scanning pattern
static int ZigZag[64] = {  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
};
static  int InvZZ[64];

static  S32    idct_tab [(IDCT_NTAB1_SIZE + IDCT_NTAB2_SIZE)][8*8][8];
                    /* LUT for all coeffs; [amplitude][zigzag position][] */
static  S32    dc_tab [MAX_DCT_INDEX];         /* LUT for DC coeff */
static  S32    hor_ac_tab [MAX_DCT_INDEX][2];  /* LUT for 1st hor AC coeff */
static  S32    vert_ac_tab [MAX_DCT_INDEX][4]; /* LUT for 1st vert AC coeff */
static  S32    intra_dc_tab[N_SYM_INDICES];    /* LUT for INTRA DC coeff */
PIXEL   clip[(CLIPMAX-CLIPMIN+1)];      /* LUT to limit to 0-255 */
static  int     even_odd_index[8*8];    /* Classify zigzag pos. as even or odd */


extern void InitReconTables( void )
{
    int level, index, quant;

//  QUANT=1 => 3,5,7,...
//  QUANT=2 => 5,9,13,...
//  QUANT=3 => 9,15,21,...
//  QUANT=4 => 11,19,27,...

    for (quant = QUANT_MIN; quant <= QUANT_MAX; quant++) {
        index = (quant + 1) / 2;
        for (level = 1; level < N_SYM_INDICES/2; level++) {
            index += quant;
            index = min( index, N_DCT_INDEX);
            Recon [quant - QUANT_MIN] [level] = index;
            Recon [quant - QUANT_MIN] [N_SYM_INDICES-level] = -index;
        }
        Recon [quant - QUANT_MIN] [0] = 0;
    }
    /* Generate zigzag table */
    for (index = 0; index < 64; index++) {
        InvZZ[ ZigZag[index] ] = index;
    }
    init_inv_dct();
    return;
}


//  Idct2 - Reconstruct DCT coeffs, perform IDCT, and clip to allowed pixel range */
//  Requires nsym > 0
extern void Idct2( SYMBOL sym[], int nsym, PIXEL x[], int xdim, S16 recon[], int clean)
{
    union {
        S16   bshort[8][8];
        S32    blong[8][4];
    } block;    /* Output from IDCT */
    int     intra;
    int     idct_class;

    idct_class = GENERAL;   /* assume the general case */

    /* look for situations involving specific involving 1,2,3 symbols */
    if (sym[0].type==0) 
        {
            switch (nsym) 
            {
                case 1:
                    idct_class = DC_ONLY;
                break;
                case 2:
                    if (sym[1].type == 0)  
                        idct_class = DC_AC_H;
                    if (sym[1].type == 1)  
                        idct_class = DC_AC_V;
                break;
                case 3:
                    if ( (sym[1].type | sym[2].type) == 0) 
                        idct_class = DC_3;
                break;
            } 
        }

    intra = YES;

    switch (idct_class)
            {
                case DC_ONLY:
                {   
                    idct2_goertzel( sym, nsym, block.blong, recon, intra, clean, DC_ONLY);
                    idct2_clip(x, xdim, block.blong, DC_ONLY);
                    break;
                }               
                case DC_AC_H:
                {   
                    idct2_goertzel( sym, nsym, block.blong, recon, intra, clean, DC_AC_H);
                    idct2_clip(x, xdim, block.blong, DC_AC_H);
                    break;
                }               
                case DC_AC_V:
                {   
                    idct2_goertzel( sym, nsym, block.blong, recon, intra, clean, DC_AC_V);
                    idct2_clip(x, xdim, block.blong, GENERAL);
                    break;
                }               
                case DC_3:
                {   
                    idct2_goertzel( sym, nsym, block.blong, recon, intra, clean, DC_3);
                    idct2_clip(x, xdim, block.blong, GENERAL);
                    break;
                }               
                default:                
                {   
                    idct2_goertzel( sym, nsym, block.blong, recon, intra, clean, GENERAL);
                    idct2_clip(x, xdim, block.blong, GENERAL);
                }               
            }               

    return;
}

//  Idct2Sum - Reconstruct DCT coeffs, perform IDCT, add to predition,
//  and clip to allowed pixel range.  Requires nsym > 0
extern void Idct2Sum( SYMBOL sym[], int nsym, PIXEL x[], int xdim, S16 recon[],int clean)
{
    union {
        S16   bshort[8][8];
        S32    blong[8][4];
    } block;    /* Output from IDCT */
    int     intra;
    int     idct_class;
    
    idct_class = GENERAL;   /* assume the general case */

    /* look for situations involving specific involving 1,2,3 symbols */
    if (sym[0].type==0) 
        {
            switch (nsym) 
            {
                case 1:
                    idct_class = DC_ONLY;
                break;
                case 2:
                    if (sym[1].type == 0)  
                        idct_class = DC_AC_H;
                    if (sym[1].type == 1)  
                        idct_class = DC_AC_V;
                break;
                case 3:
                    if ( (sym[1].type | sym[2].type) == 0) 
                        idct_class = DC_3;
                break;
            } 
        }

    intra = NO;
//    if (clean == YES || (idct2_energy_test( sym, nsym, recon) > 20)){ 
// activate this for sleazy IDCT in decoder
//    
    switch (idct_class)
            {
                case DC_ONLY:
                {   
                    idct2_goertzel( sym, nsym, block.blong, recon, intra, clean, DC_ONLY);
                    idct2sum_clip(x, xdim, block.blong, DC_ONLY);
                    break;
                }               
                case DC_AC_H:
                {   
                    idct2_goertzel( sym, nsym, block.blong, recon, intra, clean, DC_AC_H);
                    idct2sum_clip(x, xdim, block.blong, DC_AC_H);
                    break;
                }               
                case DC_AC_V:
                {   
                    idct2_goertzel( sym, nsym, block.blong, recon, intra, clean, DC_AC_V);
                    idct2sum_clip(x, xdim, block.blong, GENERAL);
                    break;
                }               
                case DC_3:
                {   
                    idct2_goertzel( sym, nsym, block.blong, recon, intra, clean, DC_3);
                    idct2sum_clip(x, xdim, block.blong, GENERAL);
                    break;
                }               
                default:                
                {   
                    idct2_goertzel( sym, nsym, block.blong, recon, intra, clean, GENERAL);
                    idct2sum_clip(x, xdim, block.blong, GENERAL);
                }               
            }
 //   }         // activate this for sleazy IDCT in decoder
    return;
}

//  Idct2_s16 - Reconstruct DCT coeffs, perform IDCT, and write as signed 16-bit values
//  Set output block to zero if nsym=0
extern void Idct2_s16( int intra, SYMBOL sym[], int nsym, S16 x[], int xdim, S16 recon[] )
{
    int i, clean = YES;
    union {
        S16   bshort[8][8];
        S32    blong[8][4];
    } block;    /* Output from IDCT */

    if (nsym == 0) {
        for (i = 0; i < 8; i++) {
            S32 * px = (S32 *)x;
            px[0] = px[1] = px[2] = px[3] = 0;
            x += xdim;
        }
        return;
    }

    idct2_goertzel( sym, nsym, block.blong, recon, intra, clean, GENERAL);

    // Shift out fractional bits
    for (i = 0; i < 8; i++) {
        x[i*xdim + 0] = block.bshort[i][0] >> FRACBITS;
        x[i*xdim + 1] = block.bshort[i][1] >> FRACBITS;
        x[i*xdim + 2] = block.bshort[i][2] >> FRACBITS;
        x[i*xdim + 3] = block.bshort[i][3] >> FRACBITS;
        x[i*xdim + 4] = block.bshort[i][7] >> FRACBITS;
        x[i*xdim + 5] = block.bshort[i][6] >> FRACBITS;
        x[i*xdim + 6] = block.bshort[i][5] >> FRACBITS;
        x[i*xdim + 7] = block.bshort[i][4] >> FRACBITS;
    }
}


//  Initialize tables for inverse DCT
//  Note: This routine has not been optimized for speed
static void init_inv_dct (void)
{
    int             i,j,m,n, index, zzpos;
    double          magn;               /* amplitude of DCT coefficient */
    static  double  bfunc[8][8][4][4];  /* DCT basis functions [vert freq][hor freq][][] */

    for (n=0; n < 8; n++) {     /* Construct 2-D basis functions */
      for (m=0; m < 8; m++) {
        for (j=0; j < 4; j++) {
          for (i=0; i < 4; i++) {
              bfunc[n][m][j][i] = SCALE_FACTOR * dctfunc(n,j) * dctfunc(m,i);
          }
        }
      }
    }
    //  Initialize table for INTRA DC coeff reconstruction */
    for (index = 0; index < N_SYM_INDICES; index++) {
        magn = 8 * index;
        //printf( "Init index = %d   magn = %f \n", index, magn);
        intra_dc_tab [index] = dct_tab_entry (magn * bfunc[0][0][0][0],
                                              magn * bfunc[0][0][0][1]);
    }
    //  128 is represented by index=255
    index = 255;
    magn = 8 * 128;
    //printf( "Init index = %d   magn = %f \n", index, magn);
    intra_dc_tab [index] = dct_tab_entry (magn * bfunc[0][0][0][0],
                                              magn * bfunc[0][0][0][1]);
    
    //  Initialize tables for DC and first two AC coeffs
    for (index = 0; index < MAX_DCT_INDEX; index++) {
        magn = 2*index + 1;
        //printf( "Init index = %d   magn = %f \n", index, magn);
        dc_tab [index] = dct_tab_entry (magn * bfunc[0][0][0][0],
                                        magn * bfunc[0][0][0][1]);
        for (i = 0; i < 2; i++) {
            hor_ac_tab [index][i]
                        = dct_tab_entry (magn * bfunc[0][1][0][2*i],
                                         magn * bfunc[0][1][0][2*i+1]);
        }
        for (j = 0; j < 4; j++) {
            vert_ac_tab [index][j]
                        = dct_tab_entry (magn * bfunc[1][0][j][0],
                                         magn * bfunc[1][0][j][1]);
        }
    }
    //  Initialize table for all coeffs
    for (index = 0; index < IDCT_NTAB1_SIZE; index++) {
        magn = 2*index + 1;
        //printf( "Init index = %d   magn = %f \n", index, magn);
        for (zzpos = 0; zzpos < 8*8; zzpos++) {
            n = InvZZ[zzpos] / 8;
            m = InvZZ[zzpos] % 8;
            for (j=0; j < 4; j++) {
              for (i=0; i < 2; i++) {
                  idct_tab [index][zzpos][2*j+i]
                    = dct_tab_entry (magn * bfunc[n][m][j][2*i],
                                    magn * bfunc[n][m][j][2*i+1]);