matrixoperator.cpp

Gabor Transformation based on a Exocortex dsp d

//Class better?
Matrix FFTGaborTransform(const HBIO &hBIO,int Orientation, int Frequency)
{
//	Console::WriteLine("   ==================== FFT Gabor Transform ===================");
//	DateTime begin = DateTime::Now;

	Matrix imageMatrix,resultMatrix;
	Complex KernelFFT2[]=new Complex[256*256];
	imageMatrix=Image2Matrix(hBIO);
	PrepareKernel(KernelFFT2,Orientation, Frequency);
	GaborResult result; //= new GaborResult;
	resultMatrix=DoTransform(imageMatrix, KernelFFT2,Orientation, Frequency, result);//image

//	Console::WriteLine();
//	Console::WriteLine(" Avg = {0,10:F3}\nDeta = {1,10:F3}", result.Avg, result.Deta);
//	DateTime end = DateTime::Now;
//	TimeSpan ts = end - begin;
//	Console::WriteLine("Total Time used: " + ts.TotalMilliseconds + " ms\r\n");
	return resultMatrix;
}

void PrepareKernel(Complex KernelFFT2[],int Orientation, int Frequency)
{
	
	CalculateKernel(KernelFFT2,Orientation, Frequency);

	Matrix dspkernelReal(256,256),dspkernelIm(256,256);
	for(int i=0;i<256;i++)
		for(int j=0;j<256;j++)
		{
			dspkernelReal(i,j)=KernelFFT2[i*256+j].Re;
			dspkernelIm(i,j)=KernelFFT2[i*256+j].Im;
		}
	SaveMatrix(dspkernelReal,"dspkernelReal");
	SaveMatrix(dspkernelIm,"dspkernelIm");


	Fourier::FFT2(KernelFFT2, 256, 256, FourierDirection::Forward );

	Matrix dspFFTkernelReal(256,256),dspFFTkernelIm(256,256);
	for(int i=0;i<256;i++)
		for(int j=0;j<256;j++)
		{
			dspFFTkernelReal(i,j)=KernelFFT2[i*256+j].Re;
			dspFFTkernelIm(i,j)=KernelFFT2[i*256+j].Im;
		}
	SaveMatrix(dspFFTkernelReal,"dspFFTkernelReal");
	SaveMatrix(dspFFTkernelIm,"dspFFTkernelIm");

}

void CalculateKernel(Complex KernelFFT2[],int Orientation, int Frequency)
{
	double real, img;
	//Complex KernelFFT2[];
	for(int x = -(GaborWidth-1)/2; x<(GaborWidth-1)/2+1; x++)
		for(int y = -(GaborHeight-1)/2; y<(GaborHeight-1)/2+1; y++)
		{
			real = KernelRealPart(x, y, Orientation, Frequency);
			img = KernelImgPart(x, y, Orientation, Frequency);
			KernelFFT2[(x+(GaborWidth-1)/2) + 256 * (y+(GaborHeight-1)/2)].Re = real;
			KernelFFT2[(x+(GaborWidth-1)/2) + 256 * (y+(GaborHeight-1)/2)].Im = img;
		}
}
double KernelRealPart(int x, int y, int Orientation, int Frequency)
{
	double U, V;
	double Sigma, Kv, Qu;
	double tmp1, tmp2;
	U = Orientation;
	V = Frequency;
	Sigma = 2 * (Math::PI) * (Math::PI);
	Kv = Math::PI * Math::Exp((-(V+2)/2)*Math::Log(2, Math::E));
	Qu = U * Math::PI  / 8;
	tmp1 = Math::Exp(-(Kv * Kv * ( x*x + y*y)/(2 * Sigma)));
	tmp2 = Math::Cos(Kv * Math::Cos(Qu) * x + Kv * Math::Sin(Qu) * y) - Math::Exp(-(Sigma/2));
	return tmp1 * tmp2 * Kv * Kv / Sigma;   
}

double KernelImgPart(int x, int y, int Orientation, int Frequency)
{
	double U, V;
	double Sigma, Kv, Qu;
	double tmp1, tmp2;
	U = Orientation;
	V = Frequency;
	Sigma = 2 * Math::PI * Math::PI;
	Kv = Math::PI * Math::Exp((-(V+2)/2)*Math::Log(2, Math::E));
	Qu = U * Math::PI  / 8;
	tmp1 = Math::Exp(-(Kv * Kv * ( x*x + y*y)/(2 * Sigma)));
	tmp2 = Math::Sin(Kv * Math::Cos(Qu) * x + Kv * Math::Sin(Qu) * y) - Math::Exp(-(Sigma/2));
	return tmp1 * tmp2 * Kv * Kv / Sigma;   
}
//*/

Matrix DoTransform(Matrix image, Complex KernelFFT2[],int Orientation, int Frequency, GaborResult result)//int image[]
{
	//准备图片FFT变换
	Complex imageFFT[] = new Complex[256 * 256];
	Complex ConvolutionResult[] = new Complex[256 * 256];
	//Complex KernelFFT2[]=new Complex[256*256];

	for(int y=(GaborWidth-1)/2; y<(GaborWidth-1)/2+128; y++)
		for(int x=(GaborWidth-1)/2; x<(GaborWidth-1)/2+128; x++)
			imageFFT[y * 256 + x].Re = image(y-(GaborWidth-1)/2,x-(GaborWidth-1)/2);

	Fourier::FFT2(imageFFT, 256, 256, FourierDirection::Forward );

	//在频域执行卷积
	for(int i=0; i<256*256; i++)
		ConvolutionResult[i] = imageFFT[i] * KernelFFT2[i] / (256 * 256); 

	Matrix dspFFTConvReal(256,256),dspFFTConvIm(256,256);
	for(int i=0;i<256;i++)
		for(int j=0;j<256;j++)
		{
			dspFFTConvReal(i,j)=ConvolutionResult[i*256+j].Re;
			dspFFTConvIm(i,j)=ConvolutionResult[i*256+j].Im;
		}
	SaveMatrix(dspFFTConvReal,"dspFFTConvReal");
	SaveMatrix(dspFFTConvIm,"dspFFTConvIm");

	Fourier::FFT2(ConvolutionResult, 256, 256, FourierDirection::Backward);

	Matrix dspConvReal(256,256),dspConvIm(256,256);
	for(int i=0;i<256;i++)
		for(int j=0;j<256;j++)
		{
			dspConvReal(i,j)=ConvolutionResult[i*256+j].Re;
			dspConvIm(i,j)=ConvolutionResult[i*256+j].Im;
		}
	SaveMatrix(dspConvReal,"dspConvReal");
	SaveMatrix(dspConvIm,"dspConvIm");

	//计算均值、方差及结果
	double Sum=0, Avg=0, Deta=0;
	double tmpModulus=0;
	//double tmpMag[] = new double(256, 256);
	Matrix tmpMag(256,256);

	for(int y=(GaborWidth-1); y<(GaborWidth-1)+128; y++)
		for(int x=(GaborWidth-1); x<(GaborWidth-1)+128; x++)
		{
			tmpModulus = ConvolutionResult[y*256+x].GetModulus();

//			Sum += tmpModulus;
			tmpMag(y,x) = tmpModulus;
		}

//		Avg = Sum / (128*128);
/*
		//计算方差
		for(int y=(GaborWidth-1); y<(GaborWidth-1)+128; y++)
			for(int x=(GaborWidth-1); x<(GaborWidth-1)+128; x++)
				Deta += (tmpMag(y,x) - Avg)*(tmpMag(y,x) - Avg);

		Deta = Deta / (128 * 128);
//*/
		//输出
		//*	
		Matrix Output(128,128);
		for(int y=(GaborWidth-1)+0; y < (GaborWidth-1)+128; y++)
		{
			for(int x=(GaborWidth-1)+0; x < (GaborWidth-1)+128; x++)
				Output(y-GaborWidth+1,x-GaborWidth+1)=tmpMag(y,x);
		//Console.Write("{0,8:F2}",tmpMag[y,x]);
		//Console.WriteLine();
		}
		//*/
		//放置结果
//		result.Avg = Avg;
//		result.Deta = Deta;
		return Output;//tmpMag
}

BOOL InvertGaussJordan(Matrix &input, Matrix &output)
{
	int col=input.ColNo();
	int *pnRow,*pnCol,i,j,k,l,u,v;
	double d=0,p=0;
	double *m_pData=new double[col*col];//m_pData[col*col]

	for(i=0;i<col;i++)
		for(j=0;j<col;j++)
			m_pData[i*col+j]=input(i,j);
			//m_pData[i+j*col]=input(i,j);
	//分配内存
	pnRow=new int[col];
	pnCol=new int[col];
	if(pnRow==NULL || pnCol==NULL)
		return FALSE;

	//消元
	for(k=0;k<col;k++)
	{
		d=0.0;
		for(i=k;i<col;i++)
		{
			for(j=k;j<col;j++)
			{
				l=i*col+j;
				p=fabs(m_pData[l]);
				if(p>d)
				{
					d=p;
					pnRow[k]=i;
					pnCol[k]=j;
				}
			}
		}//end i

		//failed
		if(d==0.0)
		{
			delete [] pnRow;
			delete [] pnCol;
			return FALSE;
		}

		if(pnRow[k]!=k)
		{
			for(j=0;j<col;j++)
			{
				u=k*col+j;
				v=pnRow[k]*col+j;
				p=m_pData[u];
				m_pData[u]=m_pData[v];
				m_pData[v]=p;
			}
		}//end if pnRow[k]

		if(pnCol[k]!=k)
		{
			for(i=0;i<col;i++)
			{
				u=i*col+k;
				v=i*col+pnCol[k];
				p=m_pData[u];
				m_pData[u]=m_pData[v];
				m_pData[v]=p;
			}
		}//end if pnCol[k]

		l=k*col+k;
		m_pData[l]=1.0/m_pData[l];
		for(j=0;j<col;j++)
		{
			if(j!=k)
			{
				u=k*col+j;
				m_pData[u]=m_pData[u]*m_pData[l];
			}
		}//end for j

		for(i=0;i<col;i++)
		{
			if(i!=k)
			{
				for(j=0;j<col;j++)
				{
					if(j!=k)
					{
						u=i*col+j;
						m_pData[u]=m_pData[u]-m_pData[i*col+k]*m_pData[k*col+j];
					}
				}
			}
		}//end for i

		for(i=0;i<col;i++)
		{
			if(i!=k)
			{
				u=i*col+k;
				m_pData[u]=-m_pData[u]*m_pData[l];
			}
		}//end for i
	}

	//调整恢复行列次序
	for(k=col-1;k>=0;k--)
	{
		if(pnCol[k]!=k)
		{
			for(j=0;j<col;j++)
			{
				u=k*col+j;
				v=pnCol[k]*col+j;
				p=m_pData[u];
				m_pData[u]=m_pData[v];
				m_pData[v]=p;
			}
		}//end if pnCol[k]

		if(pnRow[k]!=k)
		{
			for(i=0;i<col;i++)
			{
				u=i*col+k;
				v=i*col+pnRow[k];
				p=m_pData[u];
				m_pData[u]=m_pData[v];
				m_pData[v]=p;
			}
		}//end if pnRow[k]
	}
	
	for(i=0;i<col;i++)
		for(j=0;j<col;j++)
			output(i,j)=m_pData[i*col+j];
			//output(i,j)=m_pData[i+j*col];	

	//清理内存
	delete [] pnRow;
	delete [] pnCol;
	delete [] m_pData;
	return TRUE;
}

Matrix Minus(Matrix &m1,Matrix &m2)
{
	int row=m1.RowNo();
	int col=m1.ColNo();
	Matrix m(row,col);
	for(int i=0;i<row;i++)
	{
		for(int j=0;j<col;j++)
		{
			m(i,j)=m1(i,j)-m2(i,j);
			if(m(i,j)<0.000001)//100*m(i,0)==m(i,0))
				m(i,j)=0.0;
		}
	}	

	return m;
}

#pragma unmanaged