#pragma managed

#include "stdafx.h"
#include "matrixOperator.h"
//#include "fft2.h"

#using <mscorlib.dll>
#using ".\dll\Exocortex.DSP.v1.dll"

using namespace std;
using namespace System;
//using namespace System::Math;
using namespace Exocortex;
using namespace Exocortex::DSP;

const int GaborWidth = 101, GaborHeight = 101;
//#pragma push_macro("new")
//#undef new	
//Complex KernelFFT2[];//=new Complex[256*256];
//GaborResult gaborResult;
//#pragma pop_macro("new")

//FFTGaborTransform
//void FFTGaborTransform(HBIO hBIO);
void PrepareKernel(Complex KernelFFT2[],int Orientation, int Frequency);
void CalculateKernel(Complex KernelFFT2[],int Orientation, int Frequency);
double KernelRealPart(int x, int y, int Orientation, int Frequency);
double KernelImgPart(int x, int y, int Orientation, int Frequency);
Matrix DoTransform(Matrix image, Complex KernelFFT2[],int Orientation, int Frequency, GaborResult result);


Matrix Image2Matrix(const HBIO &hBIO)
{	
	LPSTR lpBIO,lpBIOBits;
	lpBIO = (LPSTR) ::GlobalLock((HGLOBAL) hBIO);
	lpBIOBits=::FindDIBBits(lpBIO);	
	int nWidth=DIBWidth(lpBIO);
	int nHeight=DIBHeight(lpBIO);
	DWORD dwBitsSize=nWidth*nHeight;

	Matrix m_matBits(nHeight,nWidth);
	int wPos,hPos;
	for(hPos=0;hPos<nHeight;hPos++)
	{
		for(wPos=0;wPos<nWidth;wPos++)
		{
			m_matBits(hPos,wPos)=*(lpBIOBits+wPos+hPos*nWidth);
			//*
			if(m_matBits(hPos,wPos)<0)
			{
				m_matBits(hPos,wPos)+=256;
			}
			//*/
		}
	}
	::GlobalUnlock((HGLOBAL) hBIO);
	return m_matBits;	
}

HBIO Matrix2Image(const Matrix &m_matBits)
{
	unsigned int nWidth,nHeight;
	nWidth=m_matBits.ColNo();
	nHeight=m_matBits.RowNo();
	DWORD dwBitsSize=nWidth*nHeight;

	BIOFILEHEADER fileHeader;			
	fileHeader.flag=5;
	fileHeader.Width=nWidth;
	fileHeader.Height=nHeight;
	fileHeader.GrayNum=255;

	// 為BIO分配內存
	HBIO hBIO=NULL;
	LPSTR pBIO;
	hBIO = (HBIO) ::GlobalAlloc(GMEM_MOVEABLE | GMEM_ZEROINIT, dwBitsSize+FileHeadLen);
	if (hBIO == 0)
	{	
		return NULL;
	}	
	
	pBIO= (LPSTR)::GlobalLock((HGLOBAL) hBIO);
	memcpy(pBIO,(char*)&fileHeader,FileHeadLen);

	int wPos,hPos;
	for(hPos=0;hPos<nHeight;hPos++)
		for(wPos=0;wPos<nWidth;wPos++)
			*(pBIO+wPos+hPos*nWidth+FileHeadLen)=m_matBits(hPos,wPos);

	::GlobalUnlock((HGLOBAL) hBIO);

	return hBIO;
}

Matrix Buffer2Matrix(const HBIO &hBIO)
{	
	LPSTR lpBIO,lpBIOBits;
	lpBIO = (LPSTR) ::GlobalLock((HGLOBAL) hBIO);
	lpBIOBits=::FindDIBBits(lpBIO);
	BIOFILEHEADER fileH;
	//FILEDATA fileData;
	memcpy((char*)&fileH,lpBIO,FileHeadLen);
	int nWidth=fileH.Width;
	int nHeight=fileH.Height;
	int dbsize=sizeof(double);
	DWORD dwBitsSize=nWidth*nHeight*dbsize;

	Matrix m_matBits(nHeight,nWidth);
	int wPos,hPos;
	for(hPos=0;hPos<nHeight;hPos++)
	{
		for(wPos=0;wPos<nWidth;wPos++)
		{
			m_matBits(hPos,wPos)=double(*((double*)(lpBIOBits+dbsize*(wPos+hPos*nWidth))));
			
		}
	}
	
	::GlobalUnlock((HGLOBAL) hBIO);
	return m_matBits;	
}

HBIO Matrix2Buffer(const Matrix &m_matBits)
{
	unsigned int nWidth,nHeight;
	nWidth=m_matBits.ColNo();
	nHeight=m_matBits.RowNo();
	int dbsize=sizeof(double);
	DWORD dwBitsSize=nWidth*nHeight*dbsize;

	BIOFILEHEADER fileHeader;			
	fileHeader.flag=5;
	fileHeader.Width=nWidth;
	fileHeader.Height=nHeight;
	fileHeader.GrayNum=255;

	// 為BIO分配內存
	HBIO hBIO=NULL;
	LPSTR pBIO;
	hBIO = (HBIO) ::GlobalAlloc(GMEM_MOVEABLE | GMEM_ZEROINIT, dwBitsSize+FileHeadLen);
	if (hBIO == 0)
	{	
		return NULL;
	}	
	
	pBIO= (LPSTR)::GlobalLock((HGLOBAL) hBIO);
	memcpy(pBIO,(char*)&fileHeader,FileHeadLen);

	int wPos,hPos;
	for(hPos=0;hPos<nHeight;hPos++)
		for(wPos=0;wPos<nWidth;wPos++)
			*((double*)(pBIO+FileHeadLen+dbsize*(wPos+hPos*nWidth)))=m_matBits(hPos,wPos);

	::GlobalUnlock((HGLOBAL) hBIO);

	return hBIO;
}


//注意此處未作調試
Matrix Convolution2D(Matrix mImage,Matrix mKernel)
{
	Matrix result(mImage);
	int nHeight=mImage.RowNo();
	int nWidth=mImage.ColNo();
	int sizeH=mKernel.RowNo();
	int sizeW=mKernel.ColNo();

	int kPosX,kPosY;
	double sum=0;
	//Matrix exatrctM(sizeH,sizeW);

	int resultPosX=sizeW/2,resultPosY=sizeH/2;

	for(resultPosX=sizeH/2;resultPosX<nHeight-sizeH/2;resultPosX++)
	{
		for(resultPosY=sizeW/2;resultPosY<nWidth-sizeW/2;resultPosY++)
		{

			//convolution
			for(kPosX=0;kPosX<sizeH;kPosX++)
			{
				for(kPosY=0;kPosY<sizeW;kPosY++)
				{
					sum=sum+mKernel(kPosX,kPosY)*mImage(resultPosX+kPosX-sizeW/2,resultPosY+kPosY-sizeH/2);
				}
				
			}

			result(resultPosX,resultPosY)=sum;//(int)
			sum=0;			
		}
	}
	
	return result;
}

void printMatrix(Matrix &m,HDC hdc)
{
	char str[10];
	strncpy(str,"\0",10);
	for(int i=0;i<m.RowNo();i++)
	{
		for(int j=0;j<m.ColNo();j++)
		{
			sprintf(str,"%4.2f ",(float)m(i,j));
			::TextOut(hdc,384*2+40*j,+i*20,str,strlen(str));
		}
	}
}


void SaveMatrix(Matrix &m,string fstr)
{	
	string sPath;
	TCHAR lpBuffer[MAX_PATH];
	GetModuleFileName(NULL,lpBuffer,MAX_PATH);
	sPath=lpBuffer;
	
	int nPos;
	nPos=sPath.rfind ("\\");
	sPath=sPath.substr (0,nPos);
	string fil=sPath+"\\";
	fil+="result\\";
	fil+=fstr;
	std::ofstream f(fil.c_str());

	for(int y=0;y<m.RowNo();y++)
	{
		for(int x=0;x<m.ColNo();x++)
		{
			
			f<<m(y,x)<<" ";
		}
		f<<endl;
	}
	return;
}

Matrix LoadMatrix(int row,int col,string fstr)
{
	string sPath;
	TCHAR lpBuffer[MAX_PATH];
	GetModuleFileName(NULL,lpBuffer,MAX_PATH);
	sPath=lpBuffer;
	
	int nPos;
	nPos=sPath.rfind ("\\");
	sPath=sPath.substr (0,nPos);
	string fil=sPath+"\\";
	fil+="result\\";
	fil+=fstr;
	std::ifstream f(fil.c_str());
	Matrix m(row,col);

	for(int y=0;y<row;y++)
	{
		for(int x=0;x<col;x++)
		{
			
			f>>m(y,x);
		}
		//f>>endl;
	}
	return m;
}


//OneDimension Convolution
vector<float> Convolution(vector<float> h,vector<float> g)
{
	int M=h.size();
	int N=g.size();
	int resultSize=M+N-1;
	vector<float> result;
	float sum=0;
	int resultPos=0,mPos=0,nPos=0;

	while(resultPos<M+N)
	{
		for(mPos=0,nPos=resultPos-mPos;mPos<M,nPos>=0;mPos++,nPos--)
		{
			if(nPos<N)
				sum=sum+h[mPos]*g[nPos];
		}
		result.push_back(sum);
		resultPos++;
		sum=0;
		
	}
	return result;
}


//#pragma managed
//*
Matrix DoFFTConv2(Matrix data,Matrix kernel)
{
	int dWidth=data.ColNo();
	int dHeight=data.RowNo();
	int kWidth=kernel.ColNo();
	int kHeight=kernel.RowNo();

	int M=dHeight+kHeight+1;
	int N=dWidth+kWidth+1;
	int size=M*N;

	Complex Kernel[]= new Complex[size];
	Complex Data[]=new Complex[size];
	Complex Result[]=new Complex[size]; 

	for(int y=0;y<kHeight;y++)
		for(int x=0;x<kWidth;x++)
			Kernel[y*N+x].Re=kernel(y,x);

	for(int y=0;y<dHeight;y++)
		for(int x=0;x<dWidth;x++)
			Data[y*N+x].Re=data(y,x);

	Fourier::FFT2(Data, M, N, FourierDirection::Forward);
	Fourier::FFT2(Kernel, M, N, FourierDirection::Forward);
//*
	for(int i=0;i<size;i++)
		Result[i]=Data[i]*Kernel[i]/size;

	Fourier::FFT2(Result, M, N, FourierDirection::Backward);
//*/	
	Matrix result(dHeight+1,dWidth+1 );
	for(int y=0; y<dHeight+1; y++)//dHeight+1
    {
        for(int x=0; x<dWidth+1 ; x++) //dWidth+1        
		  result(y,x)=Result[y*N+x].Re;
        
     }
	return result;
	
}