ffuriour.cpp

一些自己做的关于图象处理的程序（如傅立叶变换

// Ffuriour.cpp: implementation of the CFfuriour class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "photostar.h"
#include "Ffuriour.h"
#include "math.h"
COMPLEX * FastFTValue=NULL;
const double pai=3.141592653589793;

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

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

CFfuriour::CFfuriour()
{

}

CFfuriour::~CFfuriour()
{

}


COMPLEX CFfuriour::ReverceOrder(COMPLEX *pdData, int iAmount)//倒序
{
	
	int NV2=iAmount/2,NM1=iAmount-1;
	int i,j=0;
	int k;
	COMPLEX temple;
	for(i=0;i<NM1;i++)
 	{
		if(i<j)
 		{
			temple.x=pdData[j].x;
			temple.y=pdData[j].y;
	 		pdData[j].x=pdData[i].x;
			pdData[j].y=pdData[i].y;
			pdData[i].x=temple.x;
			pdData[i].y=temple.y;
	 	}
	 		k=NV2;
		while(k<=j)
		{
			j-=k;
			k=k/2;
		} 
	 	j+=k;
	 		
			
	}
	return *pdData;
		 	
}
COMPLEX CFfuriour::FastFuriourTranslate(COMPLEX *templeData, int iAmount)//傅立叶变换
{
	int m=1;
	int L=(int)(log10(iAmount)/log10(2));
	int LE,LE1,j,i,IP;
	COMPLEX U,W;
	COMPLEX temp;
	for(;m<=L;m++)
	{
	
		LE=(int)pow(2,m);
		LE1=LE/2;
		U.x=1,U.y=0;
		W.x=cos(pai/(double)LE1),W.y=-sin(pai/(double)LE1);
		
		for(j=0;j<=LE1-1;j++)
		{
			for(i=j;i<=iAmount-1;i+=LE)
			{
				IP=LE1+i;
				temp=Mul(templeData[IP],U);
				templeData[IP]=Sub(templeData[i],temp);
				templeData[i]=Add(templeData[i],temp);
				
			}
			U=Mul(U,W);    
		}
		
		
	}
	
	
	return *templeData;
}

COMPLEX CFfuriour::Add(COMPLEX a1, COMPLEX a2)//复数加
{
	COMPLEX a;
	a.x=a1.x+a2.x;
	a.y=a1.y+a2.y;
	return a;
}

COMPLEX CFfuriour::Mul(COMPLEX a1, COMPLEX a2)//复数乘
{
	COMPLEX a;
	a.x=a1.x*a2.x-a1.y*a2.y;
	a.y=a1.x*a2.y+a1.y*a2.x;
	return a;
}

COMPLEX CFfuriour::Sub(COMPLEX a1, COMPLEX a2)//复数减
{
	COMPLEX a;
	a.x=a1.x-a2.x;
	a.y=a1.y-a2.y;
	return a;
}

COMPLEX CFfuriour::Div(COMPLEX a1, COMPLEX a2)
{
	COMPLEX a;
	a.x = (a1.x*a2.x+a1.y*a2.y)/(pow(a2.x,2)+pow(a2.y,2));
    a.y = (a1.y*a2.x-a1.x*a2.y)/(pow(a2.x,2)+pow(a2.y,2));
    return a;
}

BOOL CFfuriour::OnFfuriour()
{   
	if (m_pImageObject == NULL)
	{
		return FALSE;
	}
	
	 iWidth = m_pImageObject->GetWidth();
	 iHeight = m_pImageObject->GetHeight();
	int iNewHeight=iHeight;
	int iNewWidth=iWidth;
	
	/******************************the old image*****************************/
	unsigned char *pOldBuffer, *pNewBuffer, *pOldBits, *pNewBits;
	BITMAPFILEHEADER *pOldBFH, *pNewBFH;
	BITMAPINFOHEADER *pOldBIH, *pNewBIH;
	RGBQUAD *pOldPalette, *pNewPalette;
	int nWidthBytes, nNumColors;
	
	pOldBuffer = (unsigned char *) m_pImageObject->GetDIBPointer(&nWidthBytes, m_pImageObject->GetNumBits());
	if (pOldBuffer == NULL) 
	{
		return FALSE;
	}
	
	pOldBFH = (BITMAPFILEHEADER *) pOldBuffer;
	pOldBIH = (BITMAPINFOHEADER *) &pOldBuffer[sizeof(BITMAPFILEHEADER)];
	nNumColors = m_pImageObject->GetNumColors();
	pOldPalette = (RGBQUAD *)  &pOldBuffer[sizeof(BITMAPFILEHEADER)+sizeof(BITMAPINFOHEADER)];
	pOldBits = (unsigned char *) &pOldBuffer[sizeof(BITMAPFILEHEADER)
		+sizeof(BITMAPINFOHEADER)+nNumColors*sizeof(RGBQUAD)];
	
	/*****************************the new image*****************************/
	//把长宽不是2的整数倍的图像补成2的整数倍
	double tempHeight,tempWidth;
	tempHeight = double (int (log10(iNewHeight)/log10(2)));
	tempWidth = double (int (log10(iNewWidth)/log10(2)));
	
	while (tempHeight != log10(iNewHeight)/log10(2)) 
	{
		iNewHeight++;
		tempHeight = double (int (log10(iNewHeight)/log10(2)));
	}
	while (tempWidth != log10(iNewWidth)/log10(2)) 
	{
		iNewWidth++;
		tempWidth = double (int (log10(iNewWidth)/log10(2)));
	}
	int nNewWidthBytes=iNewWidth;
	DWORD dwNewSize;
	HGLOBAL hNewDib;
	
	//Allocate a global memory block for the new image
	dwNewSize = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) 
		+ nNumColors * sizeof(RGBQUAD)
		+ nNewWidthBytes * (iNewHeight);
	hNewDib = ::GlobalAlloc(GMEM_MOVEABLE | GMEM_ZEROINIT, dwNewSize);
	if (hNewDib == NULL)
	{
		m_pImageObject->m_nLastError = IMAGELIB_MEMORY_ALLOCATION_ERROR;
		::GlobalUnlock(m_pImageObject->GetDib());
		return FALSE;
	}
	
	//Get the pointer to the new memory block
	pNewBuffer = (unsigned char *) ::GlobalLock(hNewDib);
	if (pNewBuffer == NULL)
	{
		::GlobalFree( hNewDib );
		m_pImageObject->m_nLastError = IMAGELIB_MEMORY_LOCK_ERROR;
		::GlobalUnlock(m_pImageObject->GetDib());
		return FALSE;
	}
	
	pNewBFH = (BITMAPFILEHEADER *) pNewBuffer;
	pNewBIH = (BITMAPINFOHEADER *) &pNewBuffer[sizeof(BITMAPFILEHEADER)];
	pNewPalette = (RGBQUAD *) &pNewBuffer[sizeof(BITMAPFILEHEADER)+sizeof(BITMAPINFOHEADER)];
	pNewBits = (unsigned char *) &pNewBuffer[sizeof(BITMAPFILEHEADER)+sizeof(BITMAPINFOHEADER)+nNumColors*sizeof(RGBQUAD)];
	
	//Make a copy of the old image
	memcpy(pNewBuffer,pOldBuffer,dwNewSize);
	//改变新图像的大小
	
	pNewBIH->biHeight = iNewHeight;
	while (pOldBIH->biWidth%4 !=0) //更改新图像的宽度信息时应先把旧图像的宽度信息补为4的倍数
	{							   //因为旧图像的nWidthBytes系统已经自动补为4的倍数	
		pOldBIH->biWidth++;
	}
	pNewBIH->biWidth = iNewWidth;
	while (pNewBIH->biWidth%4 !=0) //再把新图像的宽度信息补为4的倍数
	{							   	
		pNewBIH->biWidth++;
	}
	pNewBIH->biSizeImage = iNewHeight*iNewWidth;
	pNewBFH->bfSize = sizeof(BITMAPFILEHEADER)+sizeof(BITMAPINFOHEADER)+
		nNumColors*sizeof(RGBQUAD)+pNewBIH->biSizeImage;    
	
	
    /*******************************Fast Furiour Translate****************************/
	 
	    
		COMPLEX *piTemplate;
		COMPLEX *templevalue;
		COMPLEX *pTemplatevalue;
		double *ptemple;
		int i,j;
		nNewHeight=iNewHeight;
		nNewWidth=iNewWidth;
		pTemplatevalue=new COMPLEX[iNewHeight*iNewWidth];
		templevalue=new COMPLEX[iNewHeight*iNewWidth];
		ptemple=new double[iNewHeight*iNewWidth];
		
		if(iNewWidth>iNewHeight)
		{
			
			piTemplate=new COMPLEX[iNewWidth];
		}
		else
		{
			
			piTemplate=new COMPLEX[iNewHeight];
		}
		
		for(i=0;i<iNewHeight;i++)//补零且坐标轴移中
		{
					for(j=0;j<iNewWidth;j++)
					{
						if(m_filtertype==1)//傅立叶变换
						{
						
						    if(j<iWidth && i<iHeight)
							{
							     pTemplatevalue[(iNewHeight-1-j)*nNewWidthBytes+i].x=pow(-1,i+j)*pOldBits[(iHeight-1-j)*nWidthBytes+i];
						       	 pTemplatevalue[(iNewHeight-1-j)*nNewWidthBytes+i].y=0;
							}
						    else
							{
		                         pTemplatevalue[(iNewHeight-1-j)*nNewWidthBytes+i].x=0;
							     pTemplatevalue[(iNewHeight-1-j)*nNewWidthBytes+i].y=0;
							}
						}
						else//同态滤波
						{
                            if(j<iWidth && i<iHeight)
							{
								pTemplatevalue[(iNewHeight-1-j)*nNewWidthBytes+i].x=pow(-1,i+j)*log10(pOldBits[(iHeight-1-j)*nWidthBytes+i])/log10(2);
								pTemplatevalue[(iNewHeight-1-j)*nNewWidthBytes+i].y=0;
							}
							else
							{
								pTemplatevalue[(iNewHeight-1-j)*nNewWidthBytes+i].x=0;
								pTemplatevalue[(iNewHeight-1-j)*nNewWidthBytes+i].y=0;
							}
						}
					}
		}
		
	    switch(m_pImageObject->GetNumBits())
		{
	    case 8 :
			{
			if(m_translatetype==1)
			{
			 FastFTValue=new COMPLEX[iNewHeight*iNewWidth];	
			 for(i=0;i<iNewWidth;i++)//逐列傅立叶变换
			 {
				for(j=0;j<iNewHeight;j++)
				{
	                
				    piTemplate[j].x=pTemplatevalue[(iNewHeight-1-j)*nNewWidthBytes+i].x;
					piTemplate[j].y=pTemplatevalue[(iNewHeight-1-j)*nNewWidthBytes+i].y;
					
				}
				
				ReverceOrder(piTemplate,iNewHeight);//每列的元素倒序
				
				FastFuriourTranslate(piTemplate,iNewHeight);//每列傅立叶变换
                for(j=0;j<iNewHeight;j++)
				{
					templevalue[(iNewHeight-1-j)*nNewWidthBytes+i].x=piTemplate[j].x;
					templevalue[(iNewHeight-1-j)*nNewWidthBytes+i].y=piTemplate[j].y;
					

				}
			 }
			 for(j=0;j<iNewHeight;j++)//逐行傅立叶便变换
			 {
				for(i=0;i<iNewWidth;i++)
				{
					
					
					piTemplate[i].x=templevalue[(iNewHeight-1-j)*nNewWidthBytes+i].x;
                    piTemplate[i].y=templevalue[(iNewHeight-1-j)*nNewWidthBytes+i].y;
				}
				ReverceOrder(piTemplate,iNewWidth);//每行的元素倒序
				
				
				FastFuriourTranslate(piTemplate,iNewWidth);//每行傅立叶变换
				
				for(i=0;i<iNewWidth;i++)
				{
					FastFTValue[(iNewHeight-1-j)*nNewWidthBytes+i].x=piTemplate[i].x/(double)iNewWidth;
					FastFTValue[(iNewHeight-1-j)*nNewWidthBytes+i].y=piTemplate[i].y/(double)iNewHeight;
					ptemple[(iNewHeight-1-j)*nNewWidthBytes+i]=log10(1+sqrt(pow(piTemplate[i].x,2)+pow(piTemplate[i].y,2)));
				}
			
			
			 }
		     double max=0,min=255;//求最大值和最小值
			 for(j=0;j<iNewHeight;j++)
			 {
				for(i=0;i<iNewWidth;i++)
				{	
					
				    if(max<ptemple[(iNewHeight-1-j)*nNewWidthBytes+i])
					{
					    max=ptemple[(iNewHeight-1-j)*nNewWidthBytes+i];
					}
					if(min>ptemple[(iNewHeight-1-j)*nNewWidthBytes+i])
					{
						min=ptemple[(iNewHeight-1-j)*nNewWidthBytes+i];
					}
					
				}
			 }
			 for(j=0;j<iNewHeight;j++)//对新图像赋值
			 {
				for(i=0;i<iNewWidth;i++)
				{
				pNewBits[(iNewHeight-1-j)*nNewWidthBytes+i]=(int)((255*((double)(ptemple[(iNewHeight-1-j)*nNewWidthBytes+i]-min)/(double)(max-min)))+0.5);
				}
			 }
			
			}
			if (m_translatetype==2)//傅立叶逆变换
			{
				for(i = 0; i<iNewWidth; i++)
				{
					for(j = 0; j<iNewHeight; j++)
					{
						piTemplate[j].x = FastFTValue[(iNewHeight-1-j)*nNewWidthBytes + i].x;
						piTemplate[j].y = -FastFTValue[(iNewHeight-1-j)*nNewWidthBytes + i].y;
					}
					ReverceOrder(piTemplate,iNewHeight);
					FastFuriourTranslate(piTemplate,iNewHeight);
					
					for(j = 0; j<nNewHeight; j++)
					{
						FastFTValue[(nNewHeight-1-j)*nNewWidthBytes + i].x = piTemplate[j].x;
						FastFTValue[(nNewHeight-1-j)*nNewWidthBytes + i].y = piTemplate[j].y;
					}
				}
				//对行进行傅立叶变换
				
				for(j = 0; j<nNewHeight; j++)
				{
					for(i = 0; i<nNewWidth; i++)
						
					{
						piTemplate[i].x = FastFTValue[(nNewHeight-1-j)*nNewWidthBytes + i].x;
						piTemplate[i].y = FastFTValue[(nNewHeight-1-j)*nNewWidthBytes + i].y;
					}
					ReverceOrder(piTemplate,iNewWidth);
					FastFuriourTranslate(piTemplate,iNewWidth);
					
					for(i = 0; i<nNewWidth; i++)
					{
						FastFTValue[(nNewHeight-1-j)*nNewWidthBytes + i].x = piTemplate[i].x;
						FastFTValue[(nNewHeight-1-j)*nNewWidthBytes + i].y = piTemplate[i].y;
					}
				}
				
				
				//再一次移中
				for(j = 0; j<nNewHeight; j++)
				{
					for(i = 0; i<nNewWidth; i++)
					{   
						
						
						FastFTValue[(nNewHeight-1-j)*nNewWidthBytes + i].x = pow(-1,i+j)*FastFTValue[(iNewHeight-1-j)*nNewWidthBytes + i].x/(double)iNewWidth;
					}
				}				
				
				double max=0,min=255;//求最大值和最小值
				
				for(i = 0; i<nNewWidth; i++)
				{
					for(j = 0; j<nNewHeight; j++)
					{
						ptemple[(nNewHeight-1-j)*nNewWidthBytes + i] = FastFTValue[(nNewHeight-1-j)*nNewWidthBytes + i].x;
						if(ptemple[(nNewHeight-1-j)*nNewWidthBytes + i]<min)
							min = ptemple[(nNewHeight-1-j)*nNewWidthBytes + i];
						if( ptemple[(nNewHeight-1-j)*nNewWidthBytes + i]>max)
							max = ptemple[(nNewHeight-1-j)*nNewWidthBytes + i];
					}
				}
				
				for (j=0; j<nNewHeight; j++)
				{
					for (i=0; i<nNewWidth; i++)
					{
						if(m_filtertype==1)//傅立叶变换
						{
						
						  pNewBits[(nNewHeight-1-j)*nNewWidthBytes+i] = 
							(int)((255* ptemple[(nNewHeight-1-j)*nNewWidthBytes+i]-255*min)/(max-min) + 0.5);
					
						}
						else//同态滤波
						{
                          pNewBits[(nNewHeight-1-j)*nNewWidthBytes+i] =
                            (int)((255* exp(ptemple[(nNewHeight-1-j)*nNewWidthBytes+i])-255*exp(min))/(exp(max)-exp(min)) + 0.5);
						}
					}
					
					
				}
			}
			}
			break;

	    
		case 16:
			
			break;
	    
	    }
	
	//Free the memory for the old image
	::GlobalUnlock(m_pImageObject->GetDib());
	::GlobalFree(m_pImageObject->GetDib());
	
	//Update the DIB in the object pointed by m_pImaeObject
	::GlobalUnlock(hNewDib); 
	m_pImageObject->SetDib(hNewDib);
	return TRUE;