// DCT2D.cpp: implementation of the CDCT2D class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "ArisImageRoutines.h"
#include "DCT2D.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CDCT2D::CDCT2D(unsigned int x, unsigned int y)
{
	pInput = new double*[x];
	for(unsigned int i = 0; i < x; i++)
		pInput[i] = new double[y];

	nX = x;
	nY = y;
}

CDCT2D::~CDCT2D()
{
	if(pInput != NULL)
	{
		for(unsigned  int i = 0; i < nX; i++)
		{
			delete[] pInput[i];
			pInput[i] = NULL;
		}
		pInput = NULL;
	}
}


double CDCT2D::GetInputAt(unsigned int x, unsigned int y)
{
	if(x >= 0 && x < nX && y >= 0 && y < nY)
		return pInput[x][y];
	else
	{
		#ifdef USE_MESSAGE_BOX_FOR_FAILURE		
			AfxMessageBox("Index out of range!");
		#endif
		return NULL;
	}
}

bool CDCT2D::SetInputAt(unsigned int x, unsigned int y, double real/* = 0 */, double imag/* = 0 */)
{
	if(x >= 0 && x < nX && y >= 0 && y < nY)
	{
		pInput[x][y] = (real);
		return true;
	}
	else
	{
		#ifdef USE_MESSAGE_BOX_FOR_FAILURE		
			AfxMessageBox("Index out of range!");
		#endif
		return false;
	}
}

void CDCT2D::ForwardTransform()
{
	DCT2D(0);	
}
	
void CDCT2D::ReverseTransform()
{
	DCT2D(1);
}


void CDCT2D::DCT2D(int inverse)
{
	int i,j,k;

	int m = nX;
	int n = nY;
	
	double **t,**y, u0,v0,u,v,w,wm,wn;
	
	t = new double*[m];
	for(i = 0; i < m; i++)
		t[i] = new double[n];

	y = new double*[m];
	for(i = 0; i < m; i++)
		y[i] = new double[n];

	k=m; 
	if (n>m) 
		k=n;
	
	wm=Pi/2/m; 
	wn=Pi/2/n;
	
	u0=sqrt(1.0/m); 
	u=sqrt(2.0/m);   
	v0=sqrt(1.0/n); 
	v=sqrt(2.0/n); 

// preparing for column transform matrix 
	for (i=0; i<m; i++)
		for (j=0; j<m; j++)
			if (inverse) 
				t[i][j]=cos( ((i<<1)+1)*j*wm );
			else t[i][j]=cos( ((j<<1)+1)*i*wm );

// matrix multiplication for column transform
	for (j=0; j<n; j++)                 
		for (i=0; i<m; i++)              
		{ 
			y[i][j]=0;
			for (k=0; k<m; k++)
			{
				if (inverse)
				{
					if (k) 
						w=u;  
					else w=u0; 
				}
				else w=1;
				y[i][j]+=w*t[i][k]*pInput[k][j];
			}
		}

// preparing for row transform matrix  
	for (i=0; i<n; i++)
		for (j=0; j<n; j++)
		if (inverse) t[i][j]=cos( ((j<<1)+1)*i*wn );
		else t[i][j]=cos(  ((i<<1)+1)*j*wn );

// matrix multiplication for row transform
	for (i=0; i<m; i++)                 
		for (j=0; j<n; j++)             
		{ 
			pInput[i][j]=0;
			for (k=0; k<n; k++)
			{	
				if (inverse)
				{
					if (k) 
						w=v;  
					else w=v0;
				}
				else w=1;

				pInput[i][j]+=w*y[i][k]*t[k][j];
			}
		}

	if(!inverse)
    { 
		pInput[0][0]*=u0*v0;
		for (i=1; i<m; i++) 
			pInput[i][0]*=u*v0;

		for (j=1; j<n; j++) 
			pInput[0][j]*=u0*v;

		for (i=1; i<m; i++)

		for (j=1; j<n; j++)
			pInput[i][j]*=u*v;
    }  

}