idctmm32.cpp

来自「这是一组DCT和iDCT的代码」· C++ 代码 · 共 1,707 行 · 第 1/4 页

#include "StdAfx.h"
#pragma warning(once:4305 4244)

//////////////////////////////////////////////////////////////////////
//  not functional
//////////////////////////////////////////////////////////////////////
// MPEG2AVI
// -------- 
//  v0.16B34
//    performance optimization, idct_mmx32_rows() and idct_mmx32_cols() now
//    transpose their respective outputs "in-place" (saves a bit of time)
//
//  v0.16B33 initial release
//
// MMX32 iDCT algorithm  (IEEE-1180 compliant) :: idct_mmx32()
//
//  This IDCT implementation is based on Intel Application Note AP-922.
//
//  This file implements the idct algorithm with no transpose.
//  The other file (idctmm32_transpose.c) is faster, but transposes the 
//  output-matrix.  (Intel's code-listing produces a transposed output.)
//
//  ALGORITHM OVERVIEW
//  ------------------
// This was one of the harder pieces of work to code.
// Intel's app-note focuses on the numerical theory/issues of the IDCT 
// implementation, but assumes the programmer is familiar with the
// requisite mathematics, leaving the exact form of the complete IDCT 
// code-listing up to the programmer's imagination.
//
// I played around with Intel's code-fragments for quite a few hours.  
// This file is *A* working IDCT implementation, but it may not be
// the implementation Intel originally intended.  Rest assured, I've
// done everything in my power to guarantee its correctness. 
// This implementation passes all six IEEE accuracy tests by a fair margin.
//
//   My IDCT algorithm consists of 4 steps:
//
//   1) IDCT-row transformation (using the IDCT-row function) on all 8 rows
//      This yields an intermediate 8x8 matrix.
//
//   2) transpose of intermediate matrix (mandatory)
//
//   3) IDCT-row transformation (2nd time) on all 8 rows of the 
//      intermediate matrix.
//      At this point, we have the final-result, in transposed form.
//
//   4) post-transformation matrix transpose 
//      (not necessary if the input-data is already transposed, this could
//       be done during the MPEG "zig-zag" scan, but since my algorithm
//       requires at least one transpose operation, why not re-use the
//       transpose-code.)
//
//   Although the (1st) and (3rd) steps use the same basic row-transform 
//   operation, the (3rd) step uses different shift&round constants 
//   (explained later.)
//
//   Also note that the intermediate transpose (2) would not be neccessary,
//   if the subsequent operation were a iDCT-column transformation.  Since
//   we only have the iDCT-row transform, we transpose the intermediate
//   matrix and use the iDCT-row transform a 2nd time.  I suppose one
//   a faster (but more complicated) code-implementation is possible, 
//   if these steps were merged.
//
//   I had to change some constants/variables for my method to work :
//
//      As given by Intel, #SHIFT_INV_COL and #RND_INV_COL are wrong.  
//      Not surprising since I'm probably implementing the IDCT in
//      perverse fashion.  
//      round_inv_col[], which is given as "4 short" values, should have the
//      same dimensions as round_inv_row[].  The corrected variables are 
//      shown.
//
//      Intel's code defines 4 tables of constants.  My code only uses only
//      one of these tables, row#0.  
//
//   IMPLEMENTATION DETAILs
//   ----------------------
// 
//   I divided the 4-steps of my algorithm into two subroutines,
//    1) idct_mmx32_rows() - transforms 8 rows, then transpose
//    2) idct_mmx32_cols() - transforms 8 rows, then transpose
//       yields final result ("drop-in" direct replacement for INT32 IDCT)
//
//   idct_mmx32_cols() is a carbon-copy of idct_mmx32_rows(), i.e. both
//   execute a row-by-row transformations.  Only the shift&rounding 
//   coefficients differ.
//
//      In the 1st function (rows), the shift & round instructions use 
//       SHIFT_INV_ROW & round_inv_row[] (renamed to r_inv_row[])
//
//      In the 2nd function (cols)-> r_inv_col[], and
//       SHIFT_INV_COL & round_inv_col[] (renamed to r_inv_col[])
//
//   Each function contains an integrated transpose-operator, which comes
//   AFTER the primary transformation operation.  In the future, I'll optimize
//   the code to do more of the transpose-work "in-place".  Right now, I've
//   left the code as two subroutines and a main calling function, so other
//   people can read the code more easily.


//;=============================================================================
//;
//;  AP-922   http://developer.intel.com/vtune/cbts/strmsimd
//; These examples contain code fragments for first stage iDCT 8x8
//; (for rows) and first stage DCT 8x8 (for columns)
//;
//;=============================================================================
/*
mword typedef qword
qword ptr equ mword ptr */

#define BITS_INV_ACC  4 //; 4 or 5 for IEEE
  // 5 yields higher accuracy, but lessens dynamic range on the input matrix
#define SHIFT_INV_ROW (16 - BITS_INV_ACC)
// for goofy mmx implementation #define SHIFT_INV_COL (1 + BITS_INV_ACC +14 )  // changed from Intel's val)
#define SHIFT_INV_COL (1 + BITS_INV_ACC )

#define RND_INV_ROW   (1 << (SHIFT_INV_ROW-1))
#define RND_INV_COL   (1 << (SHIFT_INV_COL-1)) 
#define RND_INV_CORR  (RND_INV_COL - 1)   //; correction -1.0 and round
//#define RND_INV_ROW   (1024 * (6 - BITS_INV_ACC)) //; 1 << (SHIFT_INV_ROW-1)
//#define RND_INV_COL   (16 * (BITS_INV_ACC - 3)) //; 1 << (SHIFT_INV_COL-1)


//.data
//Align 16
const static long r_inv_row[2] = { RND_INV_ROW, RND_INV_ROW};
const static long r_inv_col[2] = {RND_INV_COL, RND_INV_COL};
const static long r_inv_corr[2] = {RND_INV_CORR, RND_INV_CORR };

//const static short r_inv_col[4] = 
//  {RND_INV_COL, RND_INV_COL, RND_INV_COL, RND_INV_COL};
//const static short r_inv_corr[4] =
//  {RND_INV_CORR, RND_INV_CORR, RND_INV_CORR, RND_INV_CORR};

/* constants for the forward DCT

//#define BITS_FRW_ACC  3 //; 2 or 3 for accuracy
//#define SHIFT_FRW_COL BITS_FRW_ACC
//#define SHIFT_FRW_ROW (BITS_FRW_ACC + 17)
//#define RND_FRW_ROW   (262144 * (BITS_FRW_ACC - 1)) //; 1 << (SHIFT_FRW_ROW-1)

const static __int64 one_corr = 0x0001000100010001;
const static long r_frw_row[2] = {RND_FRW_ROW, RND_FRW_ROW };

//const static short tg_1_16[4] = {13036, 13036, 13036, 13036 }; //tg * (2<<16) + 0.5
//const static short tg_2_16[4] = {27146, 27146, 27146, 27146 }; //tg * (2<<16) + 0.5
//const static short tg_3_16[4] = {-21746, -21746, -21746, -21746 }; //tg * (2<<16) + 0.5
//const static short cos_4_16[4] = {-19195, -19195, -19195, -19195 }; //cos * (2<<16) + 0.5
//const static short ocos_4_16[4] = {23170, 23170, 23170, 23170 }; //cos * (2<<15) + 0.5

//concatenated table, for forward DCT transformation
const static short tg_all_16[] = {
  13036, 13036, 13036, 13036,   // tg * (2<<16) + 0.5
  27146, 27146, 27146, 27146,   //tg * (2<<16) + 0.5
  -21746, -21746, -21746, -21746, // tg * (2<<16) + 0.5
  -19195, -19195, -19195, -19195, //cos * (2<<16) + 0.5
  23170, 23170, 23170, 23170 }; //cos * (2<<15) + 0.5

#define tg_1_16 (tg_all_16 + 0)
#define tg_2_16 (tg_all_16 + 8)
#define tg_3_16 (tg_all_16 + 16)
#define cos_4_16 (tg_all_16 + 24)
#define ocos_4_16 (tg_all_16 + 32)
*/
/*
;=============================================================================
;
; The first stage iDCT 8x8 - inverse DCTs of rows
;
;-----------------------------------------------------------------------------
; The 8-point inverse DCT direct algorithm
;-----------------------------------------------------------------------------
;
; static const short w[32] = {
; FIX(cos_4_16), FIX(cos_2_16), FIX(cos_4_16), FIX(cos_6_16),
; FIX(cos_4_16), FIX(cos_6_16), -FIX(cos_4_16), -FIX(cos_2_16),
; FIX(cos_4_16), -FIX(cos_6_16), -FIX(cos_4_16), FIX(cos_2_16),
; FIX(cos_4_16), -FIX(cos_2_16), FIX(cos_4_16), -FIX(cos_6_16),
; FIX(cos_1_16), FIX(cos_3_16), FIX(cos_5_16), FIX(cos_7_16),
; FIX(cos_3_16), -FIX(cos_7_16), -FIX(cos_1_16), -FIX(cos_5_16),
; FIX(cos_5_16), -FIX(cos_1_16), FIX(cos_7_16), FIX(cos_3_16),
; FIX(cos_7_16), -FIX(cos_5_16), FIX(cos_3_16), -FIX(cos_1_16) };
;
; #define DCT_8_INV_ROW(x, y)

;{
; int a0, a1, a2, a3, b0, b1, b2, b3;
;
; a0 =x[0]*w[0]+x[2]*w[1]+x[4]*w[2]+x[6]*w[3];
; a1 =x[0]*w[4]+x[2]*w[5]+x[4]*w[6]+x[6]*w[7];
; a2 = x[0] * w[ 8] + x[2] * w[ 9] + x[4] * w[10] + x[6] * w[11];
; a3 = x[0] * w[12] + x[2] * w[13] + x[4] * w[14] + x[6] * w[15];
; b0 = x[1] * w[16] + x[3] * w[17] + x[5] * w[18] + x[7] * w[19];
; b1 = x[1] * w[20] + x[3] * w[21] + x[5] * w[22] + x[7] * w[23];
; b2 = x[1] * w[24] + x[3] * w[25] + x[5] * w[26] + x[7] * w[27];
; b3 = x[1] * w[28] + x[3] * w[29] + x[5] * w[30] + x[7] * w[31];
;
; y[0] = SHIFT_ROUND ( a0 + b0 );
; y[1] = SHIFT_ROUND ( a1 + b1 );
; y[2] = SHIFT_ROUND ( a2 + b2 );
; y[3] = SHIFT_ROUND ( a3 + b3 );
; y[4] = SHIFT_ROUND ( a3 - b3 );
; y[5] = SHIFT_ROUND ( a2 - b2 );
; y[6] = SHIFT_ROUND ( a1 - b1 );
; y[7] = SHIFT_ROUND ( a0 - b0 );
;}
;
;-----------------------------------------------------------------------------
;
; In this implementation the outputs of the iDCT-1D are multiplied
; for rows 0,4 - by cos_4_16,
; for rows 1,7 - by cos_1_16,
; for rows 2,6 - by cos_2_16,
; for rows 3,5 - by cos_3_16
; and are shifted to the left for better accuracy
;
; For the constants used,
; FIX(float_const) = (short) (float_const * (1<<15) + 0.5)
;
;=============================================================================
;=============================================================================
IF _MMX ; MMX code
;=============================================================================
/*

//; Table for rows 0,4 - constants are multiplied by cos_4_16
const short tab_i_04[] = {
  16384, 16384, 16384, -16384,  // ; movq-> w06 w04 w02 w00
  21407, 8867, 8867, -21407,    // w07 w05 w03 w01
  16384, -16384, 16384, 16384,  //; w14 w12 w10 w08
  -8867, 21407, -21407, -8867,  //; w15 w13 w11 w09
  22725, 12873, 19266, -22725,  //; w22 w20 w18 w16
  19266, 4520, -4520, -12873,   //; w23 w21 w19 w17
  12873, 4520, 4520, 19266,   //; w30 w28 w26 w24
  -22725, 19266, -12873, -22725 };//w31 w29 w27 w25

//; Table for rows 1,7 - constants are multiplied by cos_1_16
const short tab_i_17[] = {
  22725, 22725, 22725, -22725,  // ; movq-> w06 w04 w02 w00
  29692, 12299, 12299, -29692,  //  ; w07 w05 w03 w01
  22725, -22725, 22725, 22725,  //; w14 w12 w10 w08
  -12299, 29692, -29692, -12299,  //; w15 w13 w11 w09
  31521, 17855, 26722, -31521,  //; w22 w20 w18 w16
  26722, 6270, -6270, -17855,   //; w23 w21 w19 w17
  17855, 6270, 6270, 26722,   //; w30 w28 w26 w24
  -31521, 26722, -17855, -31521}; // w31 w29 w27 w25

//; Table for rows 2,6 - constants are multiplied by cos_2_16
const short tab_i_26[] = {
  21407, 21407, 21407, -21407,  // ; movq-> w06 w04 w02 w00
  27969, 11585, 11585, -27969,  // ; w07 w05 w03 w01
  21407, -21407, 21407, 21407,  // ; w14 w12 w10 w08
  -11585, 27969, -27969, -11585,  //  ;w15 w13 w11 w09
  29692, 16819, 25172, -29692,  // ;w22 w20 w18 w16
  25172, 5906, -5906, -16819,   // ;w23 w21 w19 w17
  16819, 5906, 5906, 25172,     // ;w30 w28 w26 w24
  -29692, 25172, -16819, -29692}; //  ;w31 w29 w27 w25


//; Table for rows 3,5 - constants are multiplied by cos_3_16
const short tab_i_35[] = {
  19266, 19266, 19266, -19266,  //; movq-> w06 w04 w02 w00
  25172, 10426, 10426, -25172,  //; w07 w05 w03 w01
  19266, -19266, 19266, 19266,  //; w14 w12 w10 w08
  -10426, 25172, -25172, -10426,  //; w15 w13 w11 w09
  26722, 15137, 22654, -26722,  //; w22 w20 w18 w16
  22654, 5315, -5315, -15137,   //; w23 w21 w19 w17
  15137, 5315, 5315, 22654,   //; w30 w28 w26 w24
  -26722, 22654, -15137, -26722}; //; w31 w29 w27 w25
*/

// CONCATENATED TABLE, rows 0,1,2,3,4,5,6,7 (in order )
//
// In our implementation, however, we only use row0 !
//
static const short tab_i_01234567[] = {
  //row0, this row is required
  16384, 16384, 16384, -16384,  // ; movq-> w06 w04 w02 w00
  21407, 8867, 8867, -21407,    // w07 w05 w03 w01
  16384, -16384, 16384, 16384,  //; w14 w12 w10 w08
  -8867, 21407, -21407, -8867,  //; w15 w13 w11 w09
  22725, 12873, 19266, -22725,  //; w22 w20 w18 w16
  19266, 4520, -4520, -12873,   //; w23 w21 w19 w17
  12873, 4520, 4520, 19266,   //; w30 w28 w26 w24
  -22725, 19266, -12873, -22725,  //w31 w29 w27 w25

  // the rest of these rows (1-7), aren't used !

  //row1
  22725, 22725, 22725, -22725,  // ; movq-> w06 w04 w02 w00
  29692, 12299, 12299, -29692,  //  ; w07 w05 w03 w01
  22725, -22725, 22725, 22725,  //; w14 w12 w10 w08
  -12299, 29692, -29692, -12299,  //; w15 w13 w11 w09
  31521, 17855, 26722, -31521,  //; w22 w20 w18 w16
  26722, 6270, -6270, -17855,   //; w23 w21 w19 w17
  17855, 6270, 6270, 26722,   //; w30 w28 w26 w24
  -31521, 26722, -17855, -31521,  // w31 w29 w27 w25

  //row2
  21407, 21407, 21407, -21407,  // ; movq-> w06 w04 w02 w00
  27969, 11585, 11585, -27969,  // ; w07 w05 w03 w01
  21407, -21407, 21407, 21407,  // ; w14 w12 w10 w08
  -11585, 27969, -27969, -11585,  //  ;w15 w13 w11 w09
  29692, 16819, 25172, -29692,  // ;w22 w20 w18 w16
  25172, 5906, -5906, -16819,   // ;w23 w21 w19 w17
  16819, 5906, 5906, 25172,     // ;w30 w28 w26 w24
  -29692, 25172, -16819, -29692,  //  ;w31 w29 w27 w25

  //row3
  19266, 19266, 19266, -19266,  //; movq-> w06 w04 w02 w00
  25172, 10426, 10426, -25172,  //; w07 w05 w03 w01
  19266, -19266, 19266, 19266,  //; w14 w12 w10 w08
  -10426, 25172, -25172, -10426,  //; w15 w13 w11 w09
  26722, 15137, 22654, -26722,  //; w22 w20 w18 w16
  22654, 5315, -5315, -15137,   //; w23 w21 w19 w17
  15137, 5315, 5315, 22654,   //; w30 w28 w26 w24
  -26722, 22654, -15137, -26722,  //; w31 w29 w27 w25

  //row4
  16384, 16384, 16384, -16384,  // ; movq-> w06 w04 w02 w00
  21407, 8867, 8867, -21407,    // w07 w05 w03 w01
  16384, -16384, 16384, 16384,  //; w14 w12 w10 w08
  -8867, 21407, -21407, -8867,  //; w15 w13 w11 w09
  22725, 12873, 19266, -22725,  //; w22 w20 w18 w16
  19266, 4520, -4520, -12873,   //; w23 w21 w19 w17
  12873, 4520, 4520, 19266,   //; w30 w28 w26 w24
  -22725, 19266, -12873, -22725,  //w31 w29 w27 w25

  //row5
  19266, 19266, 19266, -19266,  //; movq-> w06 w04 w02 w00
  25172, 10426, 10426, -25172,  //; w07 w05 w03 w01
  19266, -19266, 19266, 19266,  //; w14 w12 w10 w08
  -10426, 25172, -25172, -10426,  //; w15 w13 w11 w09
  26722, 15137, 22654, -26722,  //; w22 w20 w18 w16
  22654, 5315, -5315, -15137,   //; w23 w21 w19 w17
  15137, 5315, 5315, 22654,   //; w30 w28 w26 w24
  -26722, 22654, -15137, -26722,  //; w31 w29 w27 w25

  //row6
  21407, 21407, 21407, -21407,  // ; movq-> w06 w04 w02 w00
  27969, 11585, 11585, -27969,  // ; w07 w05 w03 w01
  21407, -21407, 21407, 21407,  // ; w14 w12 w10 w08
  -11585, 27969, -27969, -11585,  //  ;w15 w13 w11 w09
  29692, 16819, 25172, -29692,  // ;w22 w20 w18 w16
  25172, 5906, -5906, -16819,   // ;w23 w21 w19 w17
  16819, 5906, 5906, 25172,     // ;w30 w28 w26 w24
  -29692, 25172, -16819, -29692,  //  ;w31 w29 w27 w25

  //row7
  22725, 22725, 22725, -22725,  // ; movq-> w06 w04 w02 w00
  29692, 12299, 12299, -29692,  //  ; w07 w05 w03 w01
  22725, -22725, 22725, 22725,  //; w14 w12 w10 w08
  -12299, 29692, -29692, -12299,  //; w15 w13 w11 w09
  31521, 17855, 26722, -31521,  //; w22 w20 w18 w16
  26722, 6270, -6270, -17855,   //; w23 w21 w19 w17
  17855, 6270, 6270, 26722,   //; w30 w28 w26 w24
  -31521, 26722, -17855, -31521}; // w31 w29 w27 w25


#define INP eax   // pointer to (short *blk)
#define OUT ecx   // pointer to output (temporary store space qwTemp[])
#define TABLE ebx // pointer to tab_i_01234567[]
#define round_inv_row edx
#define round_inv_col edx

#define ROW_STRIDE 16 // for 8x8 matrix transposer

// private variables and functions

//temporary storage space, 8x8 of shorts
static __int64 qwTemp[32];
static __int64 scratch1, scratch2; // scratch variables

// v0.16B34, iDCT operation has been concatenated into ONE function

//__inline static void idct_mmx32_rows( short *blk ); // transform rows
//__inline static void idct_mmx32_cols( short *blk ); // transform "columns"
  // the "column" transform actually transforms rows, it is
  // identical to the row-transform except for the ROUNDING
  // and SHIFTING coefficients.


// public interface, performs iDCT on an array of 64 shorts, blk[8][8]
// input data is overwritten, of course
void 
j_rev_dct( short *blk ) // perform complete iDCT on 8x8 blk[]
{
  // this routine performs the COMPLETE iDCT 8x8 operation on blk[]
  //
  // 1) iDCT row transform (1st pass, input data is blk[])
  //    for( i = 0; i < 8; ++ i)
  //      DCT_8_INV_ROW_2( blk[i*8], qwTemp[i] );
  //
  //    Output of step 1 -> qwTemp[].  Matrix is transposed in-place.
  //
  // 2) iDCT row transform (2nd pass, input data is qwTemp[])
  //    for( i = 0; i < 8; i=i+2)
  //      DCT_8_INV_ROW_2( qwTemp[i*8], blk[i*8] );
  //
  //    Final output -> blk[].  No need to transpose.
  //
  //
  // v0.16B34 - "pseudo" in-place output transpose 
  //  -------
  //   minor optimization
  //
  //  The implementation now processes TWO rows per loop iteration.
  //  The 1st output row (y+0) is stored into scratch variables
  //
  //  The 2nd output row (y+1) is combined with the 1st-row, and
  //  transposed on the fly.  The results are written vertically into
  //  qwTemp[]  (8 dword stores.)
  //
  //     The algorithm in C look something like this
  //
  //  for( i = 0; i < 8; i = i + 2)
  //  {
  //      DCT_8_INV_ROW_1( blk[i*8] );  --> store in scratch vars
  //
  //      DCT_8_INV_ROW_1( blk[(i+1)*8] ) --> combine with scratch vars
  //      // transpose in place
  //      // stores don't occur in this order
  //      store DWORD -> [ cols (1,0) "row 0"]; 
  //      store DWORD -> [ cols (1,0) "row 2"];
  //      store DWORD -> [ cols (1,0) "row 1"];