idct_arai.cc

PIXIL is a small footprint operating environment, complete with PDA PIM applications, a browser and

#include "video12/modules/idct_arai.hh"

#include <math.h>
#include <assert.h>
#include <iostream.h>
#include <iomanip.h>

static Pixel* clip2_0_255;
static Pixel  clip2_0_255_data[1024];

static class InitCLIP_DCT
{
public:
  InitCLIP_DCT()
    {
      clip2_0_255 = &clip2_0_255_data[512];

      for (int i=-512;i<512;i++)
        {
          int val=i;
          if (val<0) val=0;
          if (val>255) val=255;
          clip2_0_255[i]=val;
        }
    }
} dummyclass_weiour;


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

// static int c1,c2;

const int FDCT_M1Scale      = 30;
const int FDCT_PostM1_VK    = 12;
const int FDCT_PostM1_NK    = 16;
const int FDCT_PostM1_Shift = 14;

const int FDCT_M2Scale      = 30;
const int FDCT_PostM2_VK    = 16;
const int FDCT_PostM2_NK    = 15;
const int FDCT_PostM2_Shift = 31;

const int FDCT_M3Scale      = 30;
const int FDCT_PostM3_VK    = 17;
const int FDCT_PostM3_NK    = 15;
const int FDCT_PostM3_Shift = 30;



const int IDCT_M3Scale      = 33;  // DOCS !
const int IDCT_PostM3_VK    = 10;
const int IDCT_PostM3_NK    = 20;
const int IDCT_PostM3_Shift = 13;

const int IDCT_M2Scale      = 28;
const int IDCT_PostM2_VK    = 12;
const int IDCT_PostM2_NK    = 16;
const int IDCT_PostM2_Shift = 32;

const int IDCT_M1Scale      = 28;
const int IDCT_PostM1_VK    = 16;
const int IDCT_PostM1_NK    = 14;
const int IDCT_PostM1_Shift = 30;


typedef int FDCTInt;
typedef int IDCTInt;



#define SC(x,u,v) { D(x); ibuf[u][v] = ((x) >> (FDCTScale-FDCT_STAGE1_NK)); }

#define D(x,s) { cout << setbase(16) << setfill('0') << setw(16) << (x) << setbase(10) << " "; \
  ibuf[c1][c2]=x; c2++; if (c2==8) { c1++; c2=0; } cout << setfill(' '); }

//#define DESCALE(x,f)   ((x<0) ? ((x - ((1LL<<(f-1)))  ) >> (f)) : ((x + ((1LL<<(f-1)))  ) >> (f)))
#define DESCALE(x,f)    ((x + ((1LL<<(f)) -1 )  ) >> (f))

//#define DESCALE(x,f)      (x >> (f))
//#define DESCALE(x,f)   ((x<0) ? ((x - (1LL<<(f))/2  ) >> (f)) : ((x + (1LL<<(f))/2  ) >> (f)))
//#define DESCALE(x,f)    ((x + ((1LL<<(f))/2 )  ) >> (f))

#define IDESCALE(x,f)    ((x + ((1LL<<(f-1)) )  ) >> (f))
#define IDESCALENeg(x,f) ((x - ((1LL<<(f-1)) )  ) >> (f))

/* Eine 32x32 -> 64bit Multiplikation durchfuehren */
#define MUL32(a,b) (((long long)(a))*((long)(b)))

const long long FDCT_M1ScaleFact = (1LL<<FDCT_M1Scale);
const long long FDCT_M2ScaleFact = (1LL<<FDCT_M2Scale);
const long long FDCT_M3ScaleFact = (1LL<<FDCT_M3Scale);

const long long IDCT_M3ScaleFact = (1LL<<IDCT_M3Scale);
const long long IDCT_M2ScaleFact = (1LL<<IDCT_M2Scale);
const long long IDCT_M1ScaleFact = (1LL<<IDCT_M1Scale);

#ifndef PI
#define PI 3.1415926535
#endif

typedef long long DCTInt2;

static bool IsInitialized = false;

static DCTInt2  A1;
static DCTInt2 MA2;
static DCTInt2  A3;
static DCTInt2  A4;
static DCTInt2  A5;

static DCTInt2  A1b;
static DCTInt2 MA2b;
static DCTInt2  A3b;
static DCTInt2  A4b;
static DCTInt2  A5b;

static IDCTInt C2;
static IDCTInt C4;
static IDCTInt C6;
static IDCTInt Q;
static IDCTInt R;

static IDCTInt C2b;
static IDCTInt C4b;
static IDCTInt C6b;
static IDCTInt Qb;
static IDCTInt Rb;

static double OutputScaleFactor[8];  // fuer FDCT
static double InputScaleFactor[8];   // fuer IDCT

static long long d_fact[8][8];

static void Init()
{
  if (IsInitialized)
    return;

  /* Beachte: Jeder Wert wird exakt einmal zugewiesen und nicht veraendert, dadurch
     ist der Code Multithread-sicher. */

  // FDCT

  A1 = (DCTInt2)(sqrt(0.5)*FDCT_M1ScaleFact+0.5);
  MA2 = (DCTInt2)((-(cos(PI/8.0) - cos(3.0*PI/8.0)))*FDCT_M1ScaleFact-0.5);
  A3 = (DCTInt2)(sqrt(0.5)*FDCT_M1ScaleFact+0.5);
  A4 = (DCTInt2)((cos(PI/8.0) + cos(3.0*PI/8.0))*FDCT_M1ScaleFact+0.5);
  A5 = (DCTInt2)(cos(3.0*PI/8.0)*FDCT_M1ScaleFact+0.5);

  A1b = (DCTInt2)(sqrt(0.5)*FDCT_M2ScaleFact+0.5);
  MA2b = (DCTInt2)((-(cos(PI/8.0) - cos(3.0*PI/8.0)))*FDCT_M2ScaleFact-0.5);
  A3b = (DCTInt2)(sqrt(0.5)*FDCT_M2ScaleFact+0.5);
  A4b = (DCTInt2)((cos(PI/8.0) + cos(3.0*PI/8.0))*FDCT_M2ScaleFact+0.5);
  A5b = (DCTInt2)(cos(3.0*PI/8.0)*FDCT_M2ScaleFact+0.5);

  OutputScaleFactor[0] = sqrt(0.5)/2.0;
  {for (int i=1;i<8;i++) OutputScaleFactor[i] = 0.25/cos(PI*(double)i/16.0);}


  // IDCT

  C2  = (DCTInt2)(2.0*cos(PI/8.0)*IDCT_M2ScaleFact+0.5);
  C4  = (DCTInt2)(sqrt(2.0)*IDCT_M2ScaleFact+0.5);
  C6  = (DCTInt2)(2.0*cos(3.0*PI/8.0)*IDCT_M2ScaleFact+0.5);
  Q   = (DCTInt2)((-2.0*cos(PI/8.0) + 2.0*cos(3.0*PI/8.0))*IDCT_M2ScaleFact-0.5);
  R   = (DCTInt2)(( 2.0*cos(PI/8.0) + 2.0*cos(3.0*PI/8.0))*IDCT_M2ScaleFact+0.5);
  /*
    cout << "C2 " << 2.0*cos(PI/8.0) << endl;
    cout << "C4 " << sqrt(2.0) << endl;
    cout << "C6 " << 2.0*cos(3.0*PI/8.0) << endl;
    cout << "Q  " << (-2.0*cos(PI/8.0) + 2.0*cos(3.0*PI/8.0)) << endl;
    cout << "R  " << ( 2.0*cos(PI/8.0) + 2.0*cos(3.0*PI/8.0)) << endl;
  */
  C2b = (DCTInt2)(2.0*cos(PI/8.0)*IDCT_M1ScaleFact+0.5);
  C4b = (DCTInt2)(sqrt(2.0)*IDCT_M1ScaleFact+0.5);
  C6b = (DCTInt2)(2.0*cos(3.0*PI/8.0)*IDCT_M1ScaleFact+0.5);
  Qb  = (DCTInt2)((-2.0*cos(PI/8.0) + 2.0*cos(3.0*PI/8.0))*IDCT_M1ScaleFact-0.5);
  Rb  = (DCTInt2)(( 2.0*cos(PI/8.0) + 2.0*cos(3.0*PI/8.0))*IDCT_M1ScaleFact+0.5);
  
  InputScaleFactor[0] = 0.5*0.5*sqrt(2);
  {for (int i=1;i<8;i++) InputScaleFactor[i] = 0.5*cos(PI*((double)i)/16.0);}

  IsInitialized=true;


#if 1
  for (int u=0;u<8;u++)
    for (int v=0;v<8;v++)
      {
	d_fact[v][u] = (DCTInt2)(InputScaleFactor[u]*InputScaleFactor[v]*IDCT_M3ScaleFact+0.5);
      }
#else
  for (int u=0;u<8;u++)
    for (int v=0;v<8;v++)
      {
	d_fact[v][u] = (DCTInt2)(InputScaleFactor[u]*InputScaleFactor[v] * fact[v][u] *
                                 IDCT_M3ScaleFact+0.5);
      }
#endif





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


static struct DummyInit
{
  DummyInit() { Init(); }
} dummyinit;



void IDCT_Int(const short* in[8], short* out[8])
{
  IDCTInt fin[8][8];
  IDCTInt tmp[8][8];

  IDCTInt a0,a1,a2,a3,a4,a5,a6,a7;
  IDCTInt       b2,   b4,b5,b6   ;
  IDCTInt n0,n1,n2,n3            ;
  IDCTInt          m3,m4,m5,m6,m7;
   DCTInt2 tmp1,tmp2,tmp3,tmp4;


  // Zeilentransformationen

  {for (int u=0;u<8;u++)
    {
      if (             !in[1][u] && !in[2][u] && !in[3][u] &&
          !in[4][u] && !in[5][u] && !in[6][u] && !in[7][u])
        {
          tmp[0][u] = tmp[1][u] = tmp[2][u] = tmp[3][u] = 
          tmp[4][u] = tmp[5][u] = tmp[6][u] = tmp[7][u] = in[0][u]<<13;
        }
      else
        {
          {for (int v=0;v<8;v++)
            fin[v][u] = IDESCALE(MUL32(d_fact[v][u] , in[v][u]) , IDCT_PostM3_Shift);}

          a0 = fin[0][u];
          a1 = fin[4][u];
          a2 = fin[2][u]-fin[6][u];
          a3 = fin[2][u]+fin[6][u];
          a4 = fin[5][u]-fin[3][u];
          tmp1 = fin[1][u]+fin[7][u];
          tmp2 = fin[3][u]+fin[5][u];
          a5 = tmp1 - tmp2;
          a6 = fin[1][u]-fin[7][u];
          a7 = tmp1 + tmp2;
          
          tmp4 = MUL32(C6,(a4+a6));
          b2 = IDESCALE(MUL32(C4 , a2)       , IDCT_PostM2_Shift);
          b4 = IDESCALE(MUL32(Q  , a4) -tmp4 , IDCT_PostM2_Shift);
          b5 = IDESCALE(MUL32(C4 , a5)       , IDCT_PostM2_Shift);
          b6 = IDESCALE(MUL32(R  , a6) -tmp4 , IDCT_PostM2_Shift);

          a0 >>= -(IDCT_M2Scale - IDCT_PostM2_Shift);
          a1 >>= -(IDCT_M2Scale - IDCT_PostM2_Shift);
          a3 >>= -(IDCT_M2Scale - IDCT_PostM2_Shift);
          a7 >>= -(IDCT_M2Scale - IDCT_PostM2_Shift);
          
          tmp3 = b6-a7;
          n0 = tmp3-b5;
          n1 = a0-a1;
          n2 = b2-a3;
          n3 = a0+a1;
          
          m3 = n1+n2;
          m4 = n3+a3;
          m5 = n1-n2;
          m6 = n3-a3;
          m7 = b4-n0;
          
          tmp[0][u] = m4+a7;
          tmp[1][u] = m3+tmp3;
          tmp[2][u] = m5-n0;
          tmp[3][u] = m6-m7;
          tmp[4][u] = m6+m7;
          tmp[5][u] = m5+n0;
          tmp[6][u] = m3-tmp3;
          tmp[7][u] = m4-a7;
        }
    }}


  // Spaltentransformationen

  {for (int v=0;v<8;v++)
    {
      a0 = tmp[v][0];
      a1 = tmp[v][4];
      a2 = tmp[v][2]-tmp[v][6];
      a3 = tmp[v][2]+tmp[v][6];
      a4 = tmp[v][5]-tmp[v][3];
      tmp1 = tmp[v][1]+tmp[v][7];
      tmp2 = tmp[v][3]+tmp[v][5];
      a5 = tmp1 - tmp2;