医学图像处理实验一指导书-2013

3
4）系数 WN 的确定
r
1 。其中每个蝶形 运算中系数的确定步骤为： ① 将每个蝶形运算中第一个节点标号值， 即 p 值，表示成 M
第 m 级蝶形运算中涉及的系数 WN 共有 2
N r
m 1
种，即 r 0,1, ,2
m 1
位（注意 M log 2 ）二进制数； ②将此 M 位二进制数左移 M-m 位，把右边空出的位置补 零，此数即为所求 r 的二进制数。
存储单元
A(1) A (2) x (1) A (3) x (2) A (4) x (3) A (5) x (4) A (6) x (5) A (7) x (6) A (8) x (7)
自然顺序输入 x (0)
变址 倒位序
x (0) x (4) x (2) x (6) x (1) x (5) x (3) x (7)
四、实验组织运行要求
本实验为集中授课模式，根据任务要求，由教师先进行讲解演示，再由学生根据学生自 主完成实验内容，最后进行总结分析。
五、实验条件
计算机、MATLAB 软件、测试图像等。
六、实验步骤
1. 图像读取、显示与信息获取 在 Windows 下打开 MATLAB 软件，在命令窗口输入相关命令，依次实现读取图像、显示 图像及获取图像信息。 2. 编写计算图像 DFT 的程序并完成测试分析 1）新建 M 文件，参考二维 DFT 公式，编写直接计算程序，并调试、保存程序；运行编 写的程序计算图像 square.bmp 的二维 DFT, 保存并显示结果，记录运算时间。 2）新建 M 文件，参考 FFT 算法，编写一维 FFT 计算程序，并调试、保存程序；基于二 维 DFT 的可分离性，调用编写的一维 FFT 程序，计算图像 square.bmp 的二维 DFT, 保存并 显示结果，记录运算时间； 3） 在 MATLAB 的命令窗口输入 MATLAB 提供的二维 FFT 计算函数， 计算图像 square.bmp 的二维 DFT, 保存并显示结果。 4）比较自编程序与 MATLAB 提供函数的运行结果与运行时间，结合相关理论进行对比 分析。 3. 编写中值滤波、均值滤波程序并测试滤波效果 1）新建 M 文件，根据均值滤波和中值滤波原理，分别编写计算程序，并调试、保存 程序；2）调用 MATLAB 提供的高斯噪声、椒盐噪声图像和乘性噪声函数，给图像？叠加三种 不同噪声；3）调用编写的滤波程序对加噪后的图像进行处理，并调用 MATLAB 提供的对 应滤波函数作同样处理，比较结果验证自编程序的正确性；4)给图像叠加不同程度的噪声、 选用不同的滤波参数对图像进行处理， 将处理结果结合理论对两种不同滤波方法、 参数对滤 波效果的影响进行比较分析。
N-1 F(u,v)
x
即：
x
图 1 二位傅里叶变换的可分离性
u
F (u , y ) f ( x, y )e
y 0
N 1
j 2
ux N
， u 0,1, , N 1;
（2） （3）
F (u, v) F (u , y )e
x 0
N 1
j 2
vy N
， v 0,1, , N 1;
实验一：图像的傅里叶变换与空域滤波
实验学时：4 实验类型：设计 实验要求：必修
一、实验目的
通过本实验，掌握数字图像的表示、傅里叶变换的原理、快速傅里叶变换的 实现方法、图像的常见噪声类型及基本的空域滤波方法，在此基础上熟悉 MATLAB 编程语言，熟悉编程方法。
二、实验内容
1. 用 MATLAB 读取、显示图像 square.bmp 和 lena.jpg，并获取图像的全部 信息及矩阵参数； 2. 根据二维离散傅里叶变换（DFT）公式和快速傅里叶变换（FFT）算法分 别编写程序，计算图像 square.bmp 的傅里叶变换；用 MATLAB 提供的二维 DFT 函数作相同的计算，比较计算结果、所需计算时间，并与理论分析作对比； 3.用 MABLAB 提供的函数分别给图像 lena.bmp 加入高斯噪声、椒盐噪声及 乘性噪声，采用均值滤波、中值滤波对加噪声后的图像进行处理，比较不同滤波 方法的处理结果，分析不同参数设置对滤波效果的影响。
图 2 蝶形运算示意图 由图可见，某一列的任何两个节点 p 和 q 的节点变量进行蝶形运算后，得到结果为下一列 p 、 q 两节点的节点变量，而和其他节点变量无关，因而可以采用原位运算，节省存储单 元。
2
2）倒位序规律 为了让 FFT 的输出 X ( p ) 按正常顺序排列在存储单元中，即按 X(0)，X(1)，…，X(N-1) 的顺序排列，输入 x(n)不能按自然顺序存储的，而是按 x(0)，x(4)， …， x(7)的顺序存入存 储单元，称之为倒位序。造成倒位序的原因是输入 x(n)按标号 n 的偶奇的不断分组。 实际运算中，一般先按自然顺序将输入序列存入存储单元， 再通过变址运算得到倒位 N 序排列。将输入序列的自然顺序号 I 用二进制数表示，即 nM 1nM 2 n1n0 （其中 M log 2 为 ，这样，在原来自然 二进制位数） ，则其倒位序 J 对应的二进制数就是（ n0 n1 nM 2 nM 1 ） 顺序时应该放 x(I)的单元，现在倒位序后应放 x(J)。 例如, N=8 时，x(3)的标号是 I 3 ，它 的二进制数是 011，倒位序的二进制数是 110，即 I 6 ，所以原来存放在 x(011)单元的数据 现在应该存放在 x(110)内。表 1 列出了 N=8 时的自然顺序二进制数以及相应的倒位序二进 制数： 表 1 N=8 时的自然顺序二进制数以及相应的倒位序二进制数 自然顺序（I） 0 1 2 3 4 5 6 7 二进制数 000 001 010 011 100 101 110 111 倒位序二进制数 000 100 010 110 001 101 011 111 倒位序（J） 0 4 2 6 1 5 3 7
j 2
ux vy N
u 0,1, , N 1; v 0,1, , N 1;
示
（1）
基于二维 DFT 的可分离性，二维 FFT 算法可用两个一维 FFT 算法来实现，如图 1 所
0 N-1
N-1 f(x,y)
y
行变换 0 N-1
N-1 F(x,v)
v
列变换
0
N-1 v
1
2. 二维 DFT 公式与性质
设 f ( x, y ) ( x 0,1, , N 1; y 0,1, , N 1 ) 是一幅 N N 图像， 其二维离散傅里 叶变换 F (u , v) 定义为：
F (u, v) f ( x, y )e
x 0 y 0
N 1 N 1
图 3 N=8 倒位序的变址处理 3）对偶节点 X (q ) 的确定 当输入为倒位序，输出为正常顺序时，其第 m 级蝶形运算中构成一个蝶形的两节点 间的“距离”为 2m-1，则 X m ( p ) 的对偶节点 X m ( q ) 的标号可确定为：
q p 2 m 1
（6）
如 X 3 (3) 的对偶节点为 X 3 (7) 。
N
的总复乘次数为
N N N 和总复加次数为 N log 2 。 log 2 2 2 直接计算一维 DFT 运算总共需要 N 次复数乘法及 N(N-1)次复数加法。
5. 空域滤波原理
1）基本原理 空域滤波可表示为
g (i, j ) f (i, j ) h(i, j ) [ f (m, n)h(i m, j n)]
把按自然顺序存放在存储单元中的数据，换成 FFT 原位运算流图所要求的倒位序的变 址功能如图 3 所示，当 I=J 时，不必调换，当 I≠J 时， 必须将原来存放数据 x(I)的存储单 元内调入数据 x(J)，而将存放 x(J)的存储单元内调入 x(I)。为了避免把已调换过的数据再次 调换，保证只调换一次（否则又回到原状） ，只需看 J 是否比 I 小。若 J 比 I 小，则意味着 此 x(I)在前边已和 x(J)互相调换过， 不必再调换了; 只有当 J>I 时， 才将原存放 x(I)及存放 x(J) 的存储单元内的内容互换。这样就得到输入所需的倒位序列的顺序。
4
值，即
g (i, j ) Median[ f (i
m 1 n 1 m 1 n 1 m 1 n 1 , j ), f (i 1, j ),, f (i , j )] 2 2 2 2 2 2
（9） ， (i, j ) 为当前模板 其中 m, n 定义模板尺寸，决定了参与排序的像素范围与数目（即邻域） 下的中心像素坐标。
m n
（7）
其中 f ( x, y ) 为输入图像， h( x, y ) 为滤波函数（即模板系数） ，图示如下：
图 4 空域滤波的原理示意图 2）均值滤波 均值滤波公式可表示为
g (i, j )
1 M
( n ,m )s
f (m, n)
（8）
其中 m, n 定义模板尺寸， M m n 为模板系数的数目， (i, j ) 为当前模板下的中心像素。 3）中值滤波 中值滤波将模板覆盖下的像素的数值进行排序， 取中间值作为当前模板下中心像素的数
3. 一维 FFT 算法原理
1）FFT 算法的特点
将一维离散傅里叶变换定义为
F (u ) f ( x)e
x 0
N 1
j 2 ux N
ux ，其中 WN e4）
，每级（每列）计算都包含 FFT 算法由多级蝶形运算构成（级数为 M log 2 ）
4. FFT 算法与直接计算 DFT 运算量的比较
对包含 N 个数据点的一维序列，其 FFT 计算共有 M log 2 级蝶形运算， 每级都由 N/2 个蝶形运算组成，每个蝶形需要一次复乘、 二次复加，因而每级运算都需 N/2 次复乘 和 N 次复加， 从而可计算出总共的复乘和复加次数。 由此可推算出用 FFT 算法计算二维 DFT