www.pudn.com > FFT2Dsource.rar > PinYuLuBoDib.cpp

#include "stdafx.h" 
#include "windowsx.h" 
#include "math.h" 
#include "PinYuLuBoDib.h" 
#include "malloc.h" 
#include "MainFrm.h" 
#include "DSplitView.h" 
#include "Cdib.h" 
#define WIDTHBYTES(bits)    (((bits) + 31) / 32 * 4)  
 
 
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr 
#define pi 3.14159265359 
  
  
PinYuLuBoDib::PinYuLuBoDib() 
{ 
} 
PinYuLuBoDib::~PinYuLuBoDib() 
{ 
} 
 
CNumber PinYuLuBoDib::Add(CNumber c1,CNumber c2) 
{ 
	CNumber c; 
	c.re=c1.re+c2.re; 
	c.im=c1.im+c2.im; 
	return c; 
} 
CNumber PinYuLuBoDib::Sub(CNumber c1,CNumber c2) 
{ 
	CNumber c; 
	c.re=c1.re-c2.re; 
	c.im=c1.im-c2.im; 
	return c; 
} 
CNumber PinYuLuBoDib::Mul(CNumber c1,CNumber c2) 
{ 
	CNumber c; 
	c.re=c1.re*c2.re-c1.im*c2.im; 
	c.im=c1.re*c2.im+c2.re*c1.im; 
	return c; 
} 
 
/////////////////////////////////////////////////// 
//此函数实现快速傅立叶变换 
//参数t、f分别是指向时域和频域的指针，power是2的幂数 
/////////////////////////////////////////////////// 
void PinYuLuBoDib::QFC(CNumber* t,CNumber* f,int power) 
{ 
	long count;//傅立叶变换点数 
	int i,j,k,p,bfsize; 
	CNumber *w,*x,*a,*b;//复数结构类型的指针变量，其中w指向加权系数 
	double angle;//计算加权系数所用的角度值 
	count=1<byBitCount==8)	//灰度图像 
		   p_data=this->GetData();//指向原图象数据区 
	else	//24位真彩色 
		   p_data=this->GetData2(); 
	width=this->GetWidth();//得到图象宽度 
	height=this->GetHeight();//得到图象高度 
    long lLineBytes=WIDTHBYTES(width*8);//计算图象每行的字节数 
	while(w*2<=width)//计算进行傅立叶变换的宽度（2的整数次方） 
	{ 
		w*=2; 
		wp++; 
	} 
	while(h*2<=height)//计算进行傅立叶变换的高度（2的整数次方） 
	{ 
		h*=2; 
		hp++; 
	} 
	t=new CNumber[w*h];//分配存储器空间 
	f=new CNumber[w*h]; 
	for(j=0;j255) 
				temp=255; 
			 
			p=p_data+lLineBytes*(height-(jbyBitCount==8)//灰度图像 
		   p_data=this->GetData();//指向原图象数据区 
	else//24位真彩色 
		   p_data=this->GetData2(); 
	width=this->GetWidth();//得到图象宽度 
	height=this->GetHeight();//得到图象高度 
    long lLineBytes=WIDTHBYTES(width*8);//计算图象每行的字节数 
	while(w*2<=width)//计算进行傅立叶变换的宽度（2的整数次方） 
	{ 
		w*=2; 
		wp++; 
	} 
	while(h*2<=height)//计算进行傅立叶变换的高度（2的整数次方） 
	{ 
		h*=2; 
		hp++; 
	} 
	t=new CNumber[w*h];//分配存储器空间 
	f=new CNumber[w*h]; 
	for(j=0;j255) 
				temp=255; 
			p=p_data+lLineBytes*(height-(j=1;idim--)  
	{ 
		n=nn[idim]; 
		nrem=ntot/(n*nprev); 
		ip1=nprev << 1; 
		ip2=ip1*n; 
		ip3=ip2*nrem; 
		i2rev=1; 
		for (i2=1;i2<=ip2;i2+=ip1) 
		{ 
			if (i2 < i2rev)  
			{ 
				for (i1=i2;i1<=i2+ip1-2;i1+=2)  
				{ 
					for (i3=i1;i3<=ip3;i3+=ip2) 
					{ 
						i3rev=i2rev+i3-i2; 
						SWAP(data[i3],data[i3rev]); 
						SWAP(data[i3+1],data[i3rev+1]); 
					} 
				} 
			} 
			ibit=ip2 >> 1; 
			while (ibit >= ip1 && i2rev > ibit) 
			{ 
				i2rev -= ibit; 
				ibit >>= 1; 
			} 
			i2rev += ibit; 
		} 
		ifp1=ip1; 
		while (ifp1 < ip2)  
		{ 
			ifp2=ifp1 << 1; 
			theta=isign*pi*2/(ifp2/ip1); 
			wtemp=sin(0.5*theta); 
			wpr = -2.0*wtemp*wtemp; 
			wpi=sin(theta); 
			wr=1.0; 
			wi=0.0; 
			for (i3=1;i3<=ifp1;i3+=ip1) 
			{ 
				for (i1=i3;i1<=i3+ip1-2;i1+=2)  
				{ 
					for (i2=i1;i2<=ip3;i2+=ifp2)  
					{ 
						k1=i2; 
						k2=k1+ifp1; 
						tempr=wr*data[k2]-wi*data[k2+1]; 
						tempi=wr*data[k2+1]+wi*data[k2]; 
						data[k2]=data[k1]-tempr; 
						data[k2+1]=data[k1+1]-tempi; 
						data[k1] += tempr; 
						data[k1+1] += tempi; 
					} 
				} 
				wr=(wtemp=wr)*wpr-wi*wpi+wr; 
				wi=wi*wpr+wtemp*wpi+wi; 
			} 
			ifp1=ifp2; 
		} 
		nprev *= n; 
	} 
} 
 
/************************************************************************* 
* 函数：BWFilterL(int u,int v,int u1,int v1) 
*参数：u、v分别是截止频率的x、y分量值，由用户给定 
*功能：此函数用来实现图象的布特沃斯低通滤波 
/*************************************************************************/ 
 
void PinYuLuBoDib::BWFilterL(int u,int v,int n) 
{ 
	LPBYTE  p_data, p;//指向原图象数据区指针 
	int width,height;//原图象的宽度和高度        
	int i,j; 
	double max=0.0,d0,d;//中间变量 
	double *t,*H; 
	if(this->byBitCount==8)//灰度图像 
		   p_data=this->GetData();//指向原图象数据区 
	else//24位真彩色 
		   p_data=this->GetData2();//指向原图象数据区 
	width=this->GetWidth();//得到图象宽度 
	height=this->GetHeight();//得到图象高度 
    long lLineBytes=WIDTHBYTES(width*8);//计算图象每行的字节数 
	t=new double [height*lLineBytes*2+1];//分配存储器空间 
	H=new double [height*lLineBytes*2+1]; 
	d0=sqrt(u*u+v*v);//计算截止频率d0 
	for(j=0;jbyBitCount==8)//灰度图像 
		   p_data=this->GetData();//指向原图象数据区 
	else//24位真彩色 
		   p_data=this->GetData2();//指向原图象数据区 
	width=this->GetWidth();//得到图象宽度 
	height=this->GetHeight();//得到图象高度 
    long lLineBytes=WIDTHBYTES(width*8);//计算图象每行的字节数 
	t=new double [height*lLineBytes*2+1];//分配存储器空间 
	H=new double [height*lLineBytes*2+1]; 
	d0=sqrt(u*u+v*v);//计算截止频率d0 
	for(j=0;jbyBitCount==8)//灰度图像 
		p_data=this->GetData();//指向原图象数据区 
	else//24位真彩色 
		p_data=this->GetData2();//指向原图象数据区 
	width=this->GetWidth();//得到图象宽度 
	height=this->GetHeight();//得到图象高度 
    long lLineBytes=WIDTHBYTES(width*8);//计算图象每行的字节数 
	t=new double [height*lLineBytes*2+1];//分配存储器空间 
	H=new double [height*lLineBytes*2+1]; 
	d0=sqrt(u*u+v*v);//计算截止频率d0 
	for(j=0;jbyBitCount==8)//灰度图像 
		p_data=this->GetData();//指向原图象数据区 
	else//24位真彩色 
		p_data=this->GetData2();//指向原图象数据区 
	width=this->GetWidth();//得到图象宽度 
	height=this->GetHeight();//得到图象高度 
    long lLineBytes=WIDTHBYTES(width*8);//计算图象每行的字节数 
	t=new double [height*lLineBytes*2+1];//分配存储器空间 
	H=new double [height*lLineBytes*2+1]; 
	d0=sqrt(u*u+v*v);//计算截止频率d0 
	for(j=0;jbyBitCount==8)//灰度图像 
		   p_data=this->GetData();//指向原图象数据区 
	else//24位真彩色 
		   p_data=this->GetData2();//指向原图象数据区 
	width=this->GetWidth();//得到图象宽度 
	height=this->GetHeight();//得到图象高度 
    long lLineBytes=WIDTHBYTES(width*8);//计算图象每行的字节数 
	t=new double [height*lLineBytes*2+1];//分配存储器空间 
	H=new double [height*lLineBytes*2+1]; 
	d0=sqrt(u*u+v*v);//计算截止频率d0 
	d1=sqrt(u1*u1+v1*v1); 
	for(j=0;jd1) 
			{ 
				H[2*i+(2*lLineBytes)*j+1]=0.0; 
				H[2*i+(2*lLineBytes)*j+2]=0.0; 
			} 
			else 
			{ 
				H[2*i+(2*lLineBytes)*j+1]=(d-d1)/(d0-d1); 
				H[2*i+(2*lLineBytes)*j+2]=0.0; 
			}				 
		} 
	} 
    fourier(t,height,lLineBytes,1); 
	for(j=1;jbyBitCount==8)//灰度图像 
		p_data=this->GetData();//指向原图象数据区 
	else//24位真彩色 
		p_data=this->GetData2();//指向原图象数据区 
	width=this->GetWidth();//得到图象宽度 
	height=this->GetHeight();//得到图象高度 
    long lLineBytes=WIDTHBYTES(width*8);//计算图象每行的字节数 
	t=new double [height*lLineBytes*2+1];//分配存储器空间 
	H=new double [height*lLineBytes*2+1]; 
	d0=sqrt(u*u+v*v);//计算截止频率d0 
	d1=sqrt(u1*u1+v1*v1); 
	for(j=0;jd1) 
			{ 
				H[2*i+(2*lLineBytes)*j+1]=1; 
				H[2*i+(2*lLineBytes)*j+2]=0.0; 
			} 
			else 
			{ 
				H[2*i+(2*lLineBytes)*j+1]=(d-d0)/(d1-d0); 
				H[2*i+(2*lLineBytes)*j+2]=0.0; 
			}						 
		} 
	} 
    fourier(t,height,lLineBytes,1); 
	for(j=1;jbyBitCount==8)//灰度图像 
		p_data=this->GetData();//指向原图象数据区 
	else//24位真彩色 
		p_data=this->GetData2();//指向原图象数据区 
	width=this->GetWidth();//得到图象宽度 
	height=this->GetHeight();//得到图象高度 
    long lLineBytes=WIDTHBYTES(width*8);//计算图象每行的字节数 
	t=new double [height*lLineBytes*2+1];//分配存储器空间 
	H=new double [height*lLineBytes*2+1]; 
	d0=sqrt(u*u+v*v);//计算截止频率d0 
	for(j=0;jbyBitCount==8)//灰度图像 
		p_data=this->GetData();//指向原图象数据区 
	else//24位真彩色 
		p_data=this->GetData2();//指向原图象数据区 
	width=this->GetWidth();//得到图象宽度 
	height=this->GetHeight();//得到图象高度 
    long lLineBytes=WIDTHBYTES(width*8);//计算图象每行的字节数 
	t=new double [height*lLineBytes*2+1];//分配存储器空间 
	H=new double [height*lLineBytes*2+1]; 
	d0=sqrt(u*u+v*v);//计算截止频率d0 
	for(j=0;j