testdft.cpp

普通的离散图像信号的傅里叶变换和逆变换。

#include "stdafx.h"
#include "TestDft.h"
#include <math.h>
#include <iostream>
using namespace std;
#define PI 3.1415926
#define T 1
#define A 0.1
#define B 0.1
#define LOW 20
#ifdef _DEBUG
#define new DEBUG_NEW
#endif

Dft::Dft()
{
	m_source = cvLoadImage("d:/23.jpg",0);
	m_size = cvCreateImage(cvGetSize(m_source),IPL_DEPTH_8U,1);
	m_angle = cvCreateImage(cvGetSize(m_source),IPL_DEPTH_8U,1);
	m_outimage = cvCreateImage(cvGetSize(m_source),IPL_DEPTH_8U,1);

	Width = (float)m_source->width;
	Height = (float)m_source->height;
	nums.Greal = new float[Width*Height];
	nums.Gim = new float[Width*Height];
	nums.htempreal = new float[Width*Height];
	nums.htempim = new float[Width*Height];
	nums.Hreal= new float[Width*Height];
	nums.Him = new float[Width*Height];
}

Dft::~Dft()
{
	delete nums.Greal;
	delete nums.Gim;
	delete nums.htempreal;
	delete nums.htempim;
	delete nums.Hreal;
	delete nums.Him;
}
void Dft::DftMethord(IplImage *psource, IplImage *paimsize, IplImage *paimangle)
{
	tempnums.imnum = 0;
	tempnums.realnum = 0;
	nums.realnum = 0;
	nums.imnum = 0;
	float localreal = 0;
	float localim = 0;
	int Step = 0;
	Step = psource->widthStep/sizeof(uchar);
	float picsize =  Width*Height;
	int m,n,i,k; //循环变量
	float cos_value,sin_value;
	float hsinsqrare;
	float hsincos;
	float hcoefficient;      //求逆滤波用
	uchar *sourcedatas = NULL;
	sourcedatas = (uchar *)psource->imageData;
	uchar *aimdatas = NULL;
	cvSetZero(paimsize);
	aimdatas = (uchar*)paimsize->imageData;
	uchar *aimangledatas = NULL;
	aimangledatas = (uchar *)paimangle->imageData;
	float sourcerealgray = 0.0;
	float sourceimgray = 0.0;
	float newgray = 0.0;
	float atan_value = 0;
	float temp = 0;
	int index = 0;
	double data = 0;
	double ivar;
	double jvar;
	int radius = 0;
	for(m=0; m<Width; m++)                          
	{
		for(n=0; n<Height; n++)
		{
			for(i=0; i<Width; i++)
			{
				for(k=0; k<Height; k++)
				{
					sourcerealgray = sourcedatas[k*Step+i];

					cos_value=cos(2*PI*(m*i/Width+n*k/Height));
					sin_value=-sin(2*PI*(m*i/Width+n*k/Height));		       
					if(fmod(/*(float)*/(i+k),2.0) == 0)
					{
						tempnums.realnum = sourcerealgray*cos_value/picsize;
						tempnums.imnum = sourcerealgray*sin_value/picsize;
					}
					else
					{
						tempnums.realnum = (-1)*sourcerealgray*cos_value/picsize;
						tempnums.imnum = (-1)*sourcerealgray*sin_value/picsize;
					}
					nums.realnum +=tempnums.realnum;  //完成书中的实域相加
					nums.imnum += tempnums.imnum;     //虚数相加
				} 
			}
			nums.Greal[n*Step+m] = nums.realnum;//这两行是原图像的傅里叶变换。
			nums.Gim[n*Step+m] = nums.imnum;    //保存实部与虚部


			aimdatas[n*Step+m] = log(sqrt(nums.realnum*nums.realnum + nums.imnum * nums.imnum)); //幅值下面这句保证傅里叶变换正确。
			atan_value = atan(nums.imnum/nums.realnum);    //相角
			aimangledatas[n*Step+m] = (uchar)atan_value;  //给输出图像复值

			nums.realnum = 0;  //为下一个清零
			nums.imnum = 0;
		}
		cout<<m<<endl;
	}//到这里傅里叶变换结束 
	aimangledatas[n*Step+m] = (uchar)atan_value;
	cvNamedWindow("blur",1);  //这里是显示模糊图像的傅里叶变换
	cvShowImage("blur",paimsize);
	cvSaveImage("D:/blur.jpg",paimsize);
	cvNamedWindow("phase",1);
	cvShowImage("phase",paimangle);
	cvSaveImage("D:/blurphase.jpg",paimangle);
	cvWaitKey(0);
	cout<<"over";
}

void Dft::DftConvertMethord()//这个函数是傅里叶反变换
{
	IplImage *pangleimage = cvLoadImage("D:\\blurphase.jpg",0);
	IplImage* pvalueimage;
	IplImage *poutimage;

	poutimage = cvCreateImage(cvGetSize(pangleimage),IPL_DEPTH_8U,1);
	cvSetZero(poutimage);
	tempnums.imnum = 0;
	tempnums.realnum = 0;
	nums.realnum = 0;
	nums.imnum = 0;
	float Width = 0;
	float Height = 0;
	int Step = 0;
	Step = pangleimage->widthStep/sizeof(uchar);
	Width = (float)pangleimage->width;
	Height = (float)pangleimage->height;
	float picsize =  Width*Height;
	int m,n,i,k; //循环变量
	float cos_value,sin_value;

	uchar *aimoutdatas = NULL;
	aimoutdatas = (uchar *)poutimage->imageData;
	float sourcerealgray = 0.0;
	float newgray = 0.0;
	float atan_value;
	double data = 0;
	double ivar;
	double jvar;
	int radius = 0;
	for(i=0; i<Width; i++)                          
	{
		for(k=0; k<Height; k++)
		{
			for(m=0; m<Width; m++)
			{
				for(n=0; n<Height; n++)
				{

					cos_value=cos(2*PI*(i*m/Width+k*n/Height));
					sin_value=sin(2*PI*(i*m/Width+k*n/Height));	


					ivar = (m-Width/2)*(m-Width/2);
					jvar = (n-Height/2)*(n-Height/2);
					data  = ivar + jvar;
					radius= (int)sqrt(data);
					if (radius<= LOW)
					{
						tempnums.realnum = (nums.Hreal[n*Step+m]*cos_value - nums.Him[n*Step+m]*sin_value);
						tempnums.imnum = (nums.Hreal[n*Step+m]*sin_value + nums.Him[n*Step+m]*cos_value);
						nums.realnum +=tempnums.realnum;
						nums.imnum += tempnums.imnum;
					}
				} 
			}
			//atan_value=(uchar)sqrt(realsq+ imsq);如果用幅值作逆变换就不用中心变换
			atan_value = (uchar)sqrt(nums.imnum*nums.imnum +nums.realnum*nums.realnum);
			aimoutdatas[k*Step+i]= atan_value;

			nums.realnum = 0;
			nums.imnum = 0;
		}
		cout<<i<<" "<<atan_value<<endl;
	}
	cvNamedWindow("angleout",1);
	cvShowImage("angleout",poutimage);
	cvSaveImage("d:/phaseout.jpg",poutimage);
	cvWaitKey(0);
}
// 唯一的应用程序对象

CWinApp theApp;

using namespace std;
int _tmain(int argc, TCHAR* argv[], TCHAR* envp[])
{
	int nRetCode = 0;

	// 初始化MFC 并在失败时显示错误
	if (!AfxWinInit(::GetModuleHandle(NULL), NULL, ::GetCommandLine(), 0))
	{
		// TODO: 更改错误代码以符合您的需要
		_tprintf(_T("错误: MFC 初始化失败\n"));
		nRetCode = 1;
	}
	else
	{
		Dft *pDft = new Dft;
		pDft->DftMethord(pDft->m_source,pDft->m_size,pDft->m_angle);
		pDft->DftConvertMethord();
		delete pDft;
	}

	return nRetCode;
}