fftgabor.cpp

关于字符分割的程序

// FFTGabor.cpp: implementation of the FFTGabor class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "hwpre.h"
#include "FFTGabor.h"
#include "normalTran.h"
#include "savetofile.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif
#define PI 3.1415926535
//#define MathE 2.718282
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
//卷积核大小
/*
    128*128 图像时，101*101,图像扩展到256*256
	32*32 图像时， 31*31,图像扩展到 64*64
*/
FFTGabor::FFTGabor(int  GaborW,int compW)
{
	GaborWidth=GaborW;
	GaborHeight=GaborW;
	complexWidth= compW;
	complexHeight= compW;
	//KernelFFT2 = new complex<double>[32];
	//KernelFFT2 = new complex<double>[32*256*256];
	int i;
	for (i=0;i<32;i++) {
		KernelFFT2[i] = new complex<double>[compW*compW];
	}
}

FFTGabor::~FFTGabor()
{
	int i;
	for (i=0;i<32;i++) {
		delete[] KernelFFT2[i];
	}

}

void FFTGabor::PrepareKernel()
{
	int Orientation,Frequency;

	for( Orientation=0; Orientation<8; Orientation++)
		for( Frequency=0; Frequency<4; Frequency++)
		{
			savetofile fs;
			CalculateKernel(Orientation, Frequency);
			FFT2(KernelFFT2[Orientation * 4 + Frequency], complexWidth, complexHeight );
		}
}

void FFTGabor::CalculateKernel(double Orientation, double Frequency)
{
		int position = Orientation * 4 + Frequency;

		double Sigma2, Kv, Qu;
		double tmp1, tmp2, tmp3;
		double realT, imagT;
		
		Sigma2 = 2 * PI * PI;
		for(int x = -(GaborWidth-1)/2; x<(GaborWidth-1)/2+1; x++)
			for(int y = -(GaborHeight-1)/2; y<(GaborHeight-1)/2+1; y++)
			{
				Kv = PI * pow(2, (-(Frequency+2)/2));
				Qu = Orientation * PI / 8;
				tmp1 = exp(-(Kv * Kv * ( x*x + y*y)/(2 * Sigma2)));
				
				tmp2 = sin(Kv * cos(Qu) * x + Kv * sin(Qu) * y) - exp(-(Sigma2/2));
				tmp3 = cos(Kv * cos(Qu) * x + Kv * sin(Qu) * y) - exp(-(Sigma2/2));

				imagT =tmp1 * tmp2 * Kv * Kv / Sigma2;
				realT =tmp1 * tmp3 * Kv * Kv / Sigma2;


				KernelFFT2[position][(x+(GaborWidth-1)/2) + complexWidth * (y+(GaborHeight-1)/2)] = complex<double>(realT,imagT);
			}
}

void FFTGabor::GaborTransform(int * lpDIBBits, LONG lWidth, LONG lHeight, double * TransformResult)
{
	// 循环变量
	LONG	i;
	LONG	j;
	int Orientation, Frequency ;
	// 分配内存
	complex<double> *TD = new complex<double>[complexWidth * complexHeight];
	complex<double> *FD = new complex<double>[complexWidth * complexHeight];
	
	// 行
	for(i = 0; i < lHeight; i++) 
	{
		// 列
		for(j = 0; j < lWidth; j++)
		{
			// 给时域赋值
			TD[j + complexWidth * i] = complex<double>(lpDIBBits[j+lWidth * i], 0);
		}
	}
	::FFT2(TD,complexWidth , complexHeight);
	
	//::IFFT2(TD,256,256);
	//在频域执行卷积
	double scale = 1.0 / (complexWidth * complexHeight); 
	double Sum=0, Avg=0, Deta=0,Sum2 = 0;
	double tmpModulus=0;
	//double * tmpMag = new double[lWidth*lHeight] ;// [128][128] ;
	double *  tmpMag = new double[complexHeight* complexWidth] ;
	int x ,y;
	
	for( Orientation=0; Orientation<8; Orientation++){
		for( Frequency=0; Frequency<4; Frequency++)
		{
			for( i=0; i<complexWidth * complexHeight; i++)
				FD[i] = TD[i] * KernelFFT2[Orientation * 4 + Frequency][i] ; // * scale; 
			
			IFFT2(FD,complexWidth , complexHeight);
			//
			//savetofile fs;
			//fs.filesave(FD,64,64);

			//计算均值、方差及结果
			Sum=0, Avg=0, Deta=0;Sum2 = 0;
			
			//for( y=(GaborHeight/2); y<(GaborHeight/2)+lHeight; y++){
			//	for( x=(GaborWidth/2); x<(GaborWidth/2)+lWidth; x++)
			for( y=0; y<complexHeight; y++){
				for( x=0; x<complexWidth; x++)
				{
					tmpModulus = sqrt(FD[y*complexWidth+x].real() * FD[y*complexWidth+x].real() +
						FD[y*complexWidth+x].imag() * FD[y*complexWidth+x].imag());
					
					Sum += tmpModulus;
					Sum2+= tmpModulus * tmpModulus;
					tmpMag[y*complexWidth+x] = tmpModulus;
				}
			}
			Avg = Sum / (complexHeight* complexWidth);
			//计算方差
//			for( y=0; y < complexHeight; y++){
//				for( x=0; x < complexWidth; x++)
//					Deta += (tmpMag[y* complexWidth + x] - Avg)*(tmpMag[y* complexWidth + x] - Avg);
//			}
			Deta = complexWidth*complexHeight * Sum2 - Sum * Sum ;
			Deta =Deta / (complexWidth * complexHeight)/(complexWidth * complexHeight) ;
			

			TransformResult[(Orientation * 4 + Frequency)*2] = Avg;
			TransformResult[(Orientation * 4 + Frequency)*2+1] = Deta;
		}
	}
	
	delete[] tmpMag;
	delete[] FD;
	delete[] TD;
}

void FFTGabor::GaborTransform(int * lpDIBBits, LONG lWidth, LONG lHeight, int Orientation, int Frequency, FFTGaborResult * result)
{
	LONG	i;
	LONG	j;

	complex<double> *TD = new complex<double>[complexWidth * complexHeight];
	complex<double> *FD = new complex<double>[complexWidth * complexHeight];
	
	// 行
	for(i = 0; i < lHeight; i++) 
	{
		// 列
		for(j = 0; j < lWidth; j++)
		{
			// 给时域赋值
			TD[j + complexWidth * i] = complex<double>(lpDIBBits[j+lWidth * i], 0);
		}
	}
	::FFT2(TD,complexWidth , complexHeight);
	
	//::IFFT2(TD,256,256);
	//在频域执行卷积
	double scale = 1.0 / (complexWidth * complexHeight); 
	
	for( i=0; i<complexWidth * complexHeight; i++)
		FD[i] = TD[i] * KernelFFT2[Orientation * 4 + Frequency][i] ; // * scale; 
	
	IFFT2(FD,complexWidth , complexHeight);
	
	//计算均值、方差及结果
	//计算均值并找出tmpMag中的最大最小值，以便调整到0～255，用于显示输出。
	double min, max;
	double Sum=0, Avg=0, Deta=0;
	double tmpModulus=0;
	double * tmpMag = new double[lWidth*lHeight] ;// [128][128] ;
	int x ,y;
	tmpModulus = sqrt(FD[(GaborHeight/2)*complexWidth+(GaborWidth/2)].real() *
								FD[(GaborHeight/2)*complexWidth+(GaborWidth/2)].real() +
								FD[(GaborHeight/2)*complexWidth+(GaborWidth/2)].imag() *
								FD[(GaborHeight/2)*complexWidth+(GaborWidth/2)].imag());
	min = max = tmpModulus;
	for( y=(GaborHeight/2); y<(GaborHeight/2)+lHeight; y++){
		for( x=(GaborWidth/2); x<(GaborWidth/2)+lWidth; x++)
		{
			tmpModulus = sqrt(FD[y*complexWidth+x].real() * FD[y*complexWidth+x].real() +
				FD[y*complexWidth+x].imag() * FD[y*complexWidth+x].imag());
			
			if(min>tmpModulus)
				min = tmpModulus;
			if(max<tmpModulus)
				max = tmpModulus;	
			
			Sum += tmpModulus;
			tmpMag[(y-(GaborHeight/2))* lWidth + (x-(GaborWidth/2))] = tmpModulus;
		}
	}
	Avg = Sum / (lHeight * lWidth);
	//计算方差
	for( y=0; y < lHeight; y++){
		for( x=0; x < lWidth; x++)
			Deta += (tmpMag[y* lWidth + x] - Avg)*(tmpMag[y* lWidth + x] - Avg);
	}
	Deta = Deta / (lHeight * lWidth);
	
	//计算magShow以及均值
	scale = 255.0/(max-min);
	
	for( y=0; y < lHeight; y++){
		for( x=0; x < lWidth; x++)
		{
			result->magShow[y * lWidth + x] = (int)(scale*(tmpMag[y* lWidth + x]-min));
		}
	}
	result->Avg = Avg;
	result->Deta = Deta;
	
	delete[] tmpMag;

	delete[] FD;
	delete[] TD;
}

void FFTGabor::intTocomplex(complex<double> * TD,int * lpDIBBits,LONG lWidth, LONG lHeight)
{
	LONG	i;
	LONG	j;

//	complex<double> *TD = new complex<double>[complexWidth * complexHeight];
	
	// 行
	for(i = 0; i < lHeight; i++) 
	{
		// 列
		for(j = 0; j < lWidth; j++)
		{
			// 给时域赋值
			TD[j + complexWidth * i] = complex<double>(lpDIBBits[j+lWidth * i], 0);
		}
	}
	::FFT2(TD,complexWidth , complexHeight);

}
void FFTGabor::GaborTransform(complex<double> *TD, LONG lWidth, LONG lHeight, double * TransformResult)
{
	// 循环变量
	LONG	i;
//	LONG	j;
	int Orientation, Frequency ;
	// 分配内存
	complex<double> *FD = new complex<double>[complexWidth * complexHeight];
	
	double scale = 1.0 / (complexWidth * complexHeight); 
	double Sum=0, Avg=0, Deta=0;
	double tmpModulus=0;
	double * tmpMag = new double[lWidth*lHeight] ;// [128][128] ;
	int x ,y;
	
	for( Orientation=0; Orientation<8; Orientation++){
		for( Frequency=0; Frequency<4; Frequency++)
		{
			//在频域执行卷积
			for( i=0; i<complexWidth * complexHeight; i++)
				FD[i] = TD[i] * KernelFFT2[Orientation * 4 + Frequency][i] ; // * scale; 
			
			IFFT2(FD,complexWidth , complexHeight);
			
			//计算均值、方差及结果
			Sum=0, Avg=0, Deta=0;
			
			for( y=(GaborHeight/2); y<(GaborHeight/2)+lHeight; y++){
				for( x=(GaborWidth/2); x<(GaborWidth/2)+lWidth; x++)
				{
					tmpModulus = sqrt(FD[y*complexWidth+x].real() * FD[y*complexWidth+x].real() +
						FD[y*complexWidth+x].imag() * FD[y*complexWidth+x].imag());
					
					Sum += tmpModulus;
					tmpMag[(y-(GaborHeight/2)) * lWidth + (x-(GaborWidth/2))] = tmpModulus;
				}
			}
			Avg = Sum / (lHeight * lWidth);
			//计算方差
			for( y=0; y < lHeight; y++){
				for( x=0; x < lWidth; x++)
					Deta += (tmpMag[y* lWidth + x] - Avg)*(tmpMag[y* lWidth + x] - Avg);
			}
			Deta = Deta / (lHeight * lWidth);
			

			TransformResult[(Orientation * 4 + Frequency)*2] = Avg;
			TransformResult[(Orientation * 4 + Frequency)*2+1] = Deta;
		}
	}
	
	delete[] tmpMag;
	delete[] FD;
}

void FFTGabor::GaborTransform(complex<double> *TD, LONG lWidth, LONG lHeight, int Orientation, int Frequency, FFTGaborResult * result)
{
	LONG	i;
//	LONG	j;

	complex<double> *FD = new complex<double>[complexWidth * complexHeight];
	

	//在频域执行卷积
	double scale = 1.0 / (complexWidth * complexHeight); 
	
	for( i=0; i<complexWidth * complexHeight; i++)
		FD[i] = TD[i] * KernelFFT2[Orientation * 4 + Frequency][i] ; // * scale; 
	
	IFFT2(FD,complexWidth , complexHeight);
	
	//计算均值、方差及结果
	//计算均值并找出tmpMag中的最大最小值，以便调整到0～255，用于显示输出。
	double min, max;
	double Sum=0, Avg=0, Deta=0;
	double tmpModulus=0;
	double * tmpMag = new double[lWidth*lHeight] ;// [128][128] ;
	int x ,y;
	tmpModulus = sqrt(FD[(GaborHeight/2)*complexWidth+(GaborWidth/2)].real() *
								FD[(GaborHeight/2)*complexWidth+(GaborWidth/2)].real() +
								FD[(GaborHeight/2)*complexWidth+(GaborWidth/2)].imag() *
								FD[(GaborHeight/2)*complexWidth+(GaborWidth/2)].imag());
	min = max = tmpModulus;
	for( y=(GaborHeight/2); y<(GaborHeight/2)+lHeight; y++){
		for( x=(GaborWidth/2); x<(GaborWidth/2)+lWidth; x++)
		{
			tmpModulus = sqrt(FD[y*complexWidth+x].real() * FD[y*complexWidth+x].real() +
				FD[y*complexWidth+x].imag() * FD[y*complexWidth+x].imag());
			
			if(min>tmpModulus)
				min = tmpModulus;
			if(max<tmpModulus)
				max = tmpModulus;	
			
			Sum += tmpModulus;
			tmpMag[(y-(GaborHeight/2))* lWidth + (x-(GaborWidth/2))] = tmpModulus;
		}
	}
	Avg = Sum / (lHeight * lWidth);
	//计算方差
	for( y=0; y < lHeight; y++){
		for( x=0; x < lWidth; x++)
			Deta += (tmpMag[y* lWidth + x] - Avg)*(tmpMag[y* lWidth + x] - Avg);
	}
	Deta = Deta / (lHeight * lWidth);
	
	//计算magShow以及均值
	scale = 255.0/(max-min);
	
	for( y=0; y < lHeight; y++){
		for( x=0; x < lWidth; x++)
		{
			result->magShow[y * lWidth + x] = (int)(scale*(tmpMag[y* lWidth + x]-min));
		}
	}
	result->Avg = Avg;
	result->Deta = Deta;
	
	delete[] tmpMag;

	delete[] FD;
}

/************************************************************************/
/* 
	private void BeginProcess()
	{
		Image imgSrc = this.spbSrc.PicBox.Image;
		Bitmap imgDest = new Bitmap(1024, 512, PixelFormat.Format24bppRgb);
		Bitmap bmp = new Bitmap(imgSrc);

		int[,] image = new int[128,128];
		int part = 0, totalPart = (imgSrc.Width*imgSrc.Height)/(128*128);
		TransformResult = new double[totalPart, 64];

		lblMessage.Text = "正在准备 Gabor 核函数...";
		FFTGabor gabor = new FFTGabor();
		gabor.PrepareKernel();

		FFTGaborResult result = new FFTGaborResult();

		for(int i=0; i<imgSrc.Height; i+=128)
			for(int j=0; j<imgSrc.Width; j+=128)
			{
				part++;
				lblMessage.Text = "开始处理第 " + part.ToString() + " 部分，共 " + totalPart.ToString() + " 部分，请等待...";
 
				for(int y=0; y<128; y++)
					for(int x=0; x<128; x++)
						image[y,x]=ColorToGray(bmp.GetPixel(j+x,i+y));

				for(int Orientation=0; Orientation<8; Orientation++)
					for(int Frequency=0; Frequency<4; Frequency++)
					{
						gabor.GaborTransform(image, Orientation, Frequency, result);

						RunProgress.Value = (Orientation * 4 + Frequency + 1) * 100 / 32;

						for(int y=0; y<128; y++)
							for(int x=0; x<128; x++)
								imgDest.SetPixel(Orientation * 128 + x, Frequency * 128 + y, 
									Color.FromArgb(result.magShow[x,y], result.magShow[x,y], result.magShow[x,y]));

						TransformResult[part-1, (Orientation * 4 + Frequency)*2] = result.Avg;
						TransformResult[part-1, (Orientation * 4 + Frequency)*2 + 1] = result.Deta;

						if(stop)
						{
							btnStop.Enabled = false;
							btnRun.Enabled = true;

							lblMessage.Text = "用户终止了程序的运行!";
							return;
						}
					}

				this.spbDest.PicBox.Image = imgDest;
			}

		btnStop.Enabled = false;
		btnRun.Enabled = true;
		Finished = true;
		lblMessage.Text = "处理完成。";
	}

	private int ColorToGray(Color color)
	{
		//Gray = 0.299 * R + 0.587 * G + 0.114 * B
		//当然把浮点运算去掉，都用整数运算会更快：Gary = (R * 77 + G * 151 + B * 28) >> 8; 
		int gray = ((((color.R * 77) + (color.G * 151)) + (color.B * 28)) >> 8);

		if(gray<0) gray=0;
		if(gray>255) gray=255;

		return gray;
	}

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