MATLAB在遥感图像处理中的应用

整理课件1第五章MATLAB在遥感图像处理中的应用

整理课件2主要内容遥感图像概述遥感图像的读写与显示遥感图像辐射增强遥感图像几何变换遥感图像配准遥感图像滤波遥感图像分割整理课件31遥感图像处理概述基本概念遥感数字图像：是指被计算机存储、处理和使用的图像，是一种空间坐标和灰度均不连续的、用数字形式表达的遥感影像，最基本单元是像素，每个像素具有空间位置特征和属性特征。空间位置特征：是用离散的X值和Y值来表示；属性特征：常用亮度值表示。整理课件4亮度值有如下特点：（1）不同图像相同地点的亮度值不同；（2）亮度值大小由传感器所探测到的电磁辐射强度决定，入射到传感器中的电磁波被探测元件转化为电信号，经过A/D转换，成为绝对辐射亮度值R。为了便于应用R又被转换为能够表征地物的辐射亮度相对值V。

R=V*（Rmax-Rmin）/Dmax+RminRmax为探测器可检测到的最大辐射亮度；Rmin为探测器可检测到的最小辐射亮度；max为级数；R为辐射亮度值；V为像素表征的地物辐射亮度的相对值。整理课件5遥感数字图像处理：利用数字计算机或其它高速、大规模集成数字硬件，对从遥感图像信息转换来的数字电信号进行某些数字运算或处理（如去除噪声、增强、复原、分割、提取特征等），以期提高遥感图像的质量以达到人们所要求的某些预期结果。整理课件6图像的表示

完整描述图像，可以用式子：

I=f(x,y,z,λ,t)

表示一个立体的、彩色的活动图像。还有：对于静止图像，则表示为f(x,y,z,λ)；对于平面图像，则表示为f(x,y,λ)；对于单色图像，则表示为f(x,y)。有时，在传播或传送图像时，常把图像扫描成一维信号，如视频信号，这时图像便成了一维函数f(t)，称之为图像信号，而前面几个式子称为图像，以示区别。整理课件7数字图像是由被称作象素的小块区域组成的二维矩阵。将物理图象行列划分后，每个小块区域称为像素（pixel）。

–每个像素包括两个属性：位置和灰度。对于单色即灰度图像而言，每个象素的亮度用一个数值来表示，通常数值范围在0到255之间，即可用一个字节来表示，

0表示黑、255表示白，而其它表示灰度级别。物理图象及对应的数字图象整理课件8物理图像19643灰度像素数字图像采样列采样行图片像素行间隔采样列间隔灰阶黑灰白0128255整理课件9彩色图象可以用红、绿、蓝三元组的二维矩阵来表示。

–通常，三元组的每个数值也是在0到255之间，0表示相应的基色在该象素中没有，而255则代表相应的基色在该象素中取得最大值，这种情况下每个象素可用三个字节来表示。彩色图象(128x128)及其对应的数值矩阵（仅列出一部分(25x31)）整理课件10主要内容遥感图像概述遥感图像的读写与显示遥感图像辐射增强遥感图像几何变换遥感图像配准遥感图像滤波遥感图像分割整理课件112遥感图像的读写与显示对于标准格式的图像，如bmp,jpg,tif等格式的图像，可以直接利用MATLAB提供的命令imread及imwrite进行读写操作。实际中不同的遥感数据并非标准格式。如何根据数据提供者提供的格式说明文件读出相应的参数及图像数据是进行后续遥感图像处理的关键。整理课件12标准格式图像读取函数imread可以从任何Matlab支持的图像文件格式中，以任意位深度读取一幅图像。格式为：[X,MAP]=imread(FILENAME,'FMT')，其中：FILENAME－为需要读入的图像文件名称，FMT－为图像格式。图像文件格式整理课件13读取图像信息可以通过调用imfinfo函数获得与图像文件有关的信息，格式如下：INFO=imfinfo(FILENAME,'FMT')其中：返回的INFO是Matlab的一个结构体。函数size可给出一副图像的行数和列数size（f）整理课件14非标准格式遥感图像读取通常每一种遥感卫星传感器都制定了一套适合本身特性的数据存储格式标准。并且为用户提供详细的格式说明书。用户必须熟悉格式才能正确地读出遥感图像数据。例：ALOSPALSAR参数及数据读取整理课件15图像的显示Matlab的图像显示函数主要有imshow、colorbar、subimage和imagesc等。例如：imshow函数的基本语法为：

imshow(f,G)

其中，f是一个图像数组，G是显示该图像的灰度级数。整理课件161、imshow(I,n)显示灰度级为n的图像，n缺省为256。【例】按256灰度级显示I=imread('moon.tif');imshow(I,256)colorbar整理课件17【例】按8灰度级显示I=imread('moon.tif');imshow(I,8)colorbar整理课件18【例】按最大灰度范围显示I=imread('moon.tif');imshow(I,[])colorbar整理课件19【例】按指定灰度范围显示I=imread('moon.tif');imshow(I,[64,128])colorbar整理课件202、imshow(X,MAP)显示索引图像，X为数据图像矩阵，MAP为调色板。【例】显示索引图像[X,MAP]=imread('canoe.tif');imshow(X,MAP);colorbar整理课件21【例】显示伪彩色图像I=imread('moon.tif');imshow(I,jet(256))colorbar整理课件223、imshow(RGB)显示真彩色图像。【例】rgb=imread('flowers.tif');imshow(rgb);整理课件234、subimage多图显示多个调色板的图像。【例】loadtrees;[x2,map2]=imread('forest.tif');subplot(2,1,1),subimage(X,map);colorbarsubplot(2,1,2),subimage(x2,map2);colorbar

整理课件24【例】不用专用函数显示多图loadtrees;[x2,map2]=imread('forest.tif');subplot(2,1,1),imshow(x2,map2);colorbarsubplot(2,1,2),imshow(X,map);colorbar整理课件25图像的写入1、imwrite函数imwrite(I,'文件名','文件格式')－保存无调色板的图像imwrite(X,MAP,'文件名','文件格式')－保存有调色板的图像imwrite(...,'文件名')－文件名中含格式imwrite(...,'参数','值')－指定保存参数文件格式：bmp、jpg、pcx、tif。整理课件26【例】将tif图像保存为jpg图像[x,map]=imread('canoe.tif');imwrite(x,map,'canoe.jpg','JPG','Quality',75)整理课件27主要内容遥感图像概述遥感图像的读写与显示遥感图像辐射增强遥感图像几何变换遥感图像配准遥感图像滤波遥感图像分割整理课件283遥感图像辐射增强定义：将原来不清晰的图像变得清晰或突出某些特征，同时抑制一些不需要的信息的处理方法。目的：突出图像中的有用信息，扩大不同影像特征之间的差别，以便于进行判读和分析。整理课件29基于直方图变换的增强方法直方图变换是一种通过直接改变图像中像元的亮度值来改变图像的对比度，从而改善图像质量的图像处理方法。灰度直方图是灰度级的函数，描述图像中具有该灰度级的像元的个数。其横坐标是灰度级，纵坐标是像元的个数（或该灰度出现的频率）。直方图反映灰度的总体结构，但是不反映空间的分布信息。直方图变换的类型线性变换非线性变换整理课件30线性变换：根据原图像直方图来确定好拉伸变换前的灰度值区间，然后把这一灰度值区间按某一直线方程关系拉伸或压缩而成为变换后灰度值区间。拉伸后的图像灰度值范围增大，对比度改善。整理课件31线性变换按比例扩大原是灰度级的范围，将原始的相对集中的灰度值分布在0–255范围内展开。整理课件32线性变换分段线性拉伸对不同范围的灰度值进行不同的拉伸。整理课件33非线性变换变换函数是非线性的，如指数变换、对数变换等。指数变换指数变换可以对图像高值区域进行拉伸整理课件34对数变换对数变换对图像低值区域进行拉伸整理课件35直方图均衡将原图像的直方图通过变换函数变为均匀的直方图，然后按均匀直方图修改原图像，从而获得一幅灰度分布均匀的新图像。整理课件36直方图均衡计算步骤1.计算原图像f的直方图h整理课件37直方图均衡计算步骤2.求出图像f的总体像素个数N，计算每个灰度级的像素个数在整个图像中所占的比例hs。整理课件38直方图均衡计算步骤3.计算图像各灰度级的累积分布hp。整理课件39直方图均衡计算步骤4.求出新图像g的灰度值。整理课件40整理课件41因为直方图是近似的概率密度函数，所以用离散灰度级作变换一般得不到完全平坦的结果。另外，从上例可以看出，变换后的灰度级减少了，这种现象叫做“简并”现象。由于简并现象的存在，处理后的灰度级总是要减少的，这是像素灰度有限的必然结果。由于上述原因，数字图像的直方图均衡只是近似的。整理课件42主要内容遥感图像概述遥感图像的读写与显示遥感图像辐射增强遥感图像几何变换遥感图像配准遥感图像滤波遥感图像分割整理课件43遥感图像几何变换几何运算可以看成是像素在图像内的移动过程，该移动过程可以改变图像中物体对象(像素)之间的空间关系。几何运算可以是不受任何限制的，但是通常都需要做出一些限制以保持图像的外观顺序。完整的几何运算需要由两个算法来实现：空间变换算法和灰度插值算法。整理课件44空间变换主要用来保持图像中曲线的连续性和物体的连通性，一般都采用数学函数形式来描述输入、输出图像相应像素间的空间关系。空间变换的一般定义为其中，f表示输入图像，g表示输出图像，坐标(x',y')指的是空间变换后的坐标，要注意这时的坐标已经不是原来的坐标(x,y)了。a(x,y)和b(x,y)分别是图像的x和y坐标的空间变换函数。整理课件45灰度级插值主要是对空间变换后的像素赋予灰度值，使之恢复原位置处的灰度值。在几何运算中，灰度级插值是必不可少的组成部分，因为图像一般用整数位置处的像素来定义。而在几何变换中，g(x,y)的灰度值一般由处在非整数坐标上的f(x,y)的值来确定，即g中的一个像素一般对应于f中的几个像素之间的位置，反过来看也是一样，即f中的一个像素往往被映射到g中的几个像素之间的位置。整理课件46图像的几何常用函数图像的缩放imresize图像的旋转imrotate图像的剪裁imcrop图像的一般几何变换imtransform

整理课件47图像的插值运算最近邻插值方法（Nearestneighbor）双线性插值方法（Bilinear）双三次插值方法（Bicubic）整理课件48(1)图像缩放MATLAB使用imresize函数来改变一幅图像的大小,调用格式如下：B=imresize(A,M,METHOD)，其中：A－原图像；M－缩放系数；B－缩放后的图像；METHOD－插值方法，可取值‘nearest’，‘bilinear’和‘bicubic’。【例】原图像放大1.25倍[I,map]=imread('kids.tif');J=imresize(I,3);subplot(1,2,1),subimage(I,map)subplot(1,2,2),subimage(J,map)整理课件49(2)图像旋转使用imrotate函数来旋转一幅图像。调用格式如下：B=imrotate(A,ANGLE,METHOD,BBOX)，其中：A－需要旋转的图像；ANGLE－表示旋转的角度，正值为逆时针；METHOD－插值方法；BBOX－取值loose（缺省）、crop【例】图像旋转clf[I,map]=imread('kids.tif');J=imrotate(I,35,'bilinear');J1=imrotate(I,35,'bilinear','crop');subplot(2,2,1),imshow(I,map)subplot(2,2,3),imshow(J,map)subplot(2,2,4),imshow(J1,map)整理课件50(3)图像剪切使用imcrop函数可以从一幅图像中抽取一个矩形的部分。imcrop函数的调用格式如下：X2=imcrop(X,MAP,RECT)其中，X表示有待剪切的图像，不指定X时，imcrop将当前坐标轴中的图像作为待剪切的图像。MAP表示X为索引图像时的调色板，RECT定义剪切区的矩形坐标。如果调用imcrop时不指定矩形的坐标，那么当光标位于图像中时会变成十字形，可以通过拖曳鼠标的方式交互式地选择一个矩形。imcrop函数根据用户的选择绘制一个矩形，释放鼠标键后将产生一个新的图像。整理课件51仿射变换仿射变换，可以用以下函数来描述：其中，A是变形矩阵，b是平移矩阵。

(1)尺度变换变换矩阵：,S≥0整理课件52【例】clf;I=checkerboard(20,2);subplot(121);imshow(I);axison;title('原图')s=1.5;T=[s0;0s;00];tf=maketform('affine',T);I1=imtransform(I,tf,'bicubic','FillValues',0.3);subplot(122);imshow(I1);axison;title('尺度变换')整理课件53(2)伸缩变换变换矩阵：整理课件54【例】clf;I=checkerboard(20,2);subplot(121);imshow(I);axison;title('原图')t=2;T=[10;0t;00];tf=maketform('affine',T);I1=imtransform(I,tf,'bicubic','FillValues',0.3);subplot(122);imshow(I1);axison;title('伸缩变换')整理课件55(3)扭曲变换变换矩阵：整理课件56【例】clf;I=checkerboard(20,2);subplot(121);imshow(I);axison;title('原图')u=0.5;T=[1u;01;00];tf=maketform('affine',T);I1=imtransform(I,tf,'bicubic','FillValues',0.3);subplot(122);imshow(I1);axison;title('扭曲变换')整理课件57(4)旋转变换变换矩阵：整理课件58clf;I=checkerboard(20,2);subplot(1,2,1);imshow(I);title('原图')angle=20*pi/180;sc=cos(angle);ss=sin(angle);T=[sc-ss;sssc;00];tf=maketform('affine',T);I1=imtransform(I,tf,'bicubic','FillValues',0.3);subplot(122);imshow(I1);title('旋转变换')整理课件59(5)综合变换变换矩阵：整理课件60【例】clf;I=checkerboard(20,2);subplot(1,2,1);imshow(I);title('原图')s=2;As=[s0;0s];%尺度t=2;At=[10;0t];%伸缩u=1.5;Au=[1u;01];%扭曲st=30*pi/180;sc=cos(st);ss=sin(st);Ast=[sc-ss;sssc];%旋转T=[As*At*Au*Ast;35];tf=maketform('affine',T);I1=imtransform(I,tf,'bicubic','FillValues',0.3);subplot(122);imshow(I1);title('综合')整理课件61主要内容遥感图像概述遥感图像的读写与显示遥感图像辐射增强遥感图像几何变换遥感图像配准遥感图像滤波遥感图像分割整理课件62图像配准是指依据一些相似性度量决定图像间的变换参数使得从不同传感器不同视角不同时间获取的同一场景的两幅或多幅图像变换到统一坐标系下在像素层上得到最佳匹配的过程。待配准图像相对于参考图像的配准可定义为两幅图像在空间和亮度上的映射两幅图像可定义为两个二维数组分别用I1(x,y)和I2(x,y)表示他们分别是两幅图像的亮度值则两幅图像间的映射可表示为：其中：f为二维空间坐标变换，g为一维亮度或其他亮度变换整理课件63整理课件64遥感图像配准步骤第一步：读取图像orthophoto=imread('westconcordorthophoto.png');figure,imshow(orthophoto)unregistered=imread('westconcordaerial.png');figure,imshow(unregistered)整理课件65第2步：选取控制点cpselect(unregistered,orthophoto)整理课件66第3步：保存控制点到MATLAB工作空间中Filemenu->choosetheExportPointstoWorkspace第4步(可选)：控制点微调以提高配准精度cpcorr(??)整理课件67第5步：指定变换类型及其参数mytform=cp2tform(input_points,base_points,'projective');第6步：转换待配准图像registered=imtransform(unregistered,mytform);imshow(orthophoto)holdonh=imshow(registered,gray(256));set(h,'AlphaData',0.6)整理课件68主要内容遥感图像概述遥感图像的读写与显示遥感图像辐射增强遥感图像几何变换遥感图像配准遥感图像滤波

遥感图像分割整理课件69遥感图像滤波方法分类空间域滤波频率域滤波整理课件70图像的空间域滤波滤波是一种图像修正或增强技术。可以突出图像的某些特征，也可以删除另一些特征。图像滤波的本质是一种邻域操作，输出图像的任一个像素值都是通过输入图像对应的像素邻域内的像素值利用一定的算法得到的。图像的空间域滤波就是输出图像任一个像素值都是通过输入图像对应的像素邻域内的像素值线性组合得到的。整理课件71图像的空间域滤波二维卷积函数conv2滤波函数filter2图像滤波imfilter预定义滤波器fspecial

整理课件72设当前待处理像素为f(m,n)，给出一个大小为3×3的处理模板。整理课件73以模块运算系数表示即：整理课件74常用的模板：整理课件75【例】x=imread(‘cameraman.tif’);A=1/25*ones(5,5);x2=conv2(double(x),A);imshow(x)figureimshow(int8(x2))B=0.5*[00.250;0.2510.25;00.250];x3=conv2(double(x),B);figureimshow(int8(x3))C=[0-10;-14-1;0-10];x4=conv2(double(x),C);figureImshow(int8(x4))整理课件76二维卷积函数conv2

A=[

17241815235714164613202210121921311182529]h=[816357492]

C=conv2(A,h)整理课件77卷积步骤：关于卷积核的中心，旋转卷积核180度。滑动卷积核，将卷积核的中心位于图像矩阵的每一个元素。将旋转后的卷积核作为权重，乘以对应的矩阵元素求加权和整理课件78计算卷积输出(2,4)整理课件79C=conv2(A,B)C=conv2(A,B，shape)A:输入图像，B：卷积核，C：输出图像，若A大小为ma×na，B大小为mb×nb，则C大小为（ma+mb-1）×（na+nb-1）shape指定卷积运算的范围：shape=‘full’(thedefault),返回全部二维卷积结果shape=‘same’,返回与A同样大小的卷积中心部分shape=‘valid’,不考虑边界补零，返回C大小为（ma-mb+1）×（na-nb+1）整理课件80相关运算实现滤波：filter2相关运算(correlation)滑动相关核，将相关核的中心位于图像矩阵的每一个元素。将相关核作为权重，乘以对应的矩阵元素求加权和用相关运算实现滤波的函数是filter2B=filter2(h,A)整理课件81计算相关输出(2,4)整理课件82图像的线性滤波imfilterB=imfilter(A,H,option1,option2,…)A:多维图像阵列，H:多维滤波器，option1,option2…决定边缘上的处理方法，输出图像大小，采用与filter2相同的方法还是卷积的方法。整理课件83例：利用imfilter函数实现均值滤波

I=imread('coins.png');h=ones(5,5)/25;I2=imfilter(I,h);imshow(I),title('OriginalImage');figure,imshow(I2),title('FilteredImage')整理课件84imfilter输入输出的数据类型是一样的A=magic(5)h=[-101];imfilter(A,h)可以看到输出有负值，所以有时候在imfilter前用类型转换，避免这种情况。A=uint8(magic(5))imfilter(A,h)整理课件85option:相关和卷积imfilter既可以用相关，也可以用卷积实现滤波操作，缺省是相关。A=magic(5);h=[-101]imfilter(A,h)%filterusingcorrelationimfilter(A,h,'conv')

%filterusingconvolution整理课件86option:边界补零(zero-padding)和边界复制(borderreplication)边界补零(zero-padding):缺省整理课件87option:边界补零(zero-padding)和边界复制(borderreplication)边界补零(zero-padding):缺省I=imread('eight.tif');h=ones(5,5)/25;I2=imfilter(I,h);imshow(I),title('OriginalImage');figure,imshow(I2),title('FilteredImagewithBlackBorder')缺点：滤波得到的图像边缘有一个darkband整理课件88option:边界补零(zero-padding)和边界复制(borderreplication)边界复制(borderreplication)整理课件89option:边界补零(zero-padding)和边界复制(borderreplication)边界复制(borderreplication)

I3=imfilter(I,h,'replicate');figure,imshow(I3);title('FilteredImagewithBorderReplication')imfilter还有其他的边界补充选项，参考imfilter的帮助整理课件90多维滤波imfilter既可以处理多维图像，也可以处理多维滤波器。用一个二维滤波器对一个三维图像滤波，相当于对三维图像的每个平面进行二维滤波。例：用同样的滤波器对一个真彩色图像的每个颜色平面进行滤波。整理课件91%ReadinanRGBimageanddisplayitrgb=imread('peppers.png');imshow(rgb);%Filtertheimageanddisplayit.h=ones(5,5)/25;rgb2=imfilter(rgb,h);figure,imshow(rgb2)整理课件92imfilter与filter2、conv2的关系：filter2、conv2、convn：将输入转换为double类型，输出也是double的，输入总是补零（zeropadded),不支持其他的边界补充选项。imfilter：不将输入转换为double，输出只与输入同类型，有灵活的边界补充选项整理课件93主要内容空间域滤波频率域滤波整理课件94背景

Background法国数学家傅立叶(生于1768年)在1822年出版的《热分析理论》一书中指出:任何周期函数都可以表达为不同频率的正弦和或余弦和的形式,即傅立叶级数。20世纪50年代后期,快速傅立叶变换算法出现,得到了广泛的应用。傅立叶变换(FourierTransform)整理课件95傅里叶变换和频率域的介绍一维傅立叶变换及其反变换二维DFT变换及其反变换二维DFT变换性质整理课件96一维傅立叶变换及其反变换连续函数f(x)的傅立叶变换F(u):傅立叶变换F(u)的反变换:整理课件97一维傅立叶变换及其反变换离散函数f(x)(其中x，u=0,1,2,…,M-1)的傅立叶变换:F(u)的反变换的反变换:计算F(u)：1)在指数项中代入u=0，然后将所有x值相加2)u=1，重复对所有x的相加；3)对所有M个u重复此过程，得到完整的FT。整理课件98一维离散傅里叶变换离散傅里叶变换及其反变换总存在。用欧拉公式得每个F(u)

由f(x)与对应频率的正弦和余弦乘积和组成；u值决定了变换的频率成份，因此，F(u)覆盖的域(u值)称为频率域，其中每一项都被称为FT的频率分量。与f(x)的“时间域”和“时间成份”相对应。整理课件99一维离散傅里叶变换傅里叶变换将信号分成不同频率成份。类似光学中的分色棱镜把白光按波长(频率)分成不同颜色，称数学棱镜。整理课件100一维离散傅里叶变换傅立叶变换在极坐标下表示:频率谱相位谱功率谱整理课件101一维离散傅里叶变换f(x)是一门函数，如图所示，它表示为：求其傅立叶变换F(u)整理课件102一维离散傅里叶变换解：整理课件103一维离散傅里叶变换对应的傅立叶谱为：整理课件104一维离散傅里叶变换简单函数的傅里叶谱M点离散函数及其傅里叶频谱(M=1024,A=1,K=8)；对应的傅里叶频谱

曲线下面积：当x域加倍时，频率谱的高度也加倍；当函数长度加倍时，相同间隔下频谱中零点的数量也加倍。整理课件105二维DFT傅里叶变换一个图像尺寸为M×N的函数f(x,y)的离散傅立叶变换F(u,v):F(u,v)的反变换的反变换:整理课件106二维DFT傅里叶变换二维离散傅立叶变换在极坐标下表示:频率谱相位谱功率谱整理课件107二维DFT傅里叶变换(u,v)=(0,0)位置的傅里叶变换值为即f(x,y)的均值，原点(0,0)的傅里叶变换是图像的平均灰度。F(0,0)称为频率谱的直流分量(系数)，其它F(u,v)值称为交流分量(交流系数)。整理课件108离散傅里叶变换是对区间[0,M-1]中的u值表述的，变换结果是关于原点对称的两个半周期，要显示完全的周期，需要将变换的原点移到u=M/2，二维图像中心化亦是如此整理课件109简单二维函数的中心谱空间域和频率域抽样点之间的关系如下：整理课件110简单二维函数的中心谱整理课件111MATLAB中离散傅立叶(DFT)的计算：fft、fft2、fftn：分别实现一维、二维和n维离散快速傅立叶变换(DFT)ifft、ifft2、ifftn：分别实现一维、二维和n维离散快速傅立叶逆变换(IDFT)整理课件112例：图像傅立叶变换和幅值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,[-15],'InitialMagnification','fit');colormap(jet);colorbar整理课件113整理课件114上图：傅立叶变换的取样很粗，且零频率系数不像传统的那样显示在图像中心，而是显示在左上角。为获得傅立叶变换较精细的取样，采用补零的方法。3.补零和计算f的傅立叶变换F=fft2(f,256,256);F2=log(abs(F));figure,imshow(F2,[-15],'InitialMagnification','fit');colormap(jet);colorbar整理课件1154.将零频系数移动到图像中心F=fft2(f,256,256);F2=fftshift(F);figure,imshow(log(abs(F2)),[-15],'InitialMagnification','fit');colormap(jet);colorbar整理课件116频率域滤波频率域的基本性质

每个F(u,v)项包含了被指数项修正的f(x,y)的所有值：直观上将傅里叶变换和图像中的亮度变化联系起来并不困难：直流分量F(0,0)对应一幅图像的平均灰度；低频部分对应图像缓慢变化的分量；高频部分对应图像边缘和灰度级突变的部分整理课件117图为一幅集成电路的扫描电子显微镜(SEM)图像，放大将近2500倍。注意图中±45°的强边缘，和两个因热感应不足而产生的白色氧化突起图是上图的傅里叶频谱，沿着±45°方向对应上图边缘突起部分。沿垂直轴偏左部分有垂直分量，由氧化突起的上下黑白边沿形成。整理课件118频率域的基本性质

频率域的基本性质：频域的中心邻域对应图像中慢变化部分，较高的频率开始对应图像中变化较快的部分（如：物体的边缘、线条等）。整理课件119频率域中滤波步骤输入图像前处理傅里叶变换滤波函数傅里叶反变换后处理增强后的图像前处理、后处理：1.中心变换2.输入图像向其最接近的偶数维转换3.灰度级标定4.输入向浮点的转换5.输出向8比特整数的转换整理课件1201.用(-1)x+y乘以输入图像来进行中心变化。在MATLAB中是利用函数fftshift()实现。2.由(1)计算图像的DFT，即F(u,v)；3.用滤波函数H(u,v)乘以F(u,v)。

H(u,v)称为滤波器：抑制某些频率，其他频率不受影响

频率域中滤波步骤：整理课件1214.计算(3)中结果的反DFT。5.得到(4)中结果的实部。6.用(-1)x+y乘以(5)中的结果。

频率域中滤波步骤：整理课件122f=imread('cameraman.tif');PQ=paddedsize(size(f));%zeropaddedF=fft2(f,PQ(1),PQ(2));h=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)整理课件123一些基本的滤波器及其性质低通滤波器：使低频通过，高频衰减低频主要决定图像在平滑区域中总体灰度级的显示比原始图像少一些尖锐的细节部分高通滤波器：使高频通过，低频衰减高频决定图像细节部分，如边缘和噪声在平滑区域中减少灰度级变化，突出过渡（如边缘）灰度级的细节部分，使图像更加锐化。整理课件124基本的滤波器及其性质图像被模糊锐化F(0,0)=0:几乎没有平滑细节周期对称LowpassfilterHighpassfilter

整理课件125空间域滤波和频率域滤波之间的对应关系卷积定理：

空间域的乘法对应频域卷积整理课件126空间域滤波和频率域滤波之间的对应关系重要性质:根据冲击函数和卷积定理的性质，可知空间域和频率域的滤波器组成傅里叶变换对h(x,y)和H(u,v)。给出频率域滤波器H(u,v)，通过反傅里叶变换可以得到空间域相应的滤波器h(x,y)。整理课件127空间域滤波和频率域滤波之间的对应关系滤波器大小前述的所有函数均具有相同的尺寸M×N。在实际中，指定一个频率域滤波器，进行反变换后会得到一个相同尺寸的空间域滤波器。如果两个域中滤波器尺寸相同，那么通常频域中进行滤波计算更为有效，更为直观，但空域中适用更小尺寸的滤波器，更为有效。整理课件128空间域滤波和频率域滤波之间的对应关系方法：在频率域指定滤波器；做反变换；使用结果滤波器作为在空间域构建更小空间滤波模板的指导；整理课件129基于高斯函数的滤波高斯滤波器函数(低通):对应的空间域滤波器δ为高斯曲线的标准差。整理课件130基于高斯函数的滤波高斯滤波器的重要特性频域和空域高斯滤波器构成傅里叶变换对，且都是实高斯函数。处理时不用考虑复数，而且高斯曲线直观，易于操作。高斯滤波器傅里叶变换对之间有相互作用：当H(u)有很宽轮廓时(大δ值)，h(x)很窄轮廓，反之亦然；当δ趋于无穷时，H(u)趋于常函数，h(x)趋于冲击函数。整理课件131基于高斯函数的滤波器频域滤波器越窄，滤除的低频部分越多，图像越模糊。在空域中意味着滤波器越宽，模板就越大(阶数越大)。整理课件132空间域滤波和频域滤波之间的对应关系一些在空间域直接表述非常困难，甚至是不可能的增强任务，在频率域中变的非常简单；通过频率域实验选择合适的滤波器，进行反变换获得空间滤波器，实际实施通常都是在空间域进行的。整理课件133主要内容遥感图像概述遥感图像的读写与显示遥感图像辐射增强遥感图像几何变换遥感图像配准遥感图像滤波遥感图像分割整理课件134对于给定的一幅含有多个物体的数字图像，模式识别的过程由三个阶段组成，如图所示图像分割特征抽取输入图像物体图像特征矢量分类物体类型“Bar”检测出各种物体，并把他们的图像和其余景物分离对物体进行度量，即对物体进行定量分析估计输出仅仅是一种决策，确定每个物体应该归属的类别整理课件135图像分割介绍定义将图像分成各具特性的区域，并提取出感兴趣目标的技术和过程灰度、颜色、纹理对应单个区域和多个区域图像处理过渡到图像分析的关键步骤，也是一种基本的计算机视觉技术借助集合概念进行正式的定义：令R代表整个图像区域，对R的分割可看做将R分成若干个满足以下条件的非空子集(子区域)R1,R2,R3…Rn整理课件136(1)∪i=1nRi=R分割所得到的全部子区域的总和(并集)应能包括图像中所有象素，或者说分割应将图像中的每个象素都分进某1个子区域中(2)对所有的i和j，i≠j，有Ri∩Rj=φ各个子区域是互不重叠的，或者说1个象素不能同时属于2各区域(3)i=1,2…n,有P(Ri)=TRUE在分割后得到的属于同1个区域中的象素应该具有某些相同特性(4)对i≠j,有P(Ri∪Rj)=FALSE在分割后得到的属于不同区域中的象素应该具有一些不同的特性(5)对i=1,2…n,Ri是连通的区域要求同1个子区域内的象素应当是连通的分割准则应可适用于所有区域和象素分割准则应能帮助确定各区域象素有代表性的特性整理课件137基于阈值的分割－通过阈值对不同物体进行分割基于边缘的分割－先确定边缘象素，并把它们连接在一起，以构成所需的边界基于区域的分割－把各象素划归到各个物体或区域中图像分割方法分类整理课件138基于阈值的分割图像阈值分割是一种广泛应用的分割技术，利用图像中要提取的目标物与其背景在灰度特性上的差异，把图像视为具有不同灰度级的两类区域(目标和背景)的组合，选取一个合适的阈值，以确定图像中每个象素点应该属于目标还是背景区域，从而产生相应的二值图像可以大量压缩数据，减少存储容量，而且能大大简化其后的分析和处理步骤整理课件139由此可见，阈值化分割算法主要有两个步骤：1、确定需要的分割阈值2、将分割阈值与象素值比较以划分象素在利用阈值方法来分割灰度图像时一般都对图像有一定的假设。基于一定的图像模型的。最常用的模型：假设图像由具有单峰灰度分布的目标和背景组成，处于目标或背景内部相邻象素间的灰度值是高度相关的，但处于目标和背景交界处两边的象素在灰度值上有很大的差别。如果一幅图像满足这些条件，它的灰度直方图基本上可看作是由分别对应目标和背景的两个单峰直方图混合构成的。整理课件140直方图阈值分割简单直方图分割法最佳阈值整理课件141简单直方图分割法图像的灰度级范围为0,1,…l-1，设灰度级i的象素数为ni，则一幅图像的总象素N为N=∑i=0l-1ni灰度级i出现的概率定义为：pi=ni/N灰度图像的直方图反映一幅图像上灰度分布的统计特性，成为利用象素灰度作属性的分割方法的基础整理课件142Z1ZiZtZjZk暗亮PB1B2背景目标60年代中期，Prewitt提出了直方图双峰法，即如果灰度级直方图呈明显的双峰状，则选取两峰之间的谷底所对应的灰度级作为阈值。整理课件143注意：应用灰度直方图双峰法来分割图像，也需要一定的图像先验知识，因为同一个直方图可以对应若干个不同的图像，直方图只表明图像中各个灰度级上有多少个象素，并不描述这些象素的任何位置信息。该方法不适合直方图中双峰差别很大或双峰间的谷比较宽广而平坦的图像，以及单峰直方图的情况。70年代初，研究工作集中在直方图变换，但无论是直方图还是直方图变换法都仅仅考虑了直方图灰度信息而忽略了图像的空间信息整理课件144最佳阈值所谓最佳阈值是指图像中目标物与背景的分割错误最小的阈值设一幅图像只由目标物和背景组成，已知其灰度级分布概率密度分布为P1(Z)和P2(Z)，且已知目标物象素占全图象素数比为θ，因此，该图像总的灰度级概率密度分布P(Z)可用下式表示：P(Z)=θP1(Z)+(1-θ)P2(Z)假定阈值为Z，认为图像由亮背景上的暗物体所组成，即灰度小于Z的位目标物，大于Z的为背景整理课件145P1(Z)P2(Z)Zt目标物背景如图所示，如选定Zt为分割阈值，则将背景象素错认为是目标物象素的概率为：E1(Zt)=∫-∞ZtP2(Z)dZE2(Zt)=∫zt∞P1(Z)dZ将目标物象素错认为是背景象素的概率为：因此，总的错误概率E(Z)为：E(Zt)=(1-θ)E1(Zt)+θE2(Zt)整理课件146最佳阈值就是使E(Zt)为最小值时的Zt，将E(Zt)对Zt求导，并令其等于0，解出其结果为：θP1(Zt)=(1-θ)P2(Zt)设P1(Zt)和P2(Zt)均为正态分布函数，其灰度均值分别为μ1和μ2，对灰度均值得标准偏差分别为σ1和σ2，即整理课件147将上两式代入，且对两边求对数，得到：简化为：AZt2+BZt+C=0上式是Zt的一个二次方程式，有两个解，因此，要使分割误差最小，需要设置两个阈值，即上式的两个解。如果设σ2=σ12=σ22，即方差相等，则上式方程存在唯一解，即：整理课件148如果设θ=1-θ，即θ＝1/2时，E1(Zt)E2(Zt)P1(Z)P2(Z)ZtZP从前面可以看出，假如图像的目标物和背景象素灰度级概率呈正态分布，且偏差相等(σ12=σ22)，背景和目标物象素总数也相等(θ＝1/2)，则这个图像的最佳分割阈值就是目标物和背景象素灰度级两个均值得平均整理课件149基于边缘的分割先检测不连续的点，然后将点连接成边界点检测线检测边缘检测无论哪种方法，其处理过程为R=w1z1+w2z2+…+w9z9=Swkz其中zk

是与模板系数wk相联系的灰度级象素

R代表模板中心象素的值w1w2w3w6w9w8w7w4w5整理课件150边缘检测两个具有不同灰度值的相邻区域之间总存在边缘边缘是灰度值不连续的结果，这种不连续常可利用求导数方便的检测到，一般常用一阶和二阶导数来检测边缘常见的边缘剖面图有三种，如图所示整理课件151图像水平方向剖面一阶导数二阶导数阶梯状－处于图像中2个具有不同灰度值的相邻区域之间脉冲状－对应细条状的灰度值突变区域屋顶状－上升下降沿都比较缓慢整理课件152一阶导数在图像由暗变明的位置处有1个向上的阶跃，而其它位置都为0，这表明可用一阶导数的幅度值来检测边缘的存在，幅度峰值一般对应边缘位置二阶导数在一阶导数的阶跃上升区有1个向上的脉冲，而在一阶导数的阶跃下降区有1个向下的脉冲，在这两个脉冲之间有1个过0点，它的位置正对应原图像中边缘的位置，所以可用二阶导数的过0点检测边缘位置，而用二阶导数在过0点附近的符号确定边缘象素在图像边缘的暗区或明区对(a、b)而言整理课件153对(c)而言，脉冲状的剖面边缘与(a)的一阶导数形状相同，所以(c)的一阶导数形状与(a)的二阶导数形状相同，而它的2个二阶导数过0点正好分别对应脉冲的上升沿和下降沿，通过检测脉冲剖面的2个二阶导数过0点就可确定脉冲的范围对(d)而言，屋顶状边缘的剖面可看作是将脉冲边缘底部展开得到，所以它的一阶导数是将(c)脉冲剖面的一阶导数的上升沿和下降沿展开得到的，而它的二阶导数是将脉冲剖面二阶导数的上升沿和下降沿拉开得到的，通过检测屋顶状边缘剖面的一阶导数过0点，可以确定屋顶位置整理课件154主要介绍以下几种边缘检测算子（1）梯度算子（2）方向算子（3）拉普拉斯算子（4）马尔算子（5）综合正交算子（6）坎尼算子整理课件155梯度算子梯度对应一阶导数，梯度算子是一阶导数算子，在边缘灰度值过渡比较尖锐且图像中噪声比较小时，梯度算子效果好。对1个连续图像函数f(x,y)，它在位置(x,y)的梯度可表示为1个矢量：整理课件156这个矢量的幅度(也常直接简称为梯度)和方向角分别为：有时用其它方式计算幅度，如：上面各式中的偏导数都需对一个象素位置计算，在实际中，常用小区域模板进行卷积近似计算整理课件157梯度运算比较复杂一点，在数字图像梯度运算过程中，可以按图像内容试用一些近似运算，以获得既能满足要求又能使运算简单的方法。对于数字图像，可用一阶差分代替一阶微分整理课件158在数字图像中，还经常使用Robert和Sobel等算子检测边缘Robert梯度采用的是对角方向相邻两象素之差，即1-11－1Roberts模板Sobel梯度算子先做加权平均，然后再微分，即整理课件159-1111-1-1111-1-1-1Prewitt模板-1121-1-2121-1-2-1Sobel模板水平水平垂直垂直Sobel算子是常用的，而且效果较其它两种算子好整理课件160实例原图水平方向垂直方向梯度

f

|Gx|+|Gy|整理课件161Sobel算子Roberts算子Prewitt算子原图整理课件162拉普拉斯算子一阶微分是一种矢量，不但有其大小，还有方向，和标量相比较，它数据存储量大，在具有相等斜率的宽区域上，有可能将全部区域都当作边缘提取出来Lapplacian算子是一种二阶导数算子，是不依赖于边缘方向的二阶微分算子，是一个标量而不是矢量，具有旋转不变即各向同性的性质，对一个连续函数f(x,y)，它在图像中位置(x,y)的拉普拉斯值定义为：整理课件163拉普拉斯算子是无方向性的算子，它比前述计算多个方向导数算子的计算量要小，因为只需用一个模板，且不必综合各模板的值。在数字图像中，计算函数的拉普拉斯也可以借助各种模板卷积实现。这里对模板的基本要求是对应中心象素的系数应是正的，而对应中心象素邻近象素的系数应是负的，且所有系数的和应为0，这样就不会产生灰度偏移-1-1-1-14-1-1-1-14-1-1-1-1-1-1-1-18在数字图像情况下的近似为：整理课件164拉普拉斯是一种二阶导数算子，所以对图像中的噪声相当敏感。在实际中，常常在进行平滑操作地同时进行二阶微分，如图所示是具有平滑效果的二阶微分算子操作。另外它常产生双象素宽的边缘，且也不能提供边缘方向的信息。由于以上原因，拉普拉斯算子很少直接用于边缘检测，而主要用于已知边缘象素后，确定该象素是在图像的暗区或明区一边。另一方面，一阶差分算子会在较宽范围形成较大的梯度值，因此不适合于精确定位，而利用二阶差分算子过0点可以精确定位边缘111111111-4-4-4-4-4-4-4-4-4111111111111111111111111111整理课件165基于区域的分割图像分割－把图像分解为若干个有意义的子区域，而这种分解－基于物体有平滑均匀的表面，与图像中强度恒定或缓慢变化的区域相对应，即每个子区域都具有一定的均匀性质前面所讨论的边缘、阈值，没有明显使用分割定义中的均匀测度度量区域分割－直接根据事先确定的相似性准则，直接取出若干特征相近或相同象素组成区域常用的区域分割－区域增长(区域生长)、区域分裂－合并方法等整理课件166区域增长原理和步骤基本思想－将具有相似性质的象素集合起来构成区域。具体步骤－先对每个需要分割的区域找一个种子象素作为生长起点，然后将种子象素周围邻域中与种子象素有相同或相似性质的象素(根据某种事先确定的生长或相似准则来判定)合并到种子象素所在的区域中。将这些新象素当做新的种子象素继续进行上面的过程，直到再没有满足条件的象素可被包括进来，这样一个区域就长成了●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●R1R5R4R3R2整理课件167如图给出已知种子点区域生长的一个示例。1047510477015552056522564115551155511555115551155511575115771155521555225551111111111111111111111111(a)(b)(c)(d)(a)给出需要分割的图像，设已知两个种子象素(标为深浅不同的灰色方块)，现在进行区域生长采用的判断准则是：如果所考虑的象素与种子象素灰度值差的绝对值小于某个门限T，则将该象素包括进种子象素所在的区域图(b)给出T=3时区域生长的结果，整幅图被较好的分成2个区域图(c)给出T=1时区域生长的结果，有些象素无法判定图(d)给出T=6时区域生长的结果，整幅图都被分成1个区域。整理课件168从上面