//=======================================================================================
//
// encoder of still images
//
// vika1enc.c
//
// developer: Henry Guennadi Levkin
//
// uncompress file - png file
// compress file   - vik file
//
// encoding :
// vika1enc file.ppm [QPL] [QPC]
//
// result: file.vic in dir where is vika1enc
//
//=======================================================================================

#include "vika1enc.h"

#define MAX_BLK_DIM    32

//=======================================================================================
// global variables

FILE* gFp; // logging file
int gnResultSize; // size of compressed file

int gW;      // RGB image width
int gH;      // RGB image height

byte* gpInp; // input RGB image (from file)

byte* gpRed;
byte* gpGreen;
byte* gpBlue;

// luma croma parameters

int gQPLu; // luma quantization parameter

int gNBLu; // luma macroblock size
int gNTLu; // luma block transform dim

int gLuW;  // width  of Luma frame
int gLuH;  // height of Luma frame

int gQPCb; // chromaB quantization param

int gNBCb; // ChromaB macroblock size
int gNTCb; // ChromaB block transform dim

int gCbW;  // width  of ChromaB frame
int gCbH;  // height of ChromaB frame

int gQPCr; // ChromaR quantization parameter

int gNBCr; // ChromaR macroblock size
int gNTCr; // ChromaR block transform dim

int gCrW;  // width  of ChromaR frame
int gCrH;  // height of ChromaR frame

// matrices

byte* gpLu;  // luma frame
byte* gpCb;  // chroma Cb frame
byte* gpCr;  // chroma Cr frame

short* gpResidLu; // residual frame between initial and predicted Luma frames
short* gpResidCb; // residual frame between initial and predicted ChromaB frames
short* gpResidCr; // residual frame between initial and predicted ChromaR frames

byte* gpPredTypeLu;   // best prediction mode for every block in Luma frame
byte* gpPredTypeCb;   // best prediction mode for every block in ChromaB frame
byte* gpPredTypeCr;   // best prediction mode for every block in ChromaR frame

byte* gpPredLu;   // prediction frame for Luma
byte* gpPredCb;   // prediction frame for ChromaB
byte* gpPredCr;   // prediction frame for ChromaR

short* gpTQLu;    // transformed quantized residual Luma matrix
short* gpTQCb;    // transformed quantized residual ChromaB matrix
short* gpTQCr;    // transformed quantized residual ChromaR matrix

short* gpReordLu; // reordering transformed quantized residual Luma vector
short* gpReordCb; // reordering transformed quantized residual ChromaB vector
short* gpReordCr; // reordering transformed quantized residual ChromaR vector

byte* gpReconLu;  // reconstructed Luma frame
byte* gpReconCb;  // reconstructed ChromaB frame
byte* gpReconCr;  // reconstructed ChromaR frame

// for debug and reseach

int* gpSADLu; // SAD for every luma block matrix
int* gpSADCb; // SAD for every ChromaB block matrix
int* gpSADCr; // SAD for every ChromaR block matrix

// histograms

int gRunsHist[MAX_BLK_DIM*MAX_BLK_DIM];
int gFreqHist[1024];

//=======================================================================================
void Test()
{
 int qp=4;
 short a = (((1<<qp)+1)/2) - 1;
 printf("test: %d   %d\n", qp, a);
}

//=======================================================================================
// Init parameters
void InitParams()
{
 gQPLu = 31; // luma quantization parameter

 gNBLu = 8; // luma macroblock size
 gNTLu = 8; // luma block transform dim

 gLuW  = gW;  // width  of Luma frame
 gLuH  = gH;  // height of Luma frame

 gQPCb = 31; // chromaB quantization param

 gNBCb = 8; // ChromaB macroblock size
 gNTCb = 8; // ChromaB block transform dim

 gCbW  = gW ;  // width  of ChromaB frame
 gCbH  = gH ;  // height of ChromaB frame

 gQPCr = 31; // ChromaR quantization parameter

 gNBCr = 8; // ChromaR macroblock size
 gNTCr = 8; // ChromaR block transform dim

 gCrW  = gW ;  // width  of ChromaR frame
 gCrH  = gH ;  // height of ChromaR frame
}

//=======================================================================================
void AllocateMem()
{
 gpRed   = malloc(gW*gH);
 gpGreen = malloc(gW*gH);
 gpBlue  = malloc(gW*gH);

 gpLu = malloc(gLuW * gLuH);  // luma frame
 gpCb = malloc(gCbW * gCbH);  // chroma Cb frame
 gpCr = malloc(gCrW * gCrH);  // chroma Cr frame

 gpResidLu = malloc(2* gLuW * gLuH); // residual frame between initial and predicted Luma frames
 gpResidCb = malloc(2* gCbW * gCbH); // residual frame between initial and predicted ChromaB frames
 gpResidCr = malloc(2* gCrW * gCrH); // residual frame between initial and predicted ChromaR frames

 gpPredTypeLu = malloc(gLuW * gLuH / gNBLu / gNBLu); // best prediction mode for every block in Luma frame
 gpPredTypeCb = malloc(gCbW * gCbH / gNBCb / gNBCb); // best prediction mode for every block in ChromaB frame
 gpPredTypeCr = malloc(gCrW * gCrH / gNBCr / gNBCr); // best prediction mode for every block in ChromaR frame

 gpPredLu = malloc(gLuW * gLuH);   // prediction frame for Luma
 gpPredCb = malloc(gCbW * gCbH);   // prediction frame for ChromaB
 gpPredCr = malloc(gCrW * gCrH);   // prediction frame for ChromaR

 gpTQLu = malloc(2* gLuW * gLuH);  // transformed quantized residual Luma matrix
 gpTQCb = malloc(2* gCbW * gCbH);  // transformed quantized residual ChromaB matrix
 gpTQCr = malloc(2* gCrW * gCrH);  // transformed quantized residual ChromaR matrix

 gpReordLu = malloc(2* gLuW * gLuH); // reordering transformed quantized residual Luma matrix
 gpReordCb = malloc(2* gCbW * gCbH); // reordering transformed quantized residual ChromaB matrix
 gpReordCr = malloc(2* gCrW * gCrH); // reordering transformed quantized residual ChromaR matrix

 gpReconLu = malloc(gLuW * gLuH);  // reconstructed Luma frame
 gpReconCb = malloc(gCbW * gCbH);  // reconstructed ChromaB frame
 gpReconCr = malloc(gCrW * gCrH);  // reconstructed ChromaR frame

 gpSADLu = malloc(4* gLuW * gLuH / gNBLu / gNBLu); // SAD for every luma block matrix
 gpSADCb = malloc(4* gCbW * gCbH / gNBCb / gNBCb); // SAD for every ChromaB block matrix
 gpSADCr = malloc(4* gCrW * gCrH / gNBCr / gNBCr); // SAD for every ChromaR block matrix
}

//=======================================================================================
void FreeMem()
{
 free(gpSADCr);
 free(gpSADCb);
 free(gpSADLu);

 free(gpReconCr);
 free(gpReconCb);
 free(gpReconLu);

 free(gpReordCr);
 free(gpReordCb);
 free(gpReordLu);

 free(gpTQCr);
 free(gpTQCb);
 free(gpTQLu);

 free(gpPredLu);
 free(gpPredCr);
 free(gpPredCb);

 free(gpPredTypeCr);
 free(gpPredTypeCb);
 free(gpPredTypeLu);

 free(gpResidCr);
 free(gpResidCb);
 free(gpResidLu);

 free(gpCr);
 free(gpCb);
 free(gpLu);

 free(gpBlue);
 free(gpGreen);
 free(gpRed);
 free(gpInp);
}

//=======================================================================================
// transform RGB to LuCbCr
void RGBToLuCbCrTransform()
{

 GetComponentsFromRGB(gpInp, gpRed, gpGreen, gpBlue, gW, gH);

 TransformRGBtoLuCbCr(gpRed, gpGreen, gpBlue, gpLu, gpCb, gpCr, gW, gH);

 WritePGM(gpLu, gW, gH, "test_lu.pgm");
 WritePGM(gpCb, gW, gH, "test_cb.pgm");
 WritePGM(gpCr, gW, gH, "test_cr.pgm");
}

//=======================================================================================
// transform LuCbCr to RGB
void LuCbCrToRGBTransform()
{
 TransformLuCbCrToRGB(gpReconLu, gpReconCb, gpReconCr, gpRed,  gpGreen, gpBlue, gW, gH);

 MakeRGBfromComponents(gpInp, gpRed, gpGreen, gpBlue, gW, gH);

 PpmImageFileWrite("test.ppm", gpInp, gW, gH);
}

//=======================================================================================
// input:
// pMat - transformed and quantized mairix matrix 
// output:
// pVec - reordered vector
void Reordering(short* pMat, int nWidth, int nHeight, int blkDim, short* pVec)
{
 int i;

 int nBWidth = nWidth/blkDim;
 int nBHeight= nHeight/blkDim;
 int nBx, nBy;
 int xCoord, yCoord;
 short blk[MAX_BLK_DIM * MAX_BLK_DIM];
 
 int blkSize = blkDim*blkDim;
 
 int shift = 0;

 // reordering
 for(nBy=0; nBy<nBHeight; nBy++)
 {
 	 yCoord = nBy * blkDim;
 	 for(nBx=0; nBx<nBWidth; nBx++)
 	 {
 	 	 xCoord = nBx * blkDim;
 	 	 GetBlockShort(pMat, nWidth, nHeight, xCoord, yCoord, blk, blkDim, blkDim);
     for(i=0; i<blkSize; i++)
     {
     	 pVec[shift + i] = blk[reord8[i]];
     }
     shift += blkSize;
 	 }
 }
}

//=======================================================================================
void HistogramsCalcululation(int nWidth, int nHeight, short* pVec)
{
 int i;

 int nCount=0;
 int flag;
 
 double fEnt;
 int size=0;

 for(i=0; i<1024; i++)
 {
 	gFreqHist[i]=0;
 	gRunsHist[i]=0;
 }
 
 // histogram of frec
 for(i=0; i< nWidth*nHeight; i++)
 {
 	if(pVec[i]!=0) gFreqHist[pVec[i]+512]++;
 }

 // histogram of runs
 flag = 0;
 
 for(i=0; i< nWidth*nHeight; i++)
 {
 	if(pVec[i]==0)
 	{
 	  flag++;
 	}
 	else
 	{
 		gRunsHist[flag]++;
 		flag = 0;
 	}
 }

fprintf(gFp,"runs:\n");
fEnt = EntropyCalculation(150, gRunsHist);
printf("RunsEntropy=%f\n",fEnt);
nCount=0;
for(i=0; i<MAX_BLK_DIM*MAX_BLK_DIM; i++)
{
  nCount+=gRunsHist[i];
  if(gRunsHist[i]>0) fprintf(gFp,"%d   %d\n", i, gRunsHist[i]);
}
size = fEnt*nCount;
printf("nRuns=%d\n", nCount);

fprintf(gFp,"cos amplitudes:\n");
fEnt = EntropyCalculation(900, gFreqHist);
printf("CoeffsEntropy=%f\n",fEnt);
nCount=0;
for(i=300;i<700; i++)
{
  nCount+=gFreqHist[i];
  if(gFreqHist[i]>0) fprintf(gFp, "%d   %d\n", i-512, gFreqHist[i]);
}
printf("nCoefs=%d\n", nCount);
size = size + fEnt*nCount;

size = size/8;
printf("size=%d\n",size);

gnResultSize += size;
}

//=======================================================================================
// input:
// pVec - reordered vector
// output:
// nResSize - size of encoded result
// pRes - vector with encoding result
void ArithmeticEncoding(short* pMat, int nWidth, int nHeight, int* nResSize, char* pRes)
{
}

//=======================================================================================
//
// vik-file structure:
//
//VIKA1_
//V0.1
//#comments...
// 4 bytes - luma frame width
// 4 bytes - luma frame height
// 1 byte  - luma quantization parameter
// 1 byte  - luma macroblock size
// 1 byte  - luma block transform dim
// 1 byte -  Id of transform matrix for luma
// 1 byte  - luma type of reordering transform
// 1 byte  - luma quantization type
// 4 bytes - chromaB frame width
// 4 bytes - chromaB frame height
// 1 byte  - chromaB quantization parameter
// 1 byte  - chromaB macroblock size
// 1 byte  - chromaB block transform dim
// 1 byte -  Id of transform matrix for chromaB
// 1 byte  - chromaB type of reordering transform
// 1 byte  - chromaB quantization type
// 4 bytes - chromaR frame width
// 4 bytes - chromaR frame height
// 1 byte  - chromaR quantization parameter
// 1 byte  - chromaR macroblock size
// 1 byte  - chromaR block transform dim
// 1 byte -  Id of transform matrix for chromaR
// 1 byte  - chromaR type of reordering transform
// 1 byte  - chromaR quantization type
//
//
// coding tables for luma
//
// coded matrix of prediction modes for luma
//
// coded luma frame
//
// coding tables for chromaB
//
// coded matrix of prediction modes for chromaB
//
// coded chromaB frame
//
// coding tables for chromaR
//
// coded matrix of prediction modes for chromaR
//
// coded chromaR frame


//=======================================================================================
// vik-file open and write header
void VikFileOpen(char* filename, FILE** pFp)
{
}

//=======================================================================================
void VikFileSaving(fp, gLuW, gLuH, gNBLu, gpPredTypeLu, nResSize, pArray)
{
}

//=======================================================================================
// types of prediction modes:
// 0 - mean value of borders
// 1 - left
// 2 - up
// 3 - diagonal down right
// 4 - diagonal down left
//
// input:
//   blk      - block with part of frame
//   blkSize  -
//   nBx, nBy -
//   pRefImg  -
//
// return:
//   blk - predicted block
//   *pPredictType
//   SAD
int BestPrediction(byte* blk, byte blkDim, int nBx, int nBy,
                   byte* pRefImg, int nWidth, byte* pPredictType)
{
 int i,j;
 int mean;
 
 int nBlkWidth = nWidth / blkDim;

 int minSAD = 2000000000;
 byte bestMode;
 
 byte predBlk[MAX_BLK_DIM*MAX_BLK_DIM];
 byte predBlkBest[MAX_BLK_DIM*MAX_BLK_DIM];

 byte bordL[MAX_BLK_DIM];  // left border
 byte bordU[MAX_BLK_DIM];  // up border
 byte bordUR[MAX_BLK_DIM]; // up right border

 int xCoord = nBx*blkDim;
 int yCoord = nBy*blkDim;
 
 int blkSize= blkDim*blkDim;

 int nSAD = 0;

 // interior of frame
 if(nBx>0 && nBy>0)
 {
 	 // get left border
   GetBlockByte(pRefImg, gLuW, gLuH, xCoord-blkDim, yCoord, predBlk, blkDim, blkDim);
   for(i=0; i<blkDim; i++) bordL[i] = predBlk[i*blkDim + blkDim - 1];

   // get upper border
   GetBlockByte(pRefImg, gLuW, gLuH, xCoord, yCoord-blkDim, predBlk, blkDim, blkDim);
   for(i=0; i<blkDim; i++) bordU[i] = predBlk[blkDim*(blkDim-1) + i];

   // mode 0
   mean = 0;
   for(i=0; i<blkDim; i++) mean += bordL[i];
   for(i=0; i<blkDim; i++) mean += bordU[i];
   mean = (mean + blkDim -1) /2 /blkDim;

   for(i=0; i<blkSize; i++) predBlk[i] = mean;

   nSAD = CalcSAD(blk, predBlk, blkDim);
   if(nSAD<minSAD)
   {
   	 bestMode = 0;
     minSAD = nSAD;
     for(i=0; i<blkSize; i++) predBlkBest[i] = predBlk[i];
   }

   // mode 1
   for(i=0; i<blkDim; i++)
   {
   	for(j=0; j<blkDim; j++)
   	{
   		predBlk[i*blkDim+j] = bordL[i];
   	}
   }

   nSAD = CalcSAD(blk, predBlk, blkDim);
   if(nSAD<minSAD)
   {
   	 bestMode = 1;
     minSAD = nSAD;
     for(i=0; i<blkSize; i++) predBlkBest[i] = predBlk[i];
   }

   // mode 2
   for(i=0; i<blkDim; i++)
   {
   	 for(j=0; j<blkDim; j++)
   	 {
   		 predBlk[i*blkDim+j] = bordU[j];
   	 }
   }

   nSAD = CalcSAD(blk, predBlk, blkDim);
   if(nSAD<minSAD)
   {
   	 bestMode = 2;
     minSAD = nSAD;
     for(i=0; i<blkSize; i++) predBlkBest[i] = predBlk[i];
   }

   // mode 3

   for(i=1; i<blkDim; i++)
   {
   	 for(j=0; j<(blkDim-i); j++)
   	 {
   	 	 predBlk[j+(i+j)*blkDim] = bordL[i-1];
   	 	 predBlk[i+j+(j)*blkDim] = bordU[i-1];
   	 }
   }
   for(i=0; i<blkDim; i++) predBlk[i+i*blkDim] = (bordL[0] + bordU[0] + 1)/2;

   nSAD = CalcSAD(blk, predBlk, blkDim);
   if(nSAD<minSAD)
   {
   	 bestMode = 3;
     minSAD = nSAD;
     for(i=0; i<blkSize; i++) predBlkBest[i] = predBlk[i];
   }

   // mode 4
   if(nBx < (nBlkWidth-1))
   {
     // get upper border
     GetBlockByte(pRefImg, gLuW, gLuH, xCoord+blkDim, yCoord-blkDim, predBlk, blkDim, blkDim);
     for(i=0; i<blkDim; i++) bordUR[i] = predBlk[blkDim*(blkDim-1) + i];

     for(i=0; i<(blkDim-1); i++)
     {
     	 for(j=0; j<(i+1); j++)
     	 {
   	 	   predBlk[(i-j)+(j)*blkDim] = bordU[i+1];
    	 }
     }

     for(i=0; i<blkDim; i++)
     {
   	   for(j=i; j<blkDim; j++)
   	   {
   	 	   predBlk[j+(blkDim-1 + i - j)*blkDim] = bordUR[i];
   	   }
     }

     nSAD = CalcSAD(blk, predBlk, blkDim);
     if(nSAD<minSAD)
     {
   	   bestMode = 4;
       minSAD = nSAD;
       for(i=0; i<blkSize; i++) predBlkBest[i] = predBlk[i];
     }
   }
 }
 // top border of frame
 else if(nBx>0 && nBy==0)
 {
 	 // get left border
   GetBlockByte(pRefImg, gLuW, gLuH, xCoord-blkDim, yCoord, predBlk, blkDim, blkDim);
   for(i=0; i<blkDim; i++) bordL[i] = predBlk[i*blkDim + blkDim - 1];

   // mode 0
   mean = 0;
   for(i=0; i<blkDim; i++) mean += bordL[i];
   mean = (mean + blkDim /2) /blkDim;

   for(i=0; i<blkSize; i++) predBlk[i] = mean;

   nSAD = CalcSAD(blk, predBlk, blkDim);
   if(nSAD<minSAD)
   {
   	 bestMode = 0;
     minSAD = nSAD;
     for(i=0; i<blkSize; i++) predBlkBest[i] = predBlk[i];
   }

   // mode 1
   for(i=0; i<blkDim; i++)
   {
   	for(j=0; j<blkDim; j++)
   	{
   		predBlk[i*blkDim+j] = bordL[i];
   	}
   }

   nSAD = CalcSAD(blk, predBlk, blkDim);
   if(nSAD<minSAD)
   {
   	 bestMode = 1;
     minSAD = nSAD;
     for(i=0; i<blkSize; i++) predBlkBest[i] = predBlk[i];
   }

 }
 // left border of frame
 else if(nBx==0 && nBy>0)
 {
   // get upper border
   GetBlockByte(pRefImg, gLuW, gLuH, xCoord, yCoord-blkDim, predBlk, blkDim, blkDim);
   for(i=0; i<blkDim; i++) bordU[i] = predBlk[blkDim*(blkDim-1) + i];

   // mode 0
   mean = 0;
   for(i=0; i<blkDim; i++) mean += bordU[i];
   mean = (mean + blkDim /2) /blkDim;

   for(i=0; i<blkSize; i++) predBlk[i] = mean;

   nSAD = CalcSAD(blk, predBlk, blkDim);
   if(nSAD<minSAD)
   {
   	 bestMode = 0;
     minSAD = nSAD;
     for(i=0; i<blkSize; i++) predBlkBest[i] = predBlk[i];
   }

   // mode 2
   for(i=0; i<blkDim; i++)
   {
   	for(j=0; j<blkDim; j++)
   	{
   		predBlk[i*blkDim+j] = bordU[j];
   	}
   }

   nSAD = CalcSAD(blk, predBlk, blkDim);
   if(nSAD<minSAD)
   {
   	 bestMode = 2;
     minSAD = nSAD;
     for(i=0; i<blkSize; i++) predBlkBest[i] = predBlk[i];
   }

 }
 // left top block of frame
 else 
 { // nBx==0 && nBy==0
   bestMode = 0;
   mean = 127;
   for(i=0; i<blkSize; i++) predBlkBest[i] = mean;
   minSAD = CalcSAD(blk, predBlkBest, blkDim);
 }

 *pPredictType = bestMode;

 for(i=0; i<blkSize; i++) blk[i] = predBlkBest[i];

 return minSAD;
}

//=======================================================================================
// quantize transformed residual block
int Quantizer( byte QP, short* blk, int blkDim)
{
 int nn = blkDim * blkDim;
 int nZero=0;
 int i;
 
 int numer = gnQuantNumer[QP];
 int denom = gnQuantDenom[QP];

 int a = ((numer+1)/2) - 1;
 int r;

 if(QP>64) return 1;


 for(i=0; i<nn; i++) 
 {
 	r=abs(blk[i]);
 	r = (r * denom + a) / numer;
 	if(blk[i]>=0) blk[i] = r;
 	else blk[i] = -r;
 }

 return nZero;
}

//=======================================================================================
// dequantize
int Dequantizer( byte QP, short* blk, int blkDim)
{
 int nn = blkDim * blkDim;
 int nZero=0;
 int i;

 int numer = gnQuantNumer[QP];
 int denom = gnQuantDenom[QP];

 int a = ((denom+1)/2) - 1;

 for(i=0; i<nn; i++) 
 {
 	blk[i] = (blk[i] * numer + a) / denom;
 }

 return nZero;
}

//=======================================================================================
// input:
// gQPLu
// gpLu - uncompressed image
//
// output:
// gpTQLu - image after transform and quantization
// gpReconLu - reconstructed image
void BlockIntraLoopLuma()
{
  int i;
  
  short transform[MAX_BLK_DIM*MAX_BLK_DIM];
  short transformT[MAX_BLK_DIM*MAX_BLK_DIM];

  byte blk[MAX_BLK_DIM*MAX_BLK_DIM];
  byte predBlk[MAX_BLK_DIM*MAX_BLK_DIM];
  byte nPredictType;
  
  short residBlk[MAX_BLK_DIM*MAX_BLK_DIM];
  short transBlk[MAX_BLK_DIM*MAX_BLK_DIM];
  short tempBlk[MAX_BLK_DIM*MAX_BLK_DIM];

  byte reconBlk[MAX_BLK_DIM*MAX_BLK_DIM];

  int nSad;

  int nx,ny;    // block coordinates
  int nxx, nyy; // coordinates of top left corner of block

  int nBWidth  = gLuW / gNBLu;
  int nBHeight = gLuH / gNBLu;
  
  double fEntr;
  
  int hist[16];  //debug
  
  for(i=0; i<16; i++) 
  {
  	hist[i]=0;
  }

  MatrCopy(gNBLu, hadam8, transform);
  MatrTranspos(gNBLu, transform, transformT);

  for(ny=0; ny<nBHeight; ny++)
  {
  	nyy = gNBLu * ny;
    for(nx=0; nx<nBWidth; nx++)
    {
      nxx = gNBLu * nx;

      // get block from input frame
      GetBlockByte(gpLu, gLuW, gLuH, nxx, nyy, blk, gNBLu, gNBLu);
      for(i=0; i<gNBLu*gNBLu; i++) predBlk[i] = blk[i];

      // get best prediction for input block and save in gpPredLu frame
      nSad = BestPrediction(predBlk, gNBLu, nx, ny, gpReconLu, gLuW, &nPredictType);

      gpPredTypeLu[ny*nBWidth + nx] = nPredictType;
      gpSADLu[ny*nBWidth + nx]      = nSad;

      PutBlockByte(gpPredLu, gLuW, gLuH, nxx, nyy, predBlk, gNBLu, gNBLu);
      hist[nPredictType]++;
      
      // difference (block-predictedBlock) = residual block
      MatrSubtract(gNBLu, blk, predBlk, residBlk);
if(nx==ny) fprintf(gFp, "##############################################################\n");

      // transform of residual block
      MatrMult(gNBLu, transform, residBlk, tempBlk);
      MatrMult(gNBLu, tempBlk, transformT, transBlk);
      MatrDivideConst(gNBLu, transBlk, 8, 1);

      // quantize transformed residual block
      Quantizer( gQPLu, transBlk, gNBLu);
if(nx==ny) FileMatrPrintShort(gFp,gNBLu, transBlk);

      // save in frame gpTQLu
      PutBlockShort(gpTQLu, gLuW, gLuH, nxx, nyy, transBlk, gNBLu, gNBLu);

      // dequantize
      Dequantizer( gQPLu, transBlk, gNBLu);

      // reverse transformation
      MatrMult(gNBLu, transformT, transBlk, tempBlk);
      MatrMult(gNBLu, tempBlk, transform, residBlk);
      MatrDivideConst(gNBLu, residBlk, 8, 1);

      // reconstructed block = predicted block + restored residual
      MatrAdd(gNBLu, predBlk, residBlk, reconBlk);

      // save reconstructed block in frame gpReconL
      PutBlockByte(gpReconLu, gLuW, gLuH, nxx, nyy, reconBlk, gNBLu, gNBLu);
    }
  }
  for(i=0; i<5; i++) printf("%d   ",hist[i]);
  printf("\n");
  fEntr=EntropyCalculation(6, hist);
  printf("PredictionTypeEnthropy= %f \n", fEntr);

}

//=======================================================================================
// input:
// gQPCb
// gpCb - uncompressed image
//
// output:
// gpTQCb - image after transform and quantization
// gpReconCb - reconstructed image
void BlockIntraLoopChromaB()
{
  int i;
  
  short transform[MAX_BLK_DIM*MAX_BLK_DIM];
  short transformT[MAX_BLK_DIM*MAX_BLK_DIM];

  byte blk[MAX_BLK_DIM*MAX_BLK_DIM];
  byte predBlk[MAX_BLK_DIM*MAX_BLK_DIM];
  byte nPredictType;
  
  short residBlk[MAX_BLK_DIM*MAX_BLK_DIM];
  short transBlk[MAX_BLK_DIM*MAX_BLK_DIM];
  short tempBlk[MAX_BLK_DIM*MAX_BLK_DIM];

  byte reconBlk[MAX_BLK_DIM*MAX_BLK_DIM];

  int nSad;

  int nx,ny;    // block coordinates
  int nxx, nyy; // coordinates of top left corner of block

  int nBWidth  = gCbW / gNBCb;
  int nBHeight = gCbH / gNBCb;
  
  double fEntr;
  
  int hist[16];  //debug
  
  for(i=0; i<16; i++) 
  {
  	hist[i]=0;
  }

  MatrCopy(gNBCb, hadam8, transform);
  MatrTranspos(gNBCb, transform, transformT);

  for(ny=0; ny<nBHeight; ny++)
  {
  	nyy = gNBCb * ny;
    for(nx=0; nx<nBWidth; nx++)
    {
      nxx = gNBCb * nx;

      // get block from input frame
      GetBlockByte(gpCb, gCbW, gCbH, nxx, nyy, blk, gNBCb, gNBCb);
      for(i=0; i<gNBCb*gNBCb; i++) predBlk[i] = blk[i];

      // get best prediction for input block and save in gpPredCb frame
      nSad = BestPrediction(predBlk, gNBCb, nx, ny, gpReconCb, gCbW, &nPredictType);

      gpPredTypeCb[ny*nBWidth + nx] = nPredictType;
      gpSADCb[ny*nBWidth + nx]      = nSad;

      PutBlockByte(gpPredCb, gCbW, gCbH, nxx, nyy, predBlk, gNBCb, gNBCb);
      hist[nPredictType]++;
      
      // difference (block-predictedBlock) = residual block
      MatrSubtract(gNBCb, blk, predBlk, residBlk);
if(nx==ny) fprintf(gFp, "##############################################################\n");

      // transform of residual block
      MatrMult(gNBCb, transform, residBlk, tempBlk);
      MatrMult(gNBCb, tempBlk, transformT, transBlk);
      MatrDivideConst(gNBCb, transBlk, 8, 1);

      // quantize transformed residual block
      Quantizer( gQPCb, transBlk, gNBCb);
if(nx==ny) FileMatrPrintShort(gFp,gNBCb, transBlk);

      // save in frame gpTQCb
      PutBlockShort(gpTQCb, gCbW, gCbH, nxx, nyy, transBlk, gNBCb, gNBCb);

      // dequantize
      Dequantizer( gQPCb, transBlk, gNBCb);

      // reverse transformation
      MatrMult(gNBCb, transformT, transBlk, tempBlk);
      MatrMult(gNBCb, tempBlk, transform, residBlk);
      MatrDivideConst(gNBCb, residBlk, 8, 1);

      // reconstructed block = predicted block + restored residual
      MatrAdd(gNBCb, predBlk, residBlk, reconBlk);

      // save reconstructed block in frame gpReconL
      PutBlockByte(gpReconCb, gCbW, gCbH, nxx, nyy, reconBlk, gNBCb, gNBCb);
    }
  }
  for(i=0; i<5; i++) printf("%d   ",hist[i]);
  printf("\n");
  fEntr=EntropyCalculation(6, hist);
  printf("PredictionTypeEnthropy= %f \n", fEntr);

}

//=======================================================================================
// input:
// gQPCr
// gpCr - uncompressed image
//
// output:
// gpTQCr - image after transform and quantization
// gpReconCr - reconstructed image
void BlockIntraLoopChromaR()
{
  int i;
  
  short transform[MAX_BLK_DIM*MAX_BLK_DIM];
  short transformT[MAX_BLK_DIM*MAX_BLK_DIM];

  byte blk[MAX_BLK_DIM*MAX_BLK_DIM];
  byte predBlk[MAX_BLK_DIM*MAX_BLK_DIM];
  byte nPredictType;
  
  short residBlk[MAX_BLK_DIM*MAX_BLK_DIM];
  short transBlk[MAX_BLK_DIM*MAX_BLK_DIM];
  short tempBlk[MAX_BLK_DIM*MAX_BLK_DIM];

  byte reconBlk[MAX_BLK_DIM*MAX_BLK_DIM];

  int nSad;

  int nx,ny;    // block coordinates
  int nxx, nyy; // coordinates of top left corner of block

  int nBWidth  = gCrW / gNBCr;
  int nBHeight = gCrH / gNBCr;
  
  double fEntr;
  
  int hist[16];  //debug
  
  for(i=0; i<16; i++) 
  {
  	hist[i]=0;
  }

  MatrCopy(gNBCr, hadam8, transform);
  MatrTranspos(gNBCr, transform, transformT);

  for(ny=0; ny<nBHeight; ny++)
  {
  	nyy = gNBCr * ny;
    for(nx=0; nx<nBWidth; nx++)
    {
      nxx = gNBCr * nx;

      // get block from input frame
      GetBlockByte(gpCr, gCrW, gCrH, nxx, nyy, blk, gNBCr, gNBCr);
      for(i=0; i<gNBCr*gNBCr; i++) predBlk[i] = blk[i];

      // get best prediction for input block and save in gpPredCr frame
      nSad = BestPrediction(predBlk, gNBCr, nx, ny, gpReconCr, gCrW, &nPredictType);

      gpPredTypeCr[ny*nBWidth + nx] = nPredictType;
      gpSADCr[ny*nBWidth + nx]      = nSad;

      PutBlockByte(gpPredCr, gCrW, gCrH, nxx, nyy, predBlk, gNBCr, gNBCr);
      hist[nPredictType]++;
      
      // difference (block-predictedBlock) = residual block
      MatrSubtract(gNBCr, blk, predBlk, residBlk);
if(nx==ny) fprintf(gFp, "##############################################################\n");

      // transform of residual block
      MatrMult(gNBCr, transform, residBlk, tempBlk);
      MatrMult(gNBCr, tempBlk, transformT, transBlk);
      MatrDivideConst(gNBCr, transBlk, 8, 1);

      // quantize transformed residual block
      Quantizer( gQPCr, transBlk, gNBCr);
if(nx==ny) FileMatrPrintShort(gFp,gNBCr, transBlk);

      // save in frame gpTQCr
      PutBlockShort(gpTQCr, gCrW, gCrH, nxx, nyy, transBlk, gNBCr, gNBCr);

      // dequantize
      Dequantizer( gQPCr, transBlk, gNBCr);

      // reverse transformation
      MatrMult(gNBCr, transformT, transBlk, tempBlk);
      MatrMult(gNBCr, tempBlk, transform, residBlk);
      MatrDivideConst(gNBCr, residBlk, 8, 1);

      // reconstructed block = predicted block + restored residual
      MatrAdd(gNBCr, predBlk, residBlk, reconBlk);

      // save reconstructed block in frame gpReconL
      PutBlockByte(gpReconCr, gCrW, gCrH, nxx, nyy, reconBlk, gNBCr, gNBCr);
    }
  }
  for(i=0; i<5; i++) printf("%d   ",hist[i]);
  printf("\n");
  fEntr=EntropyCalculation(6, hist);
  printf("PredictionTypeEnthropy= %f \n", fEntr);

}

//=======================================================================================
int main(int nArgs,char** ppArgs)
{
 FILE* fp; // for vik-file

 char filename[256];
 
 int nResSize;
 char* pArray;

//Test();

 pArray = malloc(1000000);

 puts("======================================================");
 puts("Vika still image encoder:");
 gFp = fopen("log.txt", "w");

 switch(nArgs)
 {
 	case 1:   // help
 	  break;
 	case 2:
 	  strcpy(filename, ppArgs[1]);
 	  break;
 	case 3:
 	  strcpy(filename, ppArgs[1]);
 	  break;
 	case 4:
 	  strcpy(filename, ppArgs[1]);
 	  break;
 	case 5:
 	  strcpy(filename, ppArgs[1]);
 	  break;
 	default:
 	  break;
 }
// open and read image file with memory allocation for it
 if( PpmImageFileRead(filename, &gpInp, &gW, &gH) )
 {
   puts("image file not found!");
   return 1;
 }
 
 // Init parameters
 InitParams();
 
 // command line analysis
 switch(nArgs)
 {
 	case 3:
 	  gQPLu = atoi(ppArgs[2]);
 	  break;
 	case 4:
 	  gQPLu = atoi(ppArgs[2]);
 	  gQPCb = atoi(ppArgs[3]);
 	  break;
 	case 5:
 	  gQPLu = atoi(ppArgs[2]);
 	  gQPCb = atoi(ppArgs[3]);
 	  gQPCr = atoi(ppArgs[4]);
 	  break;
 	default:
 	  break;
 }

 // allocate memory for global frames
 AllocateMem();

 // transform RGB to LuCbCr
 RGBToLuCbCrTransform();

 // write header of vik-file
 VikFileOpen("result.vik", &fp);

 //----------------------------------------------------------------------------
 printf("------------------------\nLuma:\n");
 BlockIntraLoopLuma();

 WritePGM(gpPredLu,  gLuW, gLuH, "test_pred_lu.pgm");
 WritePGM(gpReconLu, gLuW, gLuH, "test_recon_lu.pgm");

 // reordering, , and saving in vik-file
 Reordering(gpTQLu, gLuW, gLuH, gNBLu, gpReordLu);
 HistogramsCalcululation(gLuW, gLuH, gpReordLu);
 ArithmeticEncoding(gpReordLu, gLuW, gLuH, &nResSize, pArray);
 VikFileSaving(fp, gLuW, gLuH, gNBLu, gpPredTypeLu, nResSize, pArray);
 
 //----------------------------------------------------------------------------
 printf("------------------------\nCromaB:\n");
 BlockIntraLoopChromaB();

 WritePGM(gpPredCb,  gCbW, gCbH, "test_pred_Cb.pgm");
 WritePGM(gpReconCb, gCbW, gCbH, "test_recon_Cb.pgm");

 // reordering, enthropy coding, and saving in vik-file
 Reordering(gpTQCb, gCbW, gCbH, gNBCb, gpReordCb);
 HistogramsCalcululation(gCbW, gCbH, gpReordCb);
 ArithmeticEncoding(gpReordCb, gCbW, gCbH, &nResSize, pArray);
 VikFileSaving(fp, gCbW, gCbH, gNBCb, gpPredTypeCb, nResSize, pArray);

 //----------------------------------------------------------------------------
 printf("------------------------\ncromaR:\n");
 BlockIntraLoopChromaR();

 WritePGM(gpPredCr,  gCrW, gCrH, "test_pred_Cr.pgm");
 WritePGM(gpReconCr, gCrW, gCrH, "test_recon_Cr.pgm");

 // reordering, enthropy coding, and saving in vik-file
 Reordering(gpTQCr, gCrW, gCrH, gNBCr, gpReordCr);
 HistogramsCalcululation(gCrW, gCrH, gpReordCr);
 ArithmeticEncoding(gpReordCr, gCrW, gCrH, &nResSize, pArray);
 VikFileSaving(fp, gCrW, gCrH, gNBCr, gpPredTypeCr, nResSize, pArray);

 // Decoded RGB result
 LuCbCrToRGBTransform();

 // decoded result
 printf("------------------------\ntotal size= %d\n",gnResultSize);
 printf("Luma    SNR=%f\n", CalculateSNR(gpReconLu, gpLu, gLuW, gLuH));
 printf("CromaCb SNR=%f\n", CalculateSNR(gpReconCb, gpCb, gCbW, gCbH));
 printf("CromaCr SNR=%f\n", CalculateSNR(gpReconCr, gpCr, gCrW, gCrH));

 //----------------------------------------------------------------------------
 //free all allocated memory
 FreeMem();
 free(pArray);

 fclose(gFp);

 return 0;
}