/*************************************************************
Copyright (C) 1990, 1991, 1993 Andy C. Hung, all rights reserved.
PUBLIC DOMAIN LICENSE: Stanford University Portable Video Research
Group. If you use this software, you agree to the following: This
program package is purely experimental, and is licensed "as is".
Permission is granted to use, modify, and distribute this program
without charge for any purpose, provided this license/ disclaimer
notice appears in the copies.  No warranty or maintenance is given,
either expressed or implied.  In no event shall the author(s) be
liable to you or a third party for any special, incidental,
consequential, or other damages, arising out of the use or inability
to use the program for any purpose (or the loss of data), even if we
have been advised of such possibilities.  Any public reference or
advertisement of this source code should refer to it as the Portable
Video Research Group (PVRG) code, and not by any author(s) (or
Stanford University) name.
*************************************************************/
/*
************************************************************
transform.c

This file contains the reference DCT, the zig-zag and quantization
algorithms.

************************************************************
*/

/*LABEL transform.c */

#include "globals2.h"
#include "dct2.h"
#include <math.h>

/*PUBLIC*/

//extern void ReferenceDct();
//extern void ReferenceIDct();
//extern void TransposeMatrix();
//extern void CCITTQuantize_en();
//extern void CCITTFlatQuantize_en();
extern void ICCITTFlatQuantize2();
extern void ICCITTQuantize2();
//extern void BoundDctMatrix_en();
extern void BoundIDctMatrix2();
//extern void FlatBoundQuantizeMatrix_en();
//extern void ZigzagMatrix_en();
extern void IZigzagMatrix2();
//extern void PrintMatrix();
//extern void ClearMatrix();

//static void DoubleReferenceDct1D();
//static void DoubleReferenceIDct1D();
//static void DoubleTransposeMatrix();

/*PRIVATE*/

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

static int zigzag_index2[] =
{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};


#define MakeMatrix() (int *) calloc(BLOCKSIZE,sizeof(int))
#define FixedMultiply(s,x,y)  x = ((x * y) >> s);
#define DCT_OFFSET 128


/*START*/

/*BFUNC

ReferenceDct() does a reference DCT on the input (matrix) and output2
(new matrix).

EFUNC*/

//void ReferenceDct(matrix,newmatrix)
//     int *matrix;
//     int *newmatrix;
//void ReferenceDct(int *matrix,int *newmatrix)
//{
//  BEGIN("ReferenceDct");
//  int *mptr;
//  double *sptr,*dptr;
//  double sourcematrix[BLOCKSIZE],destmatrix[BLOCKSIZE];
//
//  for(sptr=sourcematrix,mptr=matrix;mptr<matrix+BLOCKSIZE;mptr++)
//    {
//      *(sptr++) = (double) *mptr;
//    }
//  for(dptr = destmatrix,sptr=sourcematrix;
//      sptr<sourcematrix+BLOCKSIZE;sptr+=BLOCKWIDTH)
//
//    {
//      DoubleReferenceDct1D(sptr,dptr);
//      dptr+=BLOCKWIDTH;
//    }
//  DoubleTransposeMatrix(destmatrix,sourcematrix);
//  for(dptr = destmatrix,sptr=sourcematrix;
//      sptr<sourcematrix+BLOCKSIZE;sptr+=BLOCKWIDTH)
//    {
//      DoubleReferenceDct1D(sptr,dptr);
//      dptr+=BLOCKWIDTH;
//    }
//  DoubleTransposeMatrix(destmatrix,sourcematrix);
//  for(sptr = sourcematrix,mptr=newmatrix;
//      mptr<newmatrix+BLOCKSIZE;sptr++)
//    {    /* NB: Inversion on counter */
//      *(mptr++) = (int) (*sptr > 0 ? (*(sptr)+0.5):(*(sptr)-0.5));
//    }
//}

/*BFUNC

DoubleReferenceDCT1D() does a 8 point dct on an array of double
input and places the result in a double output2.

EFUNC*/

//static void DoubleReferenceDct1D(ivect,ovect)
//     double *ivect;
//     double *ovect;
//static void DoubleReferenceDct1D(double *ivect,double *ovect)
//{
//  BEGIN("DoubleReferenceDct1D");
//  double *mptr,*iptr,*optr;
//
//  for(mptr=DctMatrix2,optr=ovect;optr<ovect+BLOCKWIDTH;optr++)
//    {
//      for(*optr=0,iptr=ivect;iptr<ivect+BLOCKWIDTH;iptr++)
//	{
//	  *optr += *iptr*(*(mptr++));
//	}
//    }
//}

/*BFUNC

ReferenceIDct() is used to perform a reference 8x8 inverse dct.  It is
a balanced IDCT. It takes the input (matrix) and puts it into the
output2 (newmatrix).

EFUNC*/

//void ReferenceIDct(matrix,newmatrix)
//     int *matrix;
//     int *newmatrix;
//void ReferenceIDct(int *matrix,int *newmatrix)
//{
//  BEGIN("ReferenceIDct");
//  int *mptr;
//  double *sptr,*dptr;
//  double sourcematrix[BLOCKSIZE],destmatrix[BLOCKSIZE];
//
//  for(sptr = sourcematrix,mptr=matrix;mptr<matrix+BLOCKSIZE;mptr++)
//    {
//      *(sptr++) = (double) *mptr;
//    }
//  for(dptr = destmatrix,sptr=sourcematrix;
//      sptr<sourcematrix+BLOCKSIZE;sptr+=BLOCKWIDTH)
//    {
//      DoubleReferenceIDct1D(sptr,dptr);
//      dptr+=BLOCKWIDTH;
//    }
//  DoubleTransposeMatrix(destmatrix,sourcematrix);
//  for(dptr = destmatrix,sptr=sourcematrix;
//      sptr<sourcematrix+BLOCKSIZE;sptr+=BLOCKWIDTH)
//    {
//      DoubleReferenceIDct1D(sptr,dptr);
//      dptr+=BLOCKWIDTH;
//    }
//  DoubleTransposeMatrix(destmatrix,sourcematrix);
//  for(sptr = sourcematrix,mptr=newmatrix;mptr<newmatrix+BLOCKSIZE;sptr++)
//    {    /* NB: Inversion on counter */
//      *(mptr++) = (int) (*sptr > 0 ? (*(sptr)+0.5):(*(sptr)-0.5));
//    }
//}

/*BFUNC

DoubleReferenceIDct1D() does an 8 point inverse dct on ivect and
puts the output2 in ovect.

EFUNC*/

//static void DoubleReferenceIDct1D(ivect,ovect)
//     double *ivect;
//     double *ovect;
//static void DoubleReferenceIDct1D(double *ivect,double *ovect)
//{
//  BEGIN("DoubleReferenceIDct1D");
//  double *mptr,*iptr,*optr;
//
//  for(mptr = IDctMatrix2,optr=ovect;optr<ovect+BLOCKWIDTH;optr++)
//    {
//      for(*optr=0,iptr=ivect;iptr<ivect+BLOCKWIDTH;iptr++)
//	{
//	  *optr += *iptr*(*(mptr++));
//	}
//    }
//}

/*BFUNC

TransposeMatrix transposes an input matrix and puts the output2 in
newmatrix.

EFUNC*/

//void TransposeMatrix(matrix,newmatrix)
//     int *matrix;
//     int *newmatrix;
//void TransposeMatrix(int *matrix,int *newmatrix)
//{
//  BEGIN("TransposeMatrix");
//  int *tptr;
//
//  for(tptr=transpose_index2;tptr<transpose_index2+BLOCKSIZE;tptr++)
//    {
//      *(newmatrix++) = matrix[*tptr];
//    }
//}

/*BFUNC

DoubleTransposeMatrix transposes a double input matrix and puts the
double output2 in newmatrix.

EFUNC*/

//static void DoubleTransposeMatrix(matrix,newmatrix)
//     double *matrix;
//     double *newmatrix;
//static void DoubleTransposeMatrix(double *matrix,double *newmatrix)
//{
//  BEGIN("DoubleTransposeMatrix");
//  int *tptr;
//
//  for(tptr=transpose_index2;tptr<transpose_index2+BLOCKSIZE;tptr++)
//    {
//      *(newmatrix++) = matrix[*tptr];
//    }
//}  

/*BFUNC

CCITTQuantzie() quantizes the input matrix with a DC quantize step
and an AC quantize step.

EFUNC*/

//void CCITTQuantize_en(matrix,dcqfact,acqfact)
//     int *matrix;
//      dcqfact;
//     int acqfact;
//void CCITTQuantize_en(int *matrix,int dcqfact,int acqfact)
//{
//  BEGIN("CCITTQuantize_en");
//  int *mptr;
//
////#ifdef VERSION_1_0
////  dcqfact++;
////  acqfact++;
////  dcqfact <<=1;
////  acqfact <<=1;
////  *matrix = *matrix/dcqfact;
////  for(mptr=matrix+1;mptr<matrix+BLOCKSIZE;mptr++)
////  {
////	  *mptr = *mptr / acqfact;
////  }
////#else
//  if (dcqfact&1) /* Odd */
//    *matrix= *matrix/(dcqfact<<1);
//  else           /* Even */
//    {
//      if (*matrix>0)
//	*matrix= (*matrix+1)/(dcqfact<<1);
//      else
//	*matrix= (*matrix-1)/(dcqfact<<1);
//    }
//  if (acqfact&1) /* Odd */
//    {
//      acqfact<<=1;
//      for(mptr=matrix+1;mptr<matrix+BLOCKSIZE;mptr++)
//	*mptr = *mptr / acqfact;
//    }
//  else
//    {
//      acqfact<<=1;
//      for(mptr=matrix+1;mptr<matrix+BLOCKSIZE;mptr++)
//	{
//	  if (*mptr>0)
//	    *mptr = (*mptr+1) / acqfact;
//	  else
//	    *mptr = (*mptr-1) / acqfact;
//	}
//    }
////#endif
//}

/*BFUNC

CCITTFlatQuantize_en() quantizes the input matrix by a dc factor (flat)
and an acqfactor (thresholded).

EFUNC*/

//void CCITTFlatQuantize_en(matrix,dcqfact,acqfact)
//     int *matrix;
//     int dcqfact;
//     int acqfact;
//void CCITTFlatQuantize_en(int *matrix,int dcqfact,int acqfact)
//{
//  BEGIN("CCITTFlatQuantize_en");
//  int *mptr;
//
//
//  if (*matrix > 0) {*matrix = (*matrix + dcqfact/2)/ dcqfact;}
//  else {*matrix = (*matrix - dcqfact/2)/ dcqfact;}
//
////#ifdef VERSION_1_0
////  acqfact++;
////  acqfact <<=1;