数字图像处理实验报告

数字图像处理实验报告

实验一数字图像基本操作及灰度调整

一、实验目的

1)掌握读、写图像的基本方法。

2)掌握MATLAB语言中图像数据与信息的读取方法。

3)理解图像灰度变换处理在图像增强的作用。

4)掌握绘制灰度直方图的方法，理解灰度直方图的灰度变换及均衡化的方法。

二、实验内容与要求

1.熟悉MATLAB语言中对图像数据读取，显示等基本函数

特别需要熟悉下列命令：熟悉imread()函数、imwrite()函数、size()函数、Subplot （）函数、Figure（）函数。

1)将MA TLAB目录下work文件夹中的forest.tif图像文件读出.用到imread，imfinfo 等文件，观察一下图像数据，了解一下数字图像在MATLAB中的处理就是处理一个矩阵。将这个图像显示出来（用imshow）。尝试修改map颜色矩阵的值，再将图像显示出来，观察图像颜色的变化。

2)将MA TLAB目录下work文件夹中的b747.jpg图像文件读出，用rgb2gray()将其转化为灰度图像，记为变量B。

2.图像灰度变换处理在图像增强的作用

读入不同情况的图像，请自己编程和调用Matlab函数用常用灰度变换函数对输入图像进行灰度变换，比较相应的处理效果。

3.绘制图像灰度直方图的方法，对图像进行均衡化处理

请自己编程和调用Matlab函数完成如下实验。

1)显示B的图像及灰度直方图，可以发现其灰度值集中在一段区域，用imadjust函数将它的灰度值调整到[0，1]之间，并观察调整后的图像与原图像的差别，调整后的灰度直方图与原灰度直方图的区别。

2)对B进行直方图均衡化处理，试比较与源图的异同。

3) 对B 进行如图所示的分段线形变换处理，试比较与直方图均衡化处理的异同。

图1.1 分段线性变换函数

三、实验原理与算法分析

1. 灰度变换

灰度变换是图像增强的一种重要手段，它常用于改变图象的灰度范围及分布，是图象数字化及图象显示的重要工具。

1) 图像反转

灰度级范围为[0, L-1]的图像反转可由下式获得

r L s --=1

2) 对数运算：有时原图的动态范围太大，超出某些显示设备的允许动态范围，如直接

使用原图，则一部分细节可能丢失。解决的方法是对原图进行灰度压缩，如对数变换：

s = c log(1 + r )，c 为常数，r ≥ 0

3) 幂次变换：

0,0,≥≥=γγc cr s

4) 对比拉伸：在实际应用中,为了突出图像中感兴趣的研究对象，常常要求局部扩展

拉伸某一范围的灰度值，或对不同范围的灰度值进行不同的拉伸处理，即分段线性拉伸：

其对应的数学表达式为：

2. 直方图均衡化

灰度直方图的横坐标是灰度级，纵坐标是该灰度级出现的频度，它是图像最基本的统计特征。依据定义，在离散形式下， 用r k 代表离散灰度级，用p r (r k )代表p r (r )，并且有下式成立：

n

n r P k

k r =

)( 1,,2,1,010-=≤≤l k r k 式中：n k 为图像中出现r k 级灰度的像素数，n 是图像像素总数，而n k /n 即为频数。 直方图均衡化处理是以累积分布函数变换法为基础的直方图修正法。假定变换函数为

ωωd p r T s r r

)()(0

?==

(a) Lena 图像 (b) Lena 图像的直方图

图1.2 Lena 图像及直方图

当灰度级是离散值时，可用频数近似代替概率值，即

1,,1,010)(-=≤≤=

l k r n

n r p k k k r

式中：l 是灰度级的总数目，p r (r k )是取第k 级灰度值的概率，n k 是图像中出现第k 级灰度的次数，n 是图像中像素总数。

所以积分可以表示为下列累计分布函数(cumulative distribution function, CDF)

1,,1,01

0)

()(0

-=≤≤===∑∑

==l k r r p n

n r T s j k

j j r k

j j k k

四、实验步骤

1. 熟悉MATLAB 语言中对图像数据读取，显示等基本函数 1) 文件读取与信息显示：

load trees;

[X,map]=imread('forest.tif');

subimage(X,map);

I=imread('forest.tif');

imshow(I);

imfinfo('forest.tif');

2)map颜色矩阵的修改

[X,map]=imread('forest.tif');

map1=map+map;

subimage(X,map1);

3)灰度图像的转化

RGB=imread('b747.jpg');

B=rgb2gray(RGB);

2.图像灰度变换处理在图像增强的作用

g1=imadjust(I,[0 1],[1 0]);

g2=imcomplement(g1);

g3=im2uint8(mat2gray(log(1+double(I))));

3.绘制图像灰度直方图的方法，对图像进行均衡化处理

1)图像灰度直方图的显示与灰度调整

imhist(B);

J = imadjust(B,[],[0 1]);

imhist(J);

subimage(J);

2)对B进行直方图均衡化处理，试比较与原图的异同。

I = imread('pout.tif');

[J,T] = histeq(I);

figure,plot((0:255)/255,T);

3)对B进行如图所示的分段线形变换处理，试比较与直方图均衡化处理的异同。x1=0:0.01:0.125;

x2=0.125:0.01:0.75;

x3=0.75:0.01:1;

y1=2*x1;

y2=0.25+0.6*(x2-0.125);

y3=0.625+1.5*(x3-0.75);

x=[x1,x2,x3];

y=[y1,y2,y3];

subplot(2,2,4);

plot(x,y);

五、实验结果分析与讨论

1.熟悉MATLAB语言中对图像数据读取，显示等基本函数

1)图像文件的读出与图像数据的观察

图1.3 真彩色图像与灰度图像显示ans =

Filename:

'F:\MATLAB\R2007a\toolbox\images\imdemos\forest.tif' FileModDate: '04-Dec-2000 13:57:58' FileSize: 124888

Format: 'tif'

FormatVersion: []

Width: 447

Height: 301

BitDepth: 8

ColorType: 'indexed'

FormatSignature: [73 73 42 0]

ByteOrder: 'little-endian'

NewSubFileType: 0

BitsPerSample: 8

Compression: 'PackBits'

PhotometricInterpretation: 'RGB Palette'

StripOffsets: [17x1 double]

SamplesPerPixel: 1

RowsPerStrip: 18

StripByteCounts: [17x1 double]

XResolution: 72

YResolution: 72

ResolutionUnit: 'Inch'

Colormap: [256x3 double]

PlanarConfiguration: 'Chunky'

TileWidth: []

TileLength: []

TileOffsets: []

TileByteCounts: []

Orientation: 1

FillOrder: 1

GrayResponseUnit: 0.0100

MaxSampleValue: 255

MinSampleValue: 0

Thresholding: 1

ImageDescription: 'Carmanah Ancient Forest, British Columbia, Canada' map颜色矩阵修改后图像颜色的变化

图1.4 原图像与map矩阵值增强一倍后的图像

新的颜色矩阵值变成原文件的2倍，图像明显变亮，颜色的R、G、B值增强。

2)图像文件转化为灰度图像

图1.5 真彩色图像与转化成的灰度图像

2.图像灰度变换处理在图像增强的作用

图1.6 灰度变化增强（图像反转、求补、对数变换）

图中对图像文件进行了基本的灰度变换，包括用式s = L – 1 – r 得到的图像反转，对反转图像的求补，以及对数变换的采用。

3.绘制图像灰度直方图的方法，对图像进行均衡化处理

1)图像灰度及灰度直方图的调整

图1.7 灰度范围的调整与直方图显示

在原始图像中，直方图的组成成分集中在高灰度等级（亮）一侧，且图像灰度范围为[0,1]，故将灰度值调整到[0,1]间后直方图无明显变化。类似的，将灰度值调整到[0,0.5]时，整个图像变暗，直方图横向压缩1倍。

2)直方图均衡化

图1.8 原图像与灰度直方图均衡化结果

原图像中目标物的灰度主要集中于高亮度部分，而且象素总数较多，所占的灰度等级较少。经过直方图均衡后，目标物的所占的灰度等级得到扩展，对比度加强，使整个图像得到增强。

数字图像均衡化后，其直方图并非完全均匀分布，这是因为图像的象素个数和灰度等级均为离散值，而且均衡化使灰度级并归。因此，均衡化后，其直方图并非完全均匀分布。

3)对B进行如图所示的分段线性变换处理

图1.9 原图像与灰度直方图均衡化结果

图1.10 线性变换函数图

通过在所关心范围内为所有灰度值指定一个较高值，而为其他灰度指定一个较低值，或将所需范围变亮，分段线性变换可提高图像中特定灰度范围的亮度，常用于图像特征值的提取。

这里将原始图像位于[0.125,0.75]间的灰度值调低，放大其余的灰度值，突出显示图像低频域和高频域的部分。

六、参考文献

[美] Rafael C.Gonzalez.数字图像处理（第二版）[M].阮秋琦阮宇智，译.

北京：电子工业出版社，2003.3.

实验二 数字图像的空间域滤波和频域滤波

一、实验目的

1. 掌握图像滤波的基本定义及目的。

2. 理解空间域滤波的基本原理及方法。

3. 掌握进行图像的空域滤波的方法。

4. 掌握傅立叶变换及逆变换的基本原理方法。

5. 理解频域滤波的基本原理及方法。

6. 掌握进行图像的频域滤波的方法。

二、实验内容与要求

1. 平滑空间滤波：

1) 读出eight.tif 这幅图像，给这幅图像分别加入椒盐噪声和高斯噪声后并与前一张图 显示在同一图像窗口中。

2) 对加入噪声图像选用不同的平滑（低通）模板做运算，对比不同模板所形成的效果， 要求在同一窗口中显示。

3) 使用函数imfilter 时，分别采用不同的填充方法（或边界选项，如零填 充、’replicate ’、’symmetric ’、’circular ’）进行低通滤波，显示处理后的图像。

4) 运用for 循环，将加有椒盐噪声的图像进行10次，20次均值滤波，查看其特点,显 示均值处理后的图像（提示:利用fspecial 函数的’average’类型生成均值滤波器）。

5) 对加入椒盐噪声的图像分别采用均值滤波法，和中值滤波法对有噪声的图像做处 理，要求在同一窗口中显示结果。

6) 自己设计平滑空间滤波器，并将其对噪声图像进行处理，显示处理后的图像。 2. 锐化空间滤波

1) 读出blurry_moon.tif 这幅图像，采用3×3的拉普拉斯算子w = [ 1, 1, 1; 1 – 8 1; 1, 1, 1]对其进行滤波。

2) 编写函数w = genlaplacian(n)，自动产生任一奇数尺寸n 的拉普拉斯算子，如5×5 的拉普拉斯算子

w = [ 1 1 1 1 1 1 1 1 1 1 1 1 -24 1 1 1 1 1 1 1 1 1 1 1 1]

3) 分别采用5×5，9×9，15×15和25×25大小的拉普拉斯算子对blurry_moon.tif 进

行锐化滤波，并利用式2

(,)(,)(,)g x y f x y f x y =-?完成图像的锐化增强，观察其有何不同，要求在同一窗口中显示。

4) 采用不同的梯度算子对blurry_moon.tif 进行锐化滤波，并比较其效果。

5) 自己设计锐化空间滤波器，并将其对噪声图像进行处理，显示处理后的图像；

3.傅立叶变换

1)读出woman.tif这幅图像，对其进行快速傅立叶变换，分别显示其幅度图像和相位图像。

仅对相位部分进行傅立叶反变换后查看结果图像。

2)仅对幅度部分进行傅立叶反变换后查看结果图像。

3)将图像的傅立叶变换F置为其共轭后进行反变换，比较新生成图像与原始图像的差异。

4.平滑频域滤波

1)设计理想低通滤波器、巴特沃斯低通滤波器和高斯低通滤波器，截至频率自选，分别给出各种滤波器的透视图。

2)读出test_pattern.tif这幅图像，分别采用理想低通滤波器、巴特沃斯低通滤波器和高斯低通滤波器对其进行滤波（截至频率自选），再做反变换，观察不同的截止频率下采用不同低通滤波器得到的图像与原图像的区别，特别注意振铃效应。(提示:1)在频率域滤波同样要注意到填充问题；2）注意到(-1)x+y;)

5.锐化频域滤波

1)设计理想高通滤波器、巴特沃斯高通滤波器和高斯高通滤波器，截至频率自选，分别给出各种滤波器的透视图。

2)读出test_pattern.tif这幅图像，分别采用理想高通滤波器、巴特沃斯高通滤波器和高斯高通滤波器对其进行滤波（截至频率自选），再做反变换，观察不同的截止频率下采用不同高通滤波器得到的图像与原图像的区别。

三、实验原理与算法分析

1.空间域增强

空间域滤波是在图像空间中借助模板对图像进行领域操作，处理图像每一个像素的取值都是根据模板对输入像素相应领域内的像素值进行计算得到的。空域滤波基本上是让图像在频域空间内某个范围的分量受到抑制，同时保证其他分量不变，从而改变输出图像的频率分布，达到增强图像的目的。

空域滤波一般分为线性滤波和非线性滤波两类。线性滤波器的设计常基于对傅立叶变换的分析，非线性空域滤波器则一般直接对领域进行操作。各种空域滤波器根据功能主要分为平滑滤波器和锐化滤波器。平滑可用低通来实现，平滑的目的可分为两类：一类是模糊，目的是在提取较大的目标前去除太小的细节或将目标内的小肩端连接起来；另一类是消除噪声。锐化可用高通滤波来实现，锐化的目的是为了增强被模糊的细节。结合这两种分类方法，可将空间滤波增强分为四类：

线性平滑滤波器（低通）

非线性平滑滤波器（低通）

线性锐化滤波器（高通）

非线性锐化滤波器（高通）

空间滤波器都是基于模板卷积，其主要工作步骤是：

1)将模板在图中移动，并将模板中心与图中某个像素位置重合；

2)将模板上的系数与模板下对应的像素相乘；

3)将所有乘积相加；

4)将和（模板的输出响应）赋给图中对应模板中心位置的像素。

2.平滑滤波器

1)线性平滑滤波器

线性低通平滑滤波器也称为均值滤波器，这种滤波器的所有系数都是正数，对3×3

的模板来说，最简单的是取所有系数为1，为了保持输出图像任然在原来图像的灰度值范围内，模板与象素邻域的乘积都要除以9。

MATLAB 提供了fspecial 函数生成滤波时所用的模板，并提供filter2 函数用指定的滤波器模板对图像进行运算。函数fspecial 的语法格式为：

h=fspecial(type);

h=fspecial(type,parameters);

其中参数type 指定滤波器的种类，parameters 是与滤波器种类有关的具体参数。

表2.1 MATLAB中预定义的滤波器种类

MATLAB 提供了一个函数imnoise 来给图像增添噪声，其语法格式为：

J=imnoise(I,type);

J=imnoise(I,type,parameters);

参数type 指定噪声的种类，parameters 是与噪声种类有关的具体参数。参数的种类见表2.2。

表2.2 噪声种类及参数说明

2)非线性平滑滤波器

中值滤波器是一种常用的非线性平滑滤波器，其滤波原理与均值滤波器方法类似，但计算的非加权求和，而是把领域中的图像的象素按灰度级进行排序，然后选择改组的中间值作为输出象素值。

MATLAB 提供了medfilt2 函数来实现中值滤波，其语法格式为：

B=medfilt2(A,[m n]);

B=medfilt2(A);

其中，A 是原图象，B 是中值滤波后输出的图像。[m n]指定滤波模板的大小，默认模板为3×3。

3.锐化滤波器

图像平滑往往使图像中的边界、轮廓变得模糊，为了减少这类不利效果的影响，需要利用图像锐化技术，使图像的边缘变得清晰。

1)线性锐化滤波器

线性高通滤波器是最常用的线性锐化滤波器。这种滤波器的中心系数都是正的，而周围的系数都是负的，所有的系数之和为0。

对3×3 的模板来说，典型的系数取值为：

[-1 -1 -1;

-1 8 -1;

-1 -1 -1]

事实上这是拉普拉斯算子。语句h=-fspecial(‘laplacian’,0.5)得到的拉普拉斯算子为：

h =-0.3333 -0.3333 -0.3333

-0.3333 2.6667 -0.3333

-0.3333 -0.3333 -0.3333

2)非线性锐化滤波

邻域平均可以模糊图像，因为平均对应积分，所以利用微分可以锐化图像。图像处理中最常用的微分方法是利用梯度。常用的空域非线性锐化滤波微分算子有sobel 算子、prewitt 算子、log 算子等。

4.频域增强

频域增强是利用图像变换方法将原来的图像空间中的图像以某种形式转换到其他空间中，然后利用该空间的特有性质方便地进行图像处理，最后再转换回原来的图像空间中，从而得到处理后的图像。

频域增强的主要步骤是：

选择变换方法，将输入图像变换到频域空间。

在频域空间中，根据处理目的设计一个转移函数，并进行处理。

将所得结果用反变换得到增强的图像。

常用的频域增强方法有低通滤波和高通滤波。

5.低通滤波

图像的能量大部分集中在幅度谱的低频和中频部分，而图像的边缘和噪声对应于高频部分。因此能降低高频成分幅度的滤波器就能减弱噪声的影响。由卷积定理，在频域实现低通滤波的数学表达式：

G(u,v) =H(u,v)F(u,v)

1)理想低通滤波器（ILPF）

2)巴特沃斯低通滤波器（BLPF）

3)指数型低通滤波器（ELPF）

6.高通滤波

由于图像中的细节部分与其高频分量相对应，所以高通滤波可以对图像进行锐化处理。高通滤波与低通滤波相反，它是高频分量顺利通过，使低频分量受到削弱。高通滤波器和低通滤波器相似，其转移函数分别为:

1)理想高通滤波器（IHPF）

2)巴特沃斯高通滤波器（BLPF）

3)指数型高通滤波器（ELPF）

图像经过高通滤波处理后，会丢失许多低频信息，所以图像的平滑区基本上会消失。所以，可以采用高频加强滤波来弥补。高频加强滤波就是在设计滤波传递函数时，加上一

个大于0小于1的常数c，即：

H′(u,v) =H(u,v)+c

四、实验步骤

1.平滑空间滤波：

1)椒盐噪声和高斯噪声的加入

I=imread('eight.tif');

imshow(I);

J = imnoise(I,'salt & pepper',0.05); %noise density=0.05

K= imnoise(I,'gaussian',0.01,0.01);

2)不同的平滑（低通）模板的处理

H = fspecial('sobel');

Sobel = imfilter(I,H,'replicate');

H = fspecial('laplacian',0.4);

lap = imfilter(I,H,'replicate');

H = fspecial('gaussian',[3 3],0.5);

gaussian = imfilter(I,H,'replicate');

3)不同填充方法的滤波

originalRGB = imread('peppers.png');

h = fspecial('motion', 50, 45); %motion blurred

filteredRGB = imfilter(originalRGB, h); boundaryReplicateRGB = imfilter(originalRGB, h, 'replicate'); boundary0RGB = imfilter(originalRGB, h, 'x');

boundary0RGB = imfilter(originalRGB, h, 0); boundarysymmetricRGB = imfilter(originalRGB, h, 'symmetric'); boundarycircularRGB = imfilter(originalRGB, h, 'circular');

4)多次均值滤波

J = imnoise(I,'salt & pepper',0.05);

h=fspecial('average'); %Averaging Filtering

J1=imfilter(J,h);

for i=1:10

J1=imfilter(J,h);

end

for i=1:20

J2=imfilter(J,h);

5)均值、中值滤波

h=fspecial('average'); %Averaging Filtering

J1=imfilter(J,h);

J2=medfilt2(J); %Median Filtering

6)自行设计平滑空间滤波器

domain=[0 0 8 0 0;

0 0 8 0 0;

8 8 8 8 8;

0 0 8 0 0;

0 0 8 0 0];

K1= ordfilt2(J,5,domain);

2.平滑空间滤波：

1)3×3的拉普拉斯算子滤波

I=imread('blurry_moon.tif');

T=double(I);

subplot(1,2,1),imshow(T,[]);title('Original Image'); w =[1,1,1;

1,-8,1;

1,1,1];

K=conv2(T,w,'same');

2)奇数尺寸拉普拉斯算子随机产生函数

function w = genlaplacian(n)

%Computes the Laplacian operator

w = ones(n);

x = ceil(n/2);

w(x, x) = -1 * (n * n - 1);

3)不同尺寸拉普拉斯算子滤波以及图像增强

w1 = genlaplacian(5);

I=imread('blurry_moon.tif');

T=double(I);

K=conv2(T,w1,'same');

J=T-K;

4)不同尺寸梯度算子的锐化滤波

[I,map]=imread('blurry_moon.tif');

I=double(I);

[Gx,Gy]=gradient(I); % gradient calculation G=sqrt(Gx.*Gx+Gy.*Gy); % matrix

J1=G; % gradient1 imshow(J1,map);

J2=I; % gradient2

K=find(G>=7);

J2(K)=G(K);

imshow(J2,map);

J3=I; % gradient3

K=find(G>=7);

J3(K)=255;

imshow(J3,map);

J4=I; % gradient4

K=find(G<=7);

J4(K)=255;

imshow(J4,map);

J5=I; % gradient5

K=find(G<=7);

J5(K)=0;

Q=find(G>=7);

J5(Q)=255;

imshow(J5,map);

5)自行设计锐化空间滤波器

domain=[8 8 0 8 8;

8 8 0 8 8;

0 0 0 0 0;

8 8 0 8 8;

8 8 0 8 8];

K1= ordfilt2(J,5,domain);

3.傅立叶变换

1)图像的快速傅立叶变换，分别显示其幅度图像和相位图像。

F=imread('woman.tif');

F1=fft2(F);

F2=log(1+abs(F1)); %amplitude spectrum

F3=fftshift(F1);

imshow(log(1+abs(F3)),[]);

F4=angle(F1); %phase spectrum

2)相位部分进行傅立叶反变换。

F1=fft2(F);

i=sqrt(-1);

f2=ifft2(exp(i*angle(F1)));

imshow(real(f2),[]);

3)幅度部分进行傅立叶反变换。

f1=ifft2(abs(F1));

imshow(log(1+abs(f1)),[]);

4)将图像的傅立叶变换F置为其共轭后进行反变换

F1=fft2(F);

F2=log(1+abs(F1)); %amplitude spectrum

F3=fftshift(F1);

F4=angle(F1); %phase spectrum

F5=-F4

F6= double(F3*exp(F4)); %the complex conjugate of the fourier transform F7=ifft2(F6); %inverse fourier transform

imshow(real(F7),[]);

4.平滑频域滤波

理想低通滤波器、巴特沃斯低通滤波器和高斯低通滤波器的设计与滤波

理想低通滤波器

I=imread('test_pattern.tif');

f=double(I); % chage into double as MATLAB doesn’t suppor calculation

% of image in unsigned int type

g=fft2(f); % fourier transform

g=fftshift(g); % zero-frequency area centralized

[M,N]=size(g);

d0=100; %cutoff frequency

m=fix(M/2); n=fix(N/2);

for i=1:M

for j=1:N

d=sqrt((i-m)^2+(j-n)^2);

if(d<=d0)

h=1;

else h=0;

end

result(i,j)=h*g(i,j);

end

end

result=ifftshift(result);

J1=ifft2(result);

J2=uint8(real(J1));

imshow(J2)

巴特沃斯低通滤波器（二阶）

I=imread('test_pattern.tif');

f=double(I);

g=fft2(f);

g=fftshift(g);

[M,N]=size(g);

nn=2; % 2-grade Butterworth lowpss filter

d0=100;

m=fix(M/2); n=fix(N/2);

for i=1:M

for j=1:N

d=sqrt((i-m)^2+(j-n)^2);

h=1/(1+0.414*(d/d0)^(2*nn)); % filter transform function

%h=1./(1+(d./d0).^(2*n))

%h=exp(-(d.^2)./(2*(d0^2)));

result(i,j)=h*g(i,j);

end

end

result=ifftshift(result);

J1=ifft2(result);

J2=uint8(real(J1));

imshow(J2);

高斯低通滤波

I=imread('test_pattern.tif');

f=double(I);

g=fft2(f);

g=fftshift(g);

[M,N]=size(g);

d0=100;

m=fix(M/2); n=fix(N/2);

for i=1:M

for j=1:N

d=sqrt((i-m)^2+(j-n)^2);

h=exp(-(d.^2)./(2*(d0^2))); % gaussian filter transform

result(i,j)=h*g(i,j);

end

end

result=ifftshift(result);

J1=ifft2(result);

J2=uint8(real(J1));

5.锐化频域滤波

理想高通滤波器、巴特沃斯高通滤波器和高斯高通滤波器的设计与滤波理想高通滤波器

I=imread('test_pattern.tif');

f=double(I);

g=fft2(f);

g=fftshift(g);

[M,N]=size(g);

d0=80; %d0=15,25,80

m=fix(M/2); n=fix(N/2);

for i=1:M

for j=1:N

d=sqrt((i-m)^2+(j-n)^2);

if(d>=d0)

h=1;

else h=0;

end

result(i,j)=h*g(i,j);

end

end

result=ifftshift(result);

J1=ifft2(result);

J2=uint8(real(J1));

巴特沃斯高通滤波器

I=imread('test_pattern.tif');

f=double(I);

g=fft2(f);

g=fftshift(g);

[M,N]=size(g);

nn=2; % 2-grade nutterworth highpass filter

d0=80;

m=fix(M/2);

n=fix(N/2);

for i=1:M

for j=1:N

d=sqrt((i-m)^2+(j-n)^2);

if (d==0)

h=0;

else

h=1/(1+0.414*(d0/d)^(2*nn));%transform fuction calculation end

result(i,j)=h*g(i,j);

end

end

result=ifftshift(result);

J1=ifft2(result);

J2=uint8(real(J1));

高斯高通滤波器

I=imread('test_pattern.tif');

f=double(I);

g=fft2(f);

g=fftshift(g);

[M,N]=size(g);

d0=80;

m=fix(M/2); n=fix(N/2);

for i=1:M

for j=1:N

d=sqrt((i-m)^2+(j-n)^2);

h=1-exp(-(d.^2)./(2*(d0^2))); % gaussian filter transform

result(i,j)=h*g(i,j);

end

end

result=ifftshift(result);

J1=ifft2(result);

J2=uint8(real(J1));

五、实验结果分析与讨论

1.平滑空间滤波：

1)读出eight.tif这幅图像，给这幅图像分别加入椒盐噪声和高斯噪声后并与前一张图

显示在同一图像窗口中。