实验三-图像的傅立叶变换.doc

实验三、图像的傅立叶变换
一、实验目的
1了解图像变换的意义和手段；
2熟悉傅里叶变换的基本性质；
3熟练掌握FFT的方法及应用；
4通过实验了解二维频谱的分布特点；
5通过本实验掌握编程实现数字图像的傅立叶变换。
二、实验原理
1应用傅立叶变换进行图像处理
傅里叶变换是线性系统分析的一个有力工具，它能够定量地分析诸如数字化系统、采样点、电子放大器、卷积滤波器、噪音和显示点等的作用。通过实验培养这项技能，将有助于解决大多数图像处理问题。对任何想在工作中有效应用数字图像处理技术的人来说，把时间用在学习和掌握博里叶变换上是很有必要的。
2傅立叶（Fourier）变换的定义
　　对于二维信号，二维Fourier变换定义为：

　　二维离散傅立叶变换为：

　　图像的傅立叶变换与一维信号的傅立叶变换变换一样，有快速算法，具体参见参考书目，有关傅立叶变换的快速算法的程序不难找到。实际上，现在有实现傅立叶变换的芯片，可以实时实现傅立叶变换。
三、实验内容
1.根据二维离散Fourier变换的定义编写程序
2.实现图象的变换
3.画出图象的频谱图。
2提交实验的原始图像和结果图像。
3报告内容主要是将原理与实验现象相结合，利用原理简要解释观察到的现象。
六、思考题
1．傅里叶变换有哪些重要的性质?
2．图像的二维频谱在显示和处理时应注意什么?
1. 源程序：
// MyDIPView.cpp : implementation of the CMyDIPView class
//

#include "stdafx.h"
#include "MyDIP.h"

#include "MyDIPDoc.h"
#include "MyDIPView.h"
#include "math.h"
#include<complex>
using namespace std;
#define PI 3.14159
#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif
/*************************************************************************
*
* 函数名称：
* FFT()
*
* 参数:
* complex<double> * TD - 指向时域数组的指针
* complex<double> * FD - 指向频域数组的指针
* r －2的幂数，即迭代次数
*
* 返回值:
* 无。
*
* 说明:
* 该函数用来实现快速付立叶变换。
*
************************************************************************/
VOID WINAPI FFT(complex<double> * TD, complex<double> * FD, int r)
{
// 付立叶变换点数
LONG count;

// 循环变量
int i,j,k;

// 中间变量
int bfsize,p;

// 角度
double angle;

complex<double> *W,*X1,*X2,*X;

// 计算付立叶变换点数
count = 1 << r;

// 分配运算所需存储器
W = new complex<double>[count / 2];
X1 = new complex<double>[count];
X2 = new complex<double>[count];

// 计算加权系数
for(i = 0; i < count / 2; i++)
{
angle = -i * PI * 2 / count;
W[i] = complex<double> (cos(angle), sin(angle));
}

// 将时域点写入X1
memcpy(X1, TD, sizeof(complex<double>) * count);

// 采用蝶形算法进行快速付立叶变换
for(k = 0; k < r; k++)
{
for(j = 0; j < 1 << k; j++)
{
bfsize = 1 << (r-k);
for(i = 0; i < bfsize / 2; i++)
{
p = j * bfsize;
X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2];
X2[i + p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2]) * W[i * (1<<k)];
}
}
X = X1;
X1 = X2;
X2 = X;
}

// 重新排序
for(j = 0; j < count; j++)
{
p = 0;
for(i = 0; i < r; i++)
{
if (j&(1<<i))
{
p+=1<<(r-i-1);
}
}
FD[j]=X1[p];
}

// 释放内存
delete W;
delete X1;
delete X2;
}

/*************************************************************************
*
* 函数名称：
* Fourier()
*
* 参数:
* LPSTR lpDIBBits - 指向源DIB图像指针
* LONG lWidth - 源图像宽度（象素数）
* LONG lHeight - 源图像高度（象素数）
*
* 返回值:
* BOOL - 成功返回TRUE，否则返回FALSE。
*
* 说明:
* 该函数用来对图像进行付立叶变换。
*

BOOL WINAPI Fourier(LPSTR lpDIBBits, LONG lWidth, LONG lHeight)
{

// 指向源图像的指针
unsigned char* lpSrc;
// 中间变量
double dTemp;
// 循环变量
LONG i;
LONG j;
// 进行付立叶变换的宽度和高度（2的整数次方）
LONG w;
LONG h;
int wp;
int hp;

// 图像每行的字节数
LONG lLineBytes;
// 计算图像每行的字节数
lLineBytes = WIDTHBYTES(lWidth * 8);

// 赋初值
w = 1;
h = 1;
wp = 0;
hp = 0;

// 计算进行付立叶变换的宽度和高度（2的整数次方）
while(w * 2 <= lWidth)
{
w *= 2;
wp++;
}

while(h * 2 <= lHeight)
{
h *= 2;
hp++;
}

// 分配内存
complex<double> *TD = new complex<double>[w * h];
complex<double> *FD = new complex<double>[w * h];

// 行
for(i = 0; i < h; i++)
{
// 列
for(j = 0; j < w; j++)
{
// 指向DIB第i行，第j个象素的指针
lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;

// 给时域赋值
TD[j + w * i] = complex<double>(*(lpSrc), 0);
}
}

for(i = 0; i < h; i++)
{
// 对y方向进行快速付立叶变换
FFT(&TD[w * i], &FD[w * i], wp);
}

// 保存变换结果
for(i = 0; i < h; i++)
{
for(j = 0; j < w; j++)
{
TD[i + h * j] = FD[j + w * i];
}
}

for(i = 0; i < w; i++)
{
// 对x方向进行快速付立叶变换
FFT(&TD[i * h], &FD[i * h], hp);
}

// 行
for(i = 0; i < h; i++)
{
// 列
for(j = 0; j < w; j++)
{
// 计算频谱
dTemp = sqrt(FD[j * h + i].real() * FD[j * h + i].real() +
FD[j * h + i].imag() * FD[j * h + i].imag()) / 100;

// 判断是否超过255
if (dTemp > 255)
{
// 对于超过的，直接设置为255
dTemp = 255;
}

// 指向DIB第(i<h/2 ? i+h/2 : i-h/2)行，第(j<w/2 ? j+w/2 : j-w/2)个象素的指针
// 此处不直接取i和j，是为了将变换后的原点移到中心
//lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;
lpSrc = (unsigned char*)lpDIBBits + lLineBytes *
(lHeight - 1 - (i<h/2 ? i+h/2 : i-h/2)) + (j<w/2 ? j+w/2 : j-w/2);

// 更新源图像
* (lpSrc) = (BYTE)(dTemp);
}
}

// 删除临时变量
delete TD;
delete FD;

// 返回
return TRUE;
}

/////////////////////////////////////////////////////////////////////////////
// CMyDIPView

IMPLEMENT_DYNCREATE(CMyDIPView, CScrollView)

BEGIN_MESSAGE_MAP(CMyDIPView, CScrollView)
//{{AFX_MSG_MAP(CMyDIPView)
ON_COMMAND(ID_MENUITEM32778, OnMenuitem32778)
//}}AFX_MSG_MAP
// Standard printing commands
ON_COMMAND(ID_FILE_PRINT, CScrollView::OnFilePrint)
ON_COMMAND(ID_FILE_PRINT_DIRECT, CScrollView::OnFilePrint)
ON_COMMAND(ID_FILE_PRINT_PREVIEW, CScrollView::OnFilePrintPreview)
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CMyDIPView construction/destruction

CMyDIPView::CMyDIPView()
{
// TODO: add construction code here

}

CMyDIPView::~CMyDIPView()
{
}

BOOL CMyDIPView::PreCreateWindow(CREATESTRUCT& cs)
{
// TODO: Modify the Window class or styles here by modifying
// the CREATESTRUCT cs

return CScrollView::PreCreateWindow(cs);
}

/////////////////////////////////////////////////////////////////////////////
// CMyDIPView drawing

void CMyDIPView::OnDraw(CDC* pDC)
{
CMyDIPDoc* pDoc = GetDocument();
ASSERT_VALID(pDoc);
if(pDoc->m_hDIB == NULL)
return ;
// TODO: add draw code for native data here
int i,j;
unsigned char *lpSrc;
LPSTR lpDIB = (LPSTR) ::GlobalLock((HGLOBAL) pDoc->m_hDIB);
int cxDIB = (int) ::DIBWidth(lpDIB); // Size of DIB - x
int cyDIB = (int) ::DIBHeight(lpDIB); // Size of DIB - y
LPSTR lpDIBBits=::FindDIBBits (lpDIB);
// 计算图像每行的字节数
long lLineBytes = WIDTHBYTES(cxDIB * 8);
// 每行
for(i = 0; i < cyDIB; i++)
{
// 每列
for(j = 0; j < cxDIB; j++)
{
// 指向DIB第i行，第j个象素的指针
lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (cyDIB - 1 - i) + j;
// 计算新的灰度值
//*(lpSrc) = BYTE(255-*lpSrc);
}
}
::GlobalUnlock((HGLOBAL) pDoc->m_hDIB);
CRect rect(0,0,cxDIB,cyDIB), rcDIB(0,0,cxDIB,cyDIB);
::PaintDIB(pDC->m_hDC, &rect, pDoc->m_hDIB, &rcDIB, pDoc->m_palDIB);
}

// CMyDIPView printing

BOOL CMyDIPView::OnPreparePrinting(CPrintInfo* pInfo)
{
// default preparation
return DoPreparePrinting(pInfo);
}

void CMyDIPView::OnBeginPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)
{
// TODO: add extra initialization before printing
}

void CMyDIPView::OnEndPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)
{
// TODO: add cleanup after printing
}

/////////////////////////////////////////////////////////////////////////////
// CMyDIPView diagnostics

#ifdef _DEBUG
void CMyDIPView::AssertValid() const
{
CScrollView::AssertValid();
}

void CMyDIPView::Dump(CDumpContext& dc) const
{
CScrollView::Dump(dc);
}

CMyDIPDoc* CMyDIPView::GetDocument() // non-debug version is inline
{
ASSERT(m_pDocument->IsKindOf(RUNTIME_CLASS(CMyDIPDoc)));
return (CMyDIPDoc*)m_pDocument;
}
#endif //_DEBUG

void CMyDIPView::OnMenuitem32778()
{
// 图像付立叶变换

// 获取文档
CMyDIPDoc* pDoc = GetDocument();

// 指向DIB的指针
LPSTR lpDIB;

// 指向DIB象素指针
LPSTR lpDIBBits;

// 锁定DIB
lpDIB = (LPSTR) ::GlobalLock((HGLOBAL) pDoc->GetHDIB());

// 找到DIB图像象素起始位置
lpDIBBits = ::FindDIBBits(lpDIB);

// 判断是否是8-bpp位图（这里为了方便，只处理8-bpp位图的付立叶变换，其它的可以类推）
if (::DIBNumColors(lpDIB) != 256)
{
// 提示用户
MessageBox("目前只支持256色位图的付立叶变换！", "系统提示" ,
MB_ICONINFORMATION | MB_OK);

// 解除锁定
::GlobalUnlock((HGLOBAL) pDoc->GetHDIB());

// 返回
return;
}

// 更改光标形状
BeginWaitCursor();

// 调用Fourier()函数进行付立叶变换
if (::Fourier(lpDIBBits, ::DIBWidth(lpDIB), ::DIBHeight(lpDIB)))
{

// 设置脏标记
pDoc->SetModifiedFlag(TRUE);

// 更新视图
pDoc->UpdateAllViews(NULL);
}
else
{
// 提示用户
MessageBox("分配内存失败！", "系统提示" , MB_ICONINFORMATION | MB_OK);
}

// 解除锁定
::GlobalUnlock((HGLOBAL) pDoc->GetHDIB());

// 恢复光标
EndWaitCursor();

}
2.程序运行结果截图
实验原始图像：

一次傅里叶变换

第二次傅里叶变换

第三次傅里叶变换

3.傅立叶变换的性质及意义
傅立叶变换是数字信号处理领域一种很重要的算法。要知道傅立叶变换算法的意义，首先要了解傅立叶原理的意义。傅立叶原理表明：任何连续测量的时序或信号，都可以表示为不同频率的正弦波信号的无限叠加。而根据该原理创立的傅立叶变换算法利用直接测量到的原始信号，以累加方式来计算该信号中不同正弦波信号的频率、振幅和相位。
图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度。如：大面积的沙漠在图像中是一片灰度变化缓慢的区域，对应的频率值很低；而对于地表属性变换剧烈的边缘区域在图像中是一片灰度变化剧烈的区域，对应的频率值较高。傅立叶变换在实际中有非常明显的物理意义，设f是一个能量有限的模拟信号，则其傅立叶变换就表示f的谱。从纯粹的数学意义上看，傅立叶变换是将一个函数转换为一系列周期函数来处理的。从物理效果看，傅立叶变换是将图像从空间域转换到频率域，其逆变换是将图像从频率域转换到空间域。换句话说，傅立叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数，傅立叶逆变换是将图像的频率分布函数变换为灰度分布函数。
另外：
1、图像经过二维傅立叶变换后，其变换系数矩阵表明：
若变换矩阵Fn原点设在中心，其频谱能量集中分布在变换系数短阵的中心附近(图中阴影区)。若所用的二维傅立叶变换矩阵Fn的原点设在左上角，那么图像信号能量将集中在系数矩阵的四个角上。这是由二维傅立叶变换本身性质决定的。同时也表明一股图像能量集中低频区域。
2 、变换之后的图像在原点平移之前四角是低频，最亮，平移之后中间部分是低频，最亮，亮度大说明低频的能量大（幅角比较大）。
FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外，FFT可以将一个信号的频谱提取出来，这在频谱分析方面也是经常用的。