1、第五章第五章 MATLAB在遥感图像处理中在遥感图像处理中的应用的应用主要内容 遥感图像概述 遥感图像的读写与显示 遥感图像辐射增强 遥感图像几何变换 遥感图像配准 遥感图像滤波 遥感图像分割 1 遥感图像处理概述基本概念遥感数字图像：是指被计算机存储、处理和使用的图像，是一种空间坐标和灰度均不连续的、用数字形式表达的遥感影像，最基本单元是像素，每个像素具有空间位置特征和属性特征。空间位置特征：是用离散的X值和Y值来表示；属性特征：常用亮度值表示。亮度值有如下特点：（1）不同图像相同地点的亮度值不同；（2）亮度值大小由传感器所探测到的电磁辐射强度决定，入射到传感器中的电磁波被探测元件转化为电信
2、号，经过A/D转换，成为绝对辐射亮度值R。为了便于应用R又被转换为能够表征地物的辐射亮度相对值V。R=V*（Rmax-Rmin）/Dmax+Rmin Rmax为探测器可检测到的最大辐射亮度；Rmin为探测器可检测到的最小辐射亮度；max为级数；R为辐射亮度值；V为像素表征的地物辐射亮度的相对值。遥感数字图像处理：利用数字计算机或其它高速、大规模集成数字硬件，对从遥感图像信息转换来的数字电信号进行某些数字运算或处理（如去除噪声、增强、复原、分割、提取特征等），以期提高遥感图像的质量以达到人们所要求的某些预期结果。图像的表示 完整描述图像，可以用式子：I=f(x,y,z,t)表示一个立体的、彩色的
3、活动图像。还有：对于静止图像，则表示为f(x,y,z,)；对于平面图像，则表示为f(x,y,)；对于单色图像，则表示为f(x,y)。有时，在传播或传送图像时，常把图像扫描成一维信号，如视频信号，这时图像便成了一维函数f(t)，称之为图像信号，而前面几个式子称为图像，以示区别。数字图像数字图像是由被称作象素的小块区域组成的二维矩阵。将是由被称作象素的小块区域组成的二维矩阵。将物理图象行列划分后，每个小块区域称为物理图象行列划分后，每个小块区域称为像素像素（pixelpixel）。每个像素包括两个属性：位置和灰度。每个像素包括两个属性：位置和灰度。对于单色即对于单色即灰度图像灰度图像而言，每个象素
4、的亮度用一个数值来而言，每个象素的亮度用一个数值来表示，通常数值范围在表示，通常数值范围在0 0到到255255之间，即可用一个字节来表之间，即可用一个字节来表示，示，0 0表示黑、表示黑、255255表示白，而其它表示灰度级别。表示白，而其它表示灰度级别。物理图象及对应物理图象及对应 的数字图象的数字图象物理图像19643灰度像素数字图像采样列采样行图片像素行间隔采样列间隔灰阶黑灰白0128255彩色图象彩色图象可以用可以用红红、绿绿、蓝蓝三元组的二维矩阵来表示。三元组的二维矩阵来表示。通常，三元组的每个数值也是在通常，三元组的每个数值也是在0到到255之间，之间，0表表示相应的基色在该象素
5、中没有，而示相应的基色在该象素中没有，而255则代表相应的则代表相应的基色在该象素中取得最大值，这种情况下每个象素可基色在该象素中取得最大值，这种情况下每个象素可用三个字节来表示。用三个字节来表示。彩色图象彩色图象(128x128)及其对应的数值矩阵（仅列出一部分及其对应的数值矩阵（仅列出一部分(25x31)）(207,137,130)(220,179,163)(215,169,161)(210,179,172)(210,179,172)(207,154,146)(217,124,121)(226,144,133)(226,144,133)(224,137,124)(227,151,136)(
6、227,151,136)(226,159,142)(227,151,136)(230,170,154)(231,178,163)(231,178,163)(231,178,163)(236,187,171)(236,187,171)(239,195,176)(239,195,176)(240,205,187)(239,195,176)(231,138,123)(217,124,121)(215,169,161)(216,179,170)(216,179,170)(207,137,120)(159,51,71)(189,89,101)(216,111,110)(217,124,121)(227,
7、151,136)(227,151,136)(226,159,142)(226,159,142)(237,159,135)(237,159,135)(231,178,163)(236,187,171)(231,178,163)(236,187,171)(236,187,171)(236,187,171)(239,195,176)(239,195,176)(236,187,171)(227,133,118)(213,142,135)(216,179,170)(221,184,170)(190,89,89)(204,109,113)(204,115,118)(189,85,97)(159,60,78
8、)(136,38,65)(160,56,75)(204,109,113)(227,151,136)(226,159,142)(237,159,135)(227,151,136)主要内容 遥感图像概述 遥感图像的读写与显示 遥感图像辐射增强 遥感图像几何变换 遥感图像配准 遥感图像滤波 遥感图像分割 2 遥感图像的读写与显示 对于标准格式的图像，如bmp,jpg,tif等格式的图像，可以直接利用MATLAB提供的命令imread及imwrite进行读写操作。实际中不同的遥感数据并非标准格式。如何根据数据提供者提供的格式说明文件读出相应的参数及图像数据是进行后续遥感图像处理的关键。标准格式图像读取
9、标准格式图像读取 函数函数size可给出一副图像的行数和列数可给出一副图像的行数和列数 size（f）非标准格式遥感图像读取非标准格式遥感图像读取 通常每一种遥感卫星传感器都制定了一套适合本身特性的数据存储格式标准。并且为用户提供详细的格式说明书。用户必须熟悉格式才能正确地读出遥感图像数据。例：ALOS PALSAR参数及数据读取图像的显示 图像的写入 主要内容 遥感图像概述 遥感图像的读写与显示 遥感图像辐射增强 遥感图像几何变换 遥感图像配准 遥感图像滤波 遥感图像分割3 遥感图像辐射增强 定义：将原来不清晰的图像变得清晰或突出某些特征，同时抑制一些不需要的信息的处理方法。目的：突出图像中
10、的有用信息，扩大不同影像特征之间的差别，以便于进行判读和分析。基于直方图变换的增强方法 直方图变换是一种通过直接改变图像中像元的亮度值来改变图像的对比度，从而改善图像质量的图像处理方法。灰度直方图是灰度级的函数，描述图像中具有该灰度级的像元的个数。其横坐标是灰度级，纵坐标是像元的个数（或该灰度出现的频率）。直方图反映灰度的总体结构，但是不反映空间的分布信息。直方图变换的类型 线性变换 非线性变换线性变换：根据原图像直方图来确定好拉伸变换前的灰度值区间，然后把这一灰度值区间按某一直线方程关系拉伸或压缩而成为变换后灰度值区间。拉伸后的图像灰度值范围增大，对比度改善。线性变换 按比例扩大原是灰度级的
11、范围，将原始的相对集中的 灰度值分布在0255范围内展开。线性变换 分段线性拉伸对不同范围的灰度值进行不同的拉伸。非线性变换 变换函数是非线性的，如指数变换、对数变换等。指数变换 指数变换可以对图像高值区域进行拉伸 对数变换 对数变换对图像低值区域进行拉伸 直方图均衡 将原图像的直方图通过变换函数变为均匀的直方图，然后按均匀直方图修改原图像，从而获得一幅灰度分布均匀的新图像。直方图均衡计算步骤 1.计算原图像f的直方图h 直方图均衡计算步骤 2.求出图像f的总体像素个数N，计算每个灰度级的像素个数在整个图像中所占的比例hs。直方图均衡计算步骤 3.计算图像各灰度级的累积分布hp。直方图均衡计算
12、步骤 4.求出新图像g的灰度值。因为直方图是近似的概率密度函数，所以用离散灰度级作变换一般得不到完全平坦的结果。另外，从上例可以看出，变换后的灰度级减少了，这种现象叫做“简并”现象。由于简并现象的存在，处理后的灰度级总是要减少的，这是像素灰度有限的必然结果。由于上述原因，数字图像的直方图均衡只是近似的。主要内容 遥感图像概述 遥感图像的读写与显示 遥感图像辐射增强 遥感图像几何变换 遥感图像配准 遥感图像滤波 遥感图像分割遥感图像几何变换 空间变换主要用来保持图像中曲线的连续性和物体的连通性，一般都采用数学函数形式来描述输入、输出图像相应像素间的空间关系。空间变换的一般定义为 其中，f表示输入
13、图像，g表示输出图像，坐标(x,y)指的是空间变换后的坐标，要注意这时的坐标已经不是原来的坐标(x,y)了。a(x,y)和b(x,y)分别是图像的x和y坐标的空间变换函数。yxbyxafyxfyxg,灰度级插值主要是对空间变换后的像素赋予灰度值，使之恢复原位置处的灰度值。在几何运算中，灰度级插值是必不可少的组成部分，因为图像一般用整数位置处的像素来定义。而在几何变换中，g(x,y)的灰度值一般由处在非整数坐标上的f(x,y)的值来确定，即g中的一个像素一般对应于f中的几个像素之间的位置，反过来看也是一样，即f中的一个像素往往被映射到g中的几个像素之间的位置。图像的几何常用函数图像的几何常用函数
14、 图像的缩放图像的缩放 imresize 图像的旋转图像的旋转 imrotate 图像的剪裁图像的剪裁 imcrop 图像的一般几何变换图像的一般几何变换 imtransform 图像的插值运算图像的插值运算 最近邻插值方法（最近邻插值方法（Nearest neighbor）双线性插值方法（双线性插值方法（Bilinear）双三次插值方法（双三次插值方法（Bicubic）subplot(1,2,2),subimage(J,map)subplot(2,2,4),imshow(J1,map)X2=imcrop(X,MAP,RECT)其中，X表示有待剪切的图像，不指定X时，imcrop将当前坐标轴中
15、的图像作为待剪切的图像。MAP表示X为索引图像时的调色板，RECT定义剪切区的矩形坐标。如果调用imcrop时不指定矩形的坐标，那么当光标位于图像中时会变成十字形，可以通过拖曳鼠标的方式交互式地选择一个矩形。imcrop函数根据用户的选择绘制一个矩形，释放鼠标键后将产生一个新的图像。其中，A是变形矩阵，b是平移矩阵。变换矩阵：bAxxfSSAS00,S0 subplot(122);imshow(I1);axis on;title(尺度变换)原 图204060801020304050607080尺 度 变 换2040608010012020406080100120变换矩阵：tAt001subpl
16、ot(122);imshow(I1);axis on;title(伸缩变换)原 图204060801020304050607080伸 缩 变 换2040608020406080100120140变换矩阵：101uAusubplot(122);imshow(I1);axis on;title(扭曲变换)原 图204060801020304050607080扭 曲 变 换2040608020406080100120变换矩阵：cossinsincosrAsubplot(122);imshow(I1);title(旋转变换)原图旋转变换变换矩阵：cossinsinsincoscosSTSTUSSTST
17、USAyxDDB,subplot(122);imshow(I1);title(综合)原 图综 合主要内容 遥感图像概述 遥感图像的读写与显示 遥感图像辐射增强 遥感图像几何变换 遥感图像配准 遥感图像滤波 遥感图像分割 图像配准是指依据一些相似性度量决定图像间的变换参数使得从不同传感器不同视角不同时间获取的同一场景的两幅或多幅图像变换到统一坐标系下在像素层上得到最佳匹配的过程。待配准图像相对于参考图像的配准可定义为两幅图像在空间和亮度上的映射两幅图像可定义为两个二维数组分别用I1(x,y)和I2(x,y)表示他们分别是两幅图像的亮度值则两幅图像间的映射可表示为：其中：f为二维空间坐标变换，g为
18、一维亮度或其他亮度变换遥感图像配准步骤 第一步：读取图像 orthophoto=imread(westconcordorthophoto.png);figure,imshow(orthophoto)unregistered=imread(westconcordaerial.png);figure,imshow(unregistered)第2步：选取控制点 cpselect(unregistered,orthophoto)第3步：保存控制点到MATLAB工作空间中 File menu-choose the Export Points to Workspace 第4步(可选)：控制点微调以提高配准
19、精度 cpcorr(?)第5步：指定变换类型及其参数 mytform=cp2tform(input_points,base_points,projective);第6步：转换待配准图像 registered=imtransform(unregistered,mytform);imshow(orthophoto)hold on h=imshow(registered,gray(256);set(h,AlphaData,0.6)主要内容 遥感图像概述 遥感图像的读写与显示 遥感图像辐射增强 遥感图像几何变换 遥感图像配准 遥感图像滤波 遥感图像分割遥感图像滤波方法分类 空间域滤波 频率域滤波图像的
20、空间域滤波图像的空间域滤波滤波滤波是一种图像修正或增强技术。可以突出图像是一种图像修正或增强技术。可以突出图像的某些特征，也可以删除另一些特征。的某些特征，也可以删除另一些特征。图像滤波图像滤波的本质是一种邻域操作，输出图像的任的本质是一种邻域操作，输出图像的任一个像素值都是通过输入图像对应的像素邻域内一个像素值都是通过输入图像对应的像素邻域内的像素值利用一定的算法得到的。的像素值利用一定的算法得到的。图像的空间域滤波图像的空间域滤波就是输出图像任一个像素值都就是输出图像任一个像素值都是通过输入图像对应的像素邻域内的像素值是通过输入图像对应的像素邻域内的像素值线性线性组合组合得到的。得到的。图
21、像的空间域滤波图像的空间域滤波 二维卷积函数二维卷积函数 conv2 滤波函数滤波函数 filter2 图像滤波图像滤波 imfilter 预定义滤波器预定义滤波器 fspecial 设当前待处理像素为f(m,n)，给出一个大小为33的处理模板。以模块运算系数表示即：常用的模板：常用的模板：1111211111011H1212421211612H111101111813H0010041414141214H【例】x=imread(cameraman.tif);A=1/25*ones(5,5);x2=conv2(double(x),A);imshow(x)figureimshow(int8(x2)
22、B=0.5*0 0.25 0;0.25 1 0.25;0 0.25 0;x3=conv2(double(x),B);figureimshow(int8(x3)C=0-1 0;-1 4-1;0-1 0;x4=conv2(double(x),C);figureImshow(int8(x4)二维卷积函数二维卷积函数 conv2 A=17 24 1 8 15 23 5 7 14 16 4 6 13 20 22 10 12 19 21 3 11 18 25 2 9 h=8 1 6 3 5 7 4 9 2 C=conv2(A,h)卷积步骤：卷积步骤：1.关于卷积核的中心，旋转卷积核关于卷积核的中心，旋转卷
23、积核180度。度。2.滑动卷积核，将卷积核的中心位于图像矩阵滑动卷积核，将卷积核的中心位于图像矩阵的每一个元素。的每一个元素。3.将旋转后的卷积核作为权重，乘以对应的矩将旋转后的卷积核作为权重，乘以对应的矩阵元素阵元素4.求加权和求加权和计算卷积输出计算卷积输出(2,4)C=conv2(A,B)C=conv2(A,B，shape)A:输入图像，输入图像，B：卷积核，：卷积核，C：输出图像，若：输出图像，若A大小为大小为mana，B大小为大小为mbnb，则则C大小为（大小为（ma+mb-1）（na+nb-1）shape指定卷积运算的范围：指定卷积运算的范围：shape=full(the defa
24、ult),返回全部二维卷积结果返回全部二维卷积结果shape=same,返回与返回与A同样大小的卷积中心部分同样大小的卷积中心部分shape=valid,不考虑边界补零，返回不考虑边界补零，返回C大小为大小为（ma-mb+1）（na-nb+1）相关运算实现滤波：相关运算实现滤波：filter2相关运算相关运算(correlation)滑动相关核，将相关核的中心位于图像矩阵滑动相关核，将相关核的中心位于图像矩阵的每一个元素。的每一个元素。将相关核作为权重，乘以对应的矩阵元素将相关核作为权重，乘以对应的矩阵元素求加权和求加权和用相关运算实现滤波的函数是用相关运算实现滤波的函数是filter2B=f
25、ilter2(h,A)计算相关输出计算相关输出(2,4)图像的线性滤波图像的线性滤波 imfilter B=imfilter(A,H,option1,option2,)A:多维图像阵列，多维图像阵列，H:多维滤波器，多维滤波器，option1,option2决定边缘上的处理方法，输出图像决定边缘上的处理方法，输出图像大小，采用与大小，采用与filter2相同的方法还是卷积的方相同的方法还是卷积的方法。法。例：利用例：利用imfilter函数实现均值滤波函数实现均值滤波 I=imread(coins.png);h=ones(5,5)/25;I2=imfilter(I,h);imshow(I),t
26、itle(Original Image);figure,imshow(I2),title(Filtered Image)imfilter输入输出的数据类型是一样的输入输出的数据类型是一样的A=magic(5)h=-1 0 1;imfilter(A,h)可以看到输出有负值，所以有时候在可以看到输出有负值，所以有时候在imfilter 前用类型转换，前用类型转换，避免这种情况。避免这种情况。A=uint8(magic(5)imfilter(A,h)option:相关和卷积相关和卷积 imfilter既可以用相关，也可以用卷积实现滤波操作，既可以用相关，也可以用卷积实现滤波操作，缺省是相关。缺省是相
27、关。A=magic(5);h=-1 0 1imfilter(A,h)%filter using correlationimfilter(A,h,conv)%filter using convolutionoption:边界补零边界补零(zero-padding)和边界和边界复制复制(border replication)边界补零边界补零(zero-padding):缺省缺省option:边界补零边界补零(zero-padding)和边界和边界复制复制(border replication)边界补零边界补零(zero-padding):缺省缺省I=imread(eight.tif);h=ones
28、(5,5)/25;I2=imfilter(I,h);imshow(I),title(Original Image);figure,imshow(I2),title(Filtered Image with Black Border)缺点：滤波得到的图像边缘有一个缺点：滤波得到的图像边缘有一个dark bandoption:边界补零边界补零(zero-padding)和边界和边界复制复制(border replication)边界复制边界复制(border replication)option:边界补零边界补零(zero-padding)和边界和边界复制复制(border replication)
29、边界复制边界复制(border replication)I3=imfilter(I,h,replicate);figure,imshow(I3);title(Filtered Image with Border Replication)imfilter还有其他的边界补充选项，参考还有其他的边界补充选项，参考imfilter的帮助的帮助多维滤波多维滤波 imfilter既可以处理多维图像，也可以处理多维滤波既可以处理多维图像，也可以处理多维滤波器。用一个二维滤波器对一个三维图像滤波，相当于器。用一个二维滤波器对一个三维图像滤波，相当于对三维图像的每个平面进行二维滤波。对三维图像的每个平面进行二维
30、滤波。例例：用同样的滤波器对一个真彩色图像的每个颜色平：用同样的滤波器对一个真彩色图像的每个颜色平面进行滤波。面进行滤波。%Read in an RGB image and display itrgb=imread(peppers.png);imshow(rgb);%Filter the image and display it.h=ones(5,5)/25;rgb2=imfilter(rgb,h);figure,imshow(rgb2)imfilter与与filter2、conv2的关系：的关系：filter2、conv2、convn：将输入转换为将输入转换为double类型，输出也是类型，
31、输出也是double的，输入总是补零的，输入总是补零（zero padded),不支持其他的边界补充选项。不支持其他的边界补充选项。imfilter：不将输入转换为不将输入转换为double，输出只与输，输出只与输入同类型，有灵活的边界补充选项入同类型，有灵活的边界补充选项主要内容 空间域滤波 频率域滤波背景背景 Background 法国数学家傅立叶法国数学家傅立叶(生于生于1768年年)在在1822年年出版的出版的热分析理论热分析理论一书中指出一书中指出:任何周期任何周期函数都可以表达为不同频率的正弦和或余弦函数都可以表达为不同频率的正弦和或余弦和的形式和的形式,即傅立叶级数。即傅立叶级数
32、。20世纪世纪50年代后期年代后期,快速傅立叶变换算法出快速傅立叶变换算法出现现,得到了广泛的应用。得到了广泛的应用。傅立叶变换傅立叶变换(Fourier Transform)傅里叶变换和频率域的介绍傅里叶变换和频率域的介绍 一维傅立叶变换及其反变换一维傅立叶变换及其反变换 二维二维DFT变换及其反变换变换及其反变换 二维二维DFT变换性质变换性质一维傅立叶变换及其反变换一维傅立叶变换及其反变换 连续函数连续函数f(x)的傅立叶变换的傅立叶变换F(u):傅立叶变换傅立叶变换F(u)的反变换的反变换:dxexfuFuxj2)()(dueuFxfuxj2)()(一维傅立叶变换及其反变换一维傅立叶变
33、换及其反变换 离散函数离散函数f(x)(其中其中x，u=0,1,2,M-1)的傅立叶变换的傅立叶变换:10/2)(1)(MxMuxjexfMuF10/2)()(MxMuxjeuFxfF(u)的反变换的反变换的反变换的反变换:计算计算F(u)：1)在指数项中代入在指数项中代入u=0，然后将所有，然后将所有x 值相加值相加2)u=1，重复对所，重复对所有有x 的相加；的相加；3)对所有对所有M 个个u 重重复此过程，得到完复此过程，得到完整的整的FT。一维离散傅里叶变换一维离散傅里叶变换 离散傅里叶变换及其反变换总存在。离散傅里叶变换及其反变换总存在。用欧拉公式得用欧拉公式得sincosjej10
34、/2sin/2)cos(1)(MxMuxjMuxxfMuF每个每个F(u)由由f(x)与对应频率的正弦和余弦乘积和组成与对应频率的正弦和余弦乘积和组成；u 值决定了变换的频率成份，因此，值决定了变换的频率成份，因此，F(u)覆盖的域覆盖的域(u值值)称为频率域，其中每一项都被称为称为频率域，其中每一项都被称为FT 的频率的频率分量。与分量。与f(x)的的“时间域时间域”和和“时间成份时间成份”相对应。相对应。一维离散傅里叶变换一维离散傅里叶变换 傅里叶变换将信号分成不同频率成份。傅里叶变换将信号分成不同频率成份。类似光学中的分色棱镜把白光按波长类似光学中的分色棱镜把白光按波长(频频率率)分成不
35、同颜色，称数学棱镜。分成不同颜色，称数学棱镜。一维离散傅里叶变换一维离散傅里叶变换 傅立叶变换在极坐标下表示傅立叶变换在极坐标下表示:)()()(ujeuFuF)()()(22uIuRuF)()(arctan)(uRuIu)()()(22uIuRuP频率谱频率谱相位谱相位谱功率谱功率谱一维离散傅里叶变换一维离散傅里叶变换 f(x)是一门函数，如图所示，它表示为：是一门函数，如图所示，它表示为：)(0 0 0)(XxAxxf求其傅立叶变换求其傅立叶变换F(u)一维离散傅里叶变换一维离散傅里叶变换解：解：uXjuXjuXjuXjXuxjXuxjuxjeuXuAeeeujAeujAdxAedxexf
36、uF)sin(2 2 )()(02022一维离散傅里叶变换一维离散傅里叶变换 对应的傅立叶谱为：对应的傅立叶谱为：uXuXAXeuXuAuFuXj)sin()sin()(一维离散傅里叶变换一维离散傅里叶变换简单函数的傅里叶谱简单函数的傅里叶谱M 点离散函数及其傅里叶频点离散函数及其傅里叶频谱谱(M=1024,A=1,K=8)；对应的傅里叶频谱对应的傅里叶频谱 曲线下面积：当曲线下面积：当x 域加倍域加倍时，频率谱的高度也加倍；时，频率谱的高度也加倍；当函数长度加倍时，相同当函数长度加倍时，相同间隔下频谱中零点的数量间隔下频谱中零点的数量也加倍。也加倍。二维二维DFT傅里叶变换傅里叶变换 一个图
37、像尺寸为一个图像尺寸为MN的函数的函数f(x,y)的离散傅立叶变换的离散傅立叶变换F(u,v):F(u,v)的反变换的反变换的反变换的反变换:1010)/(2),(1),(MxNyNvyMuxjeyxfMNvuF1010)/(2),(),(MuNvNvyMuxjevuFyxf二维二维DFT傅里叶变换傅里叶变换 二维离散二维离散傅立叶变换在极坐标下表示傅立叶变换在极坐标下表示:频率谱频率谱相位谱相位谱功率谱功率谱),(),()(vujevuFuF),(),(),(22vuIvuRvuF),(),(arctan),(vuRvuIvu),(),(),(22vuIvuRvuP二维二维DFT傅里叶变换傅
38、里叶变换(u,v)=(0,0)位置的傅里叶变换值为位置的傅里叶变换值为1010),(),(1)0,0(MxNyyxfyxfMNF即即f(x,y)的均值，原点的均值，原点(0,0)的傅里叶变换是图像的的傅里叶变换是图像的平均灰度。平均灰度。F(0,0)称为频率谱的直流分量称为频率谱的直流分量(系数系数)，其它其它F(u,v)值称为交流分量值称为交流分量(交流系数交流系数)。离散傅里叶变换是对区间离散傅里叶变换是对区间0,M-1 中的中的u 值表述的，变值表述的，变换结果是关于原点对称的两个半周期，要显示完全的周换结果是关于原点对称的两个半周期，要显示完全的周期，需要将变换的原点移到期，需要将变换
39、的原点移到u=M/2，二维图像中心化亦，二维图像中心化亦是如此是如此简单二维函数的中心谱简单二维函数的中心谱 空间域和频率域抽样点之间的关系如下：空间域和频率域抽样点之间的关系如下：yNvxMu1 1简单二维函数的中心谱简单二维函数的中心谱MATLAB中离散傅立叶中离散傅立叶(DFT)的计算：的计算：fft、fft2、fftn：分别实现一维、二维和：分别实现一维、二维和n维离散快速傅立叶变换维离散快速傅立叶变换(DFT)ifft、ifft2、ifftn：分别实现一维、二维：分别实现一维、二维和和n维离散快速傅立叶逆变换维离散快速傅立叶逆变换(IDFT)例：图像傅立叶变换和幅值例：图像傅立叶变换
40、和幅值1.产生图像矩阵，包括一个矩形区域产生图像矩阵，包括一个矩形区域f=zeros(30,30);f(5:24,13:17)=1;imshow(f,InitialMagnification,fit)2.计算和显示计算和显示f的傅立叶变换的傅立叶变换F=fft2(f);F2=log(abs(F);figure,imshow(F2,-1 5,InitialMagnification,fit);colormap(jet);colorbar -1012345上图：傅立叶变换的取样很粗，且零频率系数不像上图：傅立叶变换的取样很粗，且零频率系数不像传统的那样显示在图像中心，而是显示在左上角。传统的那样显
41、示在图像中心，而是显示在左上角。为获得傅立叶变换较精细的取样，采用补零的方法。为获得傅立叶变换较精细的取样，采用补零的方法。3.补零和计算补零和计算f的傅立叶变换的傅立叶变换F=fft2(f,256,256);F2=log(abs(F);figure,imshow(F2,-1 5,InitialMagnification,fit);colormap(jet);colorbar4.将零频系数移动到图像中心将零频系数移动到图像中心F=fft2(f,256,256);F2=fftshift(F);figure,imshow(log(abs(F2),-1 5,InitialMagnification,
42、fit);colormap(jet);colorbar频率域滤波频率域滤波 频率域的基本性质频率域的基本性质 每个每个F(u,v)项包含了被指数项修正的项包含了被指数项修正的f(x,y)的所有值：的所有值：1010)/(2),(1),(MxNyNvyMuxjeyxfMNvuF直观上将傅里叶变换和图像中的亮度变化联系起来并直观上将傅里叶变换和图像中的亮度变化联系起来并不困难：不困难：直流分量直流分量F(0,0)对应一幅图像的平均灰度；对应一幅图像的平均灰度；低频部分对应图像缓慢变化的分量；低频部分对应图像缓慢变化的分量；高频部分对应图像边缘和灰度级突变的部分高频部分对应图像边缘和灰度级突变的部分
43、图为一幅集成电路的扫描电图为一幅集成电路的扫描电子显微镜子显微镜(SEM)图像，放大图像，放大将近将近2500倍。注意图中倍。注意图中45的强边缘，和两个的强边缘，和两个因热感应不足而产生的白色因热感应不足而产生的白色氧化突起氧化突起图是上图的傅里叶频谱，沿图是上图的傅里叶频谱，沿着着45方向对应上图边缘方向对应上图边缘突起部分。沿垂直轴偏左部突起部分。沿垂直轴偏左部分有垂直分量，由氧化突起分有垂直分量，由氧化突起的上下黑白边沿形成。的上下黑白边沿形成。频率域的基本性质频率域的基本性质 频率域的基本性质：频率域的基本性质：频域的中心邻域频域的中心邻域对应图像中慢变化部分，较高的频率开对应图像中
44、慢变化部分，较高的频率开始对应图像中变化较快的部分始对应图像中变化较快的部分（如：物（如：物体的边缘、线条等）。体的边缘、线条等）。频率域中滤波步骤频率域中滤波步骤输入图像前处理傅里叶变换滤波函数傅里叶反变换后处理增强后的图像前处理、后处理：前处理、后处理：1.中心变换中心变换2.输入图像向其最接近的偶数维转换输入图像向其最接近的偶数维转换3.灰度级标定灰度级标定4.输入向浮点的转换输入向浮点的转换5.输出向输出向8比特整数的转换比特整数的转换 1.用用(-1)x+y乘以输入图像来进行中心变化。在乘以输入图像来进行中心变化。在MATLAB中是利用函数中是利用函数fftshift()实现。实现。
45、2.由由(1)计算图像的计算图像的DFT，即，即F(u,v)；3.用滤波函数用滤波函数H(u,v)乘以乘以F(u,v)。H(u,v)称为滤称为滤波器：抑制某些波器：抑制某些频率，其他频率频率，其他频率不受影响不受影响 频率域中滤波步骤：频率域中滤波步骤：4.计算计算(3)中结果的反中结果的反DFT。5.得到得到(4)中结果的实部。中结果的实部。6.用用(-1)x+y乘以乘以(5)中的结果。中的结果。频率域中滤波步骤：频率域中滤波步骤：f=imread(cameraman.tif);PQ=paddedsize(size(f);%zero paddedF=fft2(f,PQ(1),PQ(2);h=
46、fspecial(laplacian,0.8)H=freqz2(h,PQ(1),PQ(2);H1=ifftshift(H);G=H1.*F;g=real(ifft2(G);gf=g(1:size(f,1),1:size(f,2);imshow(f)figureimshow(gf)一些基本的滤波器及其性质一些基本的滤波器及其性质l低通滤波器：使低频通过，高频衰减低通滤波器：使低频通过，高频衰减低频主要决定图像在平滑区域中总体灰度级的显示低频主要决定图像在平滑区域中总体灰度级的显示比原始图像少一些尖锐的细节部分比原始图像少一些尖锐的细节部分l高通滤波器：使高频通过，低频衰减高通滤波器：使高频通过，
47、低频衰减高频决定图像细节部分，如边缘和噪声高频决定图像细节部分，如边缘和噪声在平滑区域中减少灰度级变化，突出过渡（如边缘）在平滑区域中减少灰度级变化，突出过渡（如边缘）灰度级的细节部分，使图像更加锐化。灰度级的细节部分，使图像更加锐化。基本的滤波器及其性质基本的滤波器及其性质图像被模糊图像被模糊锐化锐化F(0,0)=0:几乎几乎没有平滑细节没有平滑细节周期对称Lowpass filterHighpass filter 空间域滤波和频率域滤波之间的对应关系空间域滤波和频率域滤波之间的对应关系 卷积定理：卷积定理：空间域的乘法对应频域卷积空间域的乘法对应频域卷积空间域滤波和频率域滤波之间的对应关系
48、空间域滤波和频率域滤波之间的对应关系 重要性质:),(),(),(),(),(),(),(),(),(),(yxHyxhyxHyxyxhyxvuHvuFyxhyxf根据冲击函数和卷积定理的性质，可知空间域和频率根据冲击函数和卷积定理的性质，可知空间域和频率域的滤波器组成傅里叶变换对域的滤波器组成傅里叶变换对h(x,y)和和H(u,v)。给出频。给出频率域滤波器率域滤波器H(u,v)，通过反傅里叶变换可以得到空间域，通过反傅里叶变换可以得到空间域相应的滤波器相应的滤波器h(x,y)。空间域滤波和频率域滤波之间的对应关系空间域滤波和频率域滤波之间的对应关系 滤波器大小滤波器大小前述的所有函数均具有
49、相同的尺寸前述的所有函数均具有相同的尺寸MN。在。在实际中，指定一个频率域滤波器，进行反变实际中，指定一个频率域滤波器，进行反变换后会得到一个相同尺寸的空间域滤波器。换后会得到一个相同尺寸的空间域滤波器。如果两个域中滤波器尺寸相同，那么通常频如果两个域中滤波器尺寸相同，那么通常频域中进行滤波计算更为有效，更为直观，但域中进行滤波计算更为有效，更为直观，但空域中适用更小尺寸的滤波器，更为有效。空域中适用更小尺寸的滤波器，更为有效。空间域滤波和频率域滤波之间的对应关系空间域滤波和频率域滤波之间的对应关系 方法：方法：在频率域指定滤波器；在频率域指定滤波器；做反变换；做反变换；使用结果滤波器作为在空
50、间域构建更小空使用结果滤波器作为在空间域构建更小空间滤波模板的指导；间滤波模板的指导；基于高斯函数的滤波基于高斯函数的滤波 高斯滤波器函数高斯滤波器函数(低通低通):222/)(uAeuH对应的空间域滤波器对应的空间域滤波器22222)(xAexh 为高斯曲线的标准差。为高斯曲线的标准差。基于高斯函数的滤波基于高斯函数的滤波高斯滤波器的重要特性高斯滤波器的重要特性频域和空域高斯滤波器构成傅里叶变换对，且频域和空域高斯滤波器构成傅里叶变换对，且都是实高斯函数。处理时不用考虑复数，而且高都是实高斯函数。处理时不用考虑复数，而且高斯曲线直观，易于操作。斯曲线直观，易于操作。高斯滤波器傅里叶变换对之间
