#include <stdio.h>
#include "dct.h"

/* This code uses the algorithm from Page 52 of Pennebaker's JPEG book 
   Reference to the jpeg groups implementation of this algorithm was also
   used */


#define DCTSIZE 8


#ifndef USE_FLOAT_PRECISION
const short FIX_0_382683433 = 98;
const short FIX_0_541196100 = 139;
const short FIX_0_707106781 = 181;
const short FIX_1_306562965 = 334;
const short FIX_1_082392200 = 277;
const short FIX_1_414213562 = 362;
const short FIX_1_847759065 = 473;
const short FIX_2_613125930 = 669;

#else
const float FIX_0_382683433 = 0.382683433;
const float FIX_0_541196100 = 0.541196100;
const float FIX_0_707106781 = 0.707106781;
const float FIX_1_306562965 = 1.306562965;
const float FIX_1_082392200 = 1.082392200;
const float FIX_1_414213562 = 1.414213562;
const float FIX_1_847759065 = 1.847759065;
const float FIX_2_613125930 = 2.613125930;

#endif

/* To get the correct DCT values we have to use
   Yp / (8 * dscale) where dscale is the output of dctscale.m */


#ifdef FLOAT_QNT
static float  quantmat[8][8] = 
{
	{ 8.0000, 11.0963, 10.4525,  9.4070,  8.0000, 6.2856, 4.3296, 2.2072},
	{11.0963, 15.3910, 14.4980, 13.0479, 11.0963, 8.7183, 6.0053, 3.0615},
	{10.4525, 14.4980, 13.6569, 12.2908, 10.4525, 8.2125, 5.6569, 2.8838},
	{ 9.4070, 13.0479, 12.2908, 11.0615,  9.4070, 7.3910, 5.0910, 2.5954},
	{ 8.0000, 11.0963, 10.4525,  9.4070,  8.0000, 6.2856, 4.3296, 2.2072},
	{ 6.2856,  8.7183,  8.2125,  7.3910,  6.2856, 4.9385, 3.4017, 1.7342},
	{ 4.3296,  6.0053,  5.6569,  5.0910,  4.3296, 3.4017, 2.3431, 1.1945},
	{ 2.2072,  3.0615,  2.8838,  2.5954,  2.2072, 1.7342, 1.1945, 0.6090}
};

#else

/* Using printqnt.m:  Each element is 4096 or 2^12 times the numbers
   above */
static int int_quantmat[8][8] = 
{
	{32768,   45451,   42813,   38531,   32768,   25746,   17734,    9041},
	{45451,   63042,   59384,   53444,   45451,   35710,   24598,   12540}, 
	{42813,   59384,   55938,   50343,   42813,   33638,   23170,   11812}, 
	{38531,   53444,   50343,   45308,   38531,   30274,   20853,   10631}, 
	{32768,   45451,   42813,   38531,   32768,   25746,   17734,    9041}, 
	{25746,   35710,   33638,   30274,   25746,   20228,   13933,    7103}, 
	{17734,   24598,   23170,   20853,   17734,   13933,    9598,    4893}, 
	{ 9041,   12540,   11812,   10631,    9041,    7103,    4893,    2494} 

};
#endif

void qnt(short *data, int width, int dcq, int acq)
{
	short *dptr = data;
#ifdef FLOAT_QNT
	float *qptr = (float *)quantmat;
#else
	int  *qptr = (int *)int_quantmat;
#endif
	short tmp = *data;

	for(int j=8; j>0; j--) {
		for(int i=8; i>0; i--) {
#ifdef FLOAT_QNT
  			*dptr = short( float(*dptr) / (*qptr * float(acq)));
#else
  			*dptr = short( (int(*dptr) << 12) / (*qptr * acq));
#endif
			dptr++;
			qptr++;
		}

		dptr += width-DCTSIZE;
	}
	/* Put back DC value */
#ifdef FLOAT_QNT
	data[0] =  short( float(tmp) / (quantmat[0][0] * float(dcq)));
#else
	data[0] =  short( (int(tmp) << 12) / (int_quantmat[0][0] * dcq));
#endif

}


void iqnt(short *data, int width, int dcq, int acq)
{
	short *dptr = data;
#ifdef FLOAT_QNT
	float *qptr = (float *)quantmat;
#else
	int *qptr = (int *)int_quantmat;
#endif
	short tmp = *data;

	for(int j=8; j>0; j--) {
		for(int i=8; i>0; i--) {
#ifdef FLOAT_QNT
			*dptr = short( float(*dptr) * (*qptr * float(acq)));
#else
			*dptr = short( (int(*dptr) * (*qptr) * acq) >> 12 );
#endif
			dptr++;
			qptr++;
		}

		dptr += width-DCTSIZE;
	}
	/* Put back DC value */
#ifdef FLOAT_QNT
	data[0] =  short( float(tmp) * (quantmat[0][0] * float(dcq)));
#else
	data[0] =  short( (int(tmp) * int_quantmat[0][0] * dcq) >> 12 );
#endif
	
}



#ifdef USE_FLOAT_PRECISION
inline short MULTIPLY(short var, float num)
{
	return short( float(var) * num);
}
#else
inline short MULTIPLY(short var, short num)
{
	return short( (int(var) * int(num)) >> 8);
}
#endif



/* width is to extract data from a matrix which may be larger than 8x8 */
void dct(short *data, int width)
{
	short tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
	short tmp10, tmp11, tmp12, tmp13;
	short z1, z2, z3, z4, z5, z11, z13;
	short *dataptr;
	int ctr;

	/* Pass 1: process rows. */

	dataptr = data;
	for (ctr = DCTSIZE-1; ctr >= 0; ctr--) {
		tmp0 = dataptr[0] + dataptr[7];
		tmp7 = dataptr[0] - dataptr[7];
		tmp1 = dataptr[1] + dataptr[6];
		tmp6 = dataptr[1] - dataptr[6];
		tmp2 = dataptr[2] + dataptr[5];
		tmp5 = dataptr[2] - dataptr[5];
		tmp3 = dataptr[3] + dataptr[4];
		tmp4 = dataptr[3] - dataptr[4];
    
		/* Even part */
    
		tmp10 = tmp0 + tmp3;	/* phase 2 */
		tmp13 = tmp0 - tmp3;
		tmp11 = tmp1 + tmp2;
		tmp12 = tmp1 - tmp2;
    
		dataptr[0] = tmp10 + tmp11; /* phase 3 */
		dataptr[4] = tmp10 - tmp11;
   
		z1 = MULTIPLY(tmp12 + tmp13, FIX_0_707106781); /* c4 */
		dataptr[2] = tmp13 + z1;	/* phase 5 */
		dataptr[6] = tmp13 - z1;
    
		/* Odd part */

		tmp10 = tmp4 + tmp5;	/* phase 2 */
		tmp11 = tmp5 + tmp6;
		tmp12 = tmp6 + tmp7;

		z5 = MULTIPLY(tmp10 - tmp12, FIX_0_382683433); /* c6 */
		z2 = MULTIPLY(tmp10, FIX_0_541196100) + z5; /* c2-c6 */
		z4 = MULTIPLY(tmp12, FIX_1_306562965) + z5; /* c2+c6 */
		z3 = MULTIPLY(tmp11, FIX_0_707106781); /* c4 */


		z11 = tmp7 + z3;		/* phase 5 */
		z13 = tmp7 - z3;

		dataptr[5] = z13 + z2;	/* phase 6 */
		dataptr[3] = z13 - z2;
		dataptr[1] = z11 + z4;
		dataptr[7] = z11 - z4;

		dataptr += width;

	}

	/* Pass 2: process columns. */

	dataptr = data;
	for (ctr = DCTSIZE-1; ctr >= 0; ctr--) {

		tmp0 = dataptr[width*0] + dataptr[width*7];
		tmp7 = dataptr[width*0] - dataptr[width*7];
		tmp1 = dataptr[width*1] + dataptr[width*6];
		tmp6 = dataptr[width*1] - dataptr[width*6];
		tmp2 = dataptr[width*2] + dataptr[width*5];
		tmp5 = dataptr[width*2] - dataptr[width*5];
		tmp3 = dataptr[width*3] + dataptr[width*4];
		tmp4 = dataptr[width*3] - dataptr[width*4];

		/* Even part */
    
		tmp10 = tmp0 + tmp3;	/* phase 2 */
		tmp13 = tmp0 - tmp3;
		tmp11 = tmp1 + tmp2;
		tmp12 = tmp1 - tmp2;


		dataptr[width*0] = tmp10 + tmp11; /* phase 3 */
		dataptr[width*4] = tmp10 - tmp11;

		z1 = MULTIPLY(tmp12 + tmp13, FIX_0_707106781); /* c4 */

		dataptr[width*2] = tmp13 + z1; /* phase 5 */
		dataptr[width*6] = tmp13 - z1;

		/* Odd part */

		tmp10 = tmp4 + tmp5;	/* phase 2 */
		tmp11 = tmp5 + tmp6;
		tmp12 = tmp6 + tmp7;


		z5 = MULTIPLY(tmp10 - tmp12, FIX_0_382683433); /* c6 */
		z2 = MULTIPLY(tmp10, FIX_0_541196100) + z5; /* c2-c6 */
		z4 = MULTIPLY(tmp12, FIX_1_306562965) + z5; /* c2+c6 */
		z3 = MULTIPLY(tmp11, FIX_0_707106781); /* c4 */

		z11 = tmp7 + z3;		/* phase 5 */
		z13 = tmp7 - z3;

		dataptr[width*5] = z13 + z2; /* phase 6 */
		dataptr[width*3] = z13 - z2;
		dataptr[width*1] = z11 + z4;
		dataptr[width*7] = z11 - z4;

			


		dataptr++;			/* advance pointer to next column */
	}

	



}



void idct(short *data, int width)
{
  short tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
  short tmp10, tmp11, tmp12, tmp13;
  short z5, z10, z11, z12, z13;
  
  short *inptr, *outptr;
  short workspace[64];
  short *wsptr;
  int ctr;
  

  /* Pass 1 */
  inptr = data;
  wsptr = workspace;


  for (ctr = DCTSIZE; ctr > 0; ctr--) {

	  /* Shortens things if input coefficients are 0 */
    if ((inptr[width*1] | inptr[width*2] | inptr[width*3] |
	 inptr[width*4] | inptr[width*5] | inptr[width*6] |
	 inptr[width*7]) == 0) {

      /* AC terms all zero */
      short dcval =  inptr[width*0]; 


      wsptr[DCTSIZE*0] = dcval;
      wsptr[DCTSIZE*1] = dcval;
      wsptr[DCTSIZE*2] = dcval;
      wsptr[DCTSIZE*3] = dcval;
      wsptr[DCTSIZE*4] = dcval;
      wsptr[DCTSIZE*5] = dcval;
      wsptr[DCTSIZE*6] = dcval;
      wsptr[DCTSIZE*7] = dcval;
      
      inptr++;			/* advance pointers to next column */
      /* quantptr++; */
      wsptr++;
      continue;
    }
    

    /* Even part */
    tmp0 = inptr[width*0];
    tmp1 = inptr[width*2];
    tmp2 = inptr[width*4];
    tmp3 = inptr[width*6];


    tmp10 = tmp0 + tmp2;	/* phase 3 */
    tmp11 = tmp0 - tmp2;

    tmp13 = tmp1 + tmp3;	/* phases 5-3 */
    tmp12 = MULTIPLY(tmp1 - tmp3, FIX_1_414213562) - tmp13; /* 2*c4 */

    tmp0 = tmp10 + tmp13;	/* phase 2 */
    tmp3 = tmp10 - tmp13;
    tmp1 = tmp11 + tmp12;
    tmp2 = tmp11 - tmp12;
    

    /* Odd part */
    tmp4 = inptr[width*1];
    tmp5 = inptr[width*3];
    tmp6 = inptr[width*5];
    tmp7 = inptr[width*7];


    z13 = tmp6 + tmp5;		/* phase 6 */
    z10 = tmp6 - tmp5;
    z11 = tmp4 + tmp7;
    z12 = tmp4 - tmp7;

    tmp7 = z11 + z13;		/* phase 5 */
    tmp11 = MULTIPLY(z11 - z13, FIX_1_414213562); /* 2*c4 */

    z5 = MULTIPLY(z10 + z12, FIX_1_847759065); /* 2*c2 */
    tmp10 = MULTIPLY(z12, FIX_1_082392200) - z5; /* 2*(c2-c6) */
    tmp12 = MULTIPLY(z10, - FIX_2_613125930) + z5; /* -2*(c2+c6) */

    tmp6 = tmp12 - tmp7;	/* phase 2 */
    tmp5 = tmp11 - tmp6;
    tmp4 = tmp10 + tmp5;

    wsptr[DCTSIZE*0] = (int) (tmp0 + tmp7);
    wsptr[DCTSIZE*7] = (int) (tmp0 - tmp7);
    wsptr[DCTSIZE*1] = (int) (tmp1 + tmp6);
    wsptr[DCTSIZE*6] = (int) (tmp1 - tmp6);
    wsptr[DCTSIZE*2] = (int) (tmp2 + tmp5);
    wsptr[DCTSIZE*5] = (int) (tmp2 - tmp5);
    wsptr[DCTSIZE*4] = (int) (tmp3 + tmp4);
    wsptr[DCTSIZE*3] = (int) (tmp3 - tmp4);

    inptr++;			/* advance pointers to next column */
    /* quantptr++; */
    wsptr++;
  }





  /* Pass 2: process rows from work array, store into output array. */
  /* Note that we must descale the results by a factor of 8 == 2**3, */
  /* and also undo the PASS1_BITS scaling. */

  wsptr = workspace;
  outptr = data;
  for (ctr = 0; ctr < DCTSIZE; ctr++) {


    /* Even part */

    tmp10 = ((short) wsptr[0] + (short) wsptr[4]);
    tmp11 = ((short) wsptr[0] - (short) wsptr[4]);

    tmp13 = ((short) wsptr[2] + (short) wsptr[6]);
    tmp12 = MULTIPLY((short) wsptr[2] - (short) wsptr[6], FIX_1_414213562)
	    - tmp13;

    tmp0 = tmp10 + tmp13;
    tmp3 = tmp10 - tmp13;
    tmp1 = tmp11 + tmp12;
    tmp2 = tmp11 - tmp12;

    /* Odd part */

    z13 = (short) wsptr[5] + (short) wsptr[3];
    z10 = (short) wsptr[5] - (short) wsptr[3];
    z11 = (short) wsptr[1] + (short) wsptr[7];
    z12 = (short) wsptr[1] - (short) wsptr[7];

    tmp7 = z11 + z13;		/* phase 5 */
    tmp11 = MULTIPLY(z11 - z13, FIX_1_414213562); /* 2*c4 */

    z5 = MULTIPLY(z10 + z12, FIX_1_847759065); /* 2*c2 */
    tmp10 = MULTIPLY(z12, FIX_1_082392200) - z5; /* 2*(c2-c6) */
    tmp12 = MULTIPLY(z10, - FIX_2_613125930) + z5; /* -2*(c2+c6) */

    tmp6 = tmp12 - tmp7;	/* phase 2 */
    tmp5 = tmp11 - tmp6;
    tmp4 = tmp10 + tmp5;

    /* Final output stage: scale down by a factor of 8 and range-limit */


    outptr[0] = (tmp0 + tmp7) >>  6;
    outptr[7] = (tmp0 - tmp7) >>  6;
    outptr[1] = (tmp1 + tmp6) >>  6;
    outptr[6] = (tmp1 - tmp6) >>  6;
    outptr[2] = (tmp2 + tmp5) >>  6;
    outptr[5] = (tmp2 - tmp5) >>  6;
    outptr[4] = (tmp3 + tmp4) >>  6;
    outptr[3] = (tmp3 - tmp4) >>  6;


    wsptr += DCTSIZE;		/* advance pointer to next row */
    outptr += width;
  }



}