科学计算与MATLAB语言：第五章 MATLAB在遥感图像处理中的应用

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

2、。为了便于应用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)，称之为图像信号，而前面几个式子称为图像，以示区别。数字图像是由被称作象素的小块区域组成的二维矩阵。将物理图象行列划分后，每个小块区域称为像素（pixel）。 每个像素包括两个属性：位置和灰度。对于单色即灰度图像而言，每个象素的亮度用一个数值来表示，通常数值范围在0到255之间，即可用一个字节来表示， 0表示黑、255表示白，而其它表示灰度级别。 物理图象及对应 的数字图象物理图像

4、19643灰度像素数字图像采样列采样行图片像素行间隔采样列间隔灰阶黑灰白0128255彩色图象可以用红、绿、蓝三元组的二维矩阵来表示。 通常，三元组的每个数值也是在0到255之间，0表示相应的基色在该象素中没有，而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) (22

5、6,144,133) (224,137,124) (227,151,136) (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,

6、120) (159, 51, 71) (189, 89,101) (216,111,110) (217,124,121) (227,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) (

7、216,179,170) (221,184,170) (190, 89, 89) (204,109,113) (204,115,118) (189, 85, 97) (159, 60, 78) (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等格式的图像，可以直接利用MATL

8、AB提供的命令imread及imwrite进行读写操作。实际中不同的遥感数据并非标准格式。如何根据数据提供者提供的格式说明文件读出相应的参数及图像数据是进行后续遥感图像处理的关键。标准格式图像读取函数imread可以从任何Matlab支持的图像文件格式中，以任意位深度读取一幅图像。格式为：X,MAP=imread(FILENAME,FMT)，其中：FILENAME为需要读入的图像文件名称，FMT为图像格式。 图像文件格式读取图像信息可以通过调用imfinfo函数获得与图像文件有关的信息，格式如下：INFO=imfinfo(FILENAME,FMT)其中：返回的INFO是Matlab的一个结构体

9、。函数size可给出一副图像的行数和列数size（f） 非标准格式遥感图像读取通常每一种遥感卫星传感器都制定了一套适合本身特性的数据存储格式标准。并且为用户提供详细的格式说明书。用户必须熟悉格式才能正确地读出遥感图像数据。例：ALOS PALSAR参数及数据读取图像的显示Matlab的图像显示函数主要有imshow、colorbar、subimage和imagesc等。 例如：imshow函数的基本语法为： imshow(f,G) 其中，f是一个图像数组，G是显示该图像的灰度级数。1、imshow(I,n)显示灰度级为n的图像，n缺省为256。【例】按256灰度级显示I=imread(moon

10、.tif);imshow(I,256)colorbar 【例】按8灰度级显示I=imread(moon.tif);imshow(I,8)colorbar 【例】按最大灰度范围显示I=imread(moon.tif);imshow(I,)colorbar 【例】按指定灰度范围显示I=imread(moon.tif);imshow(I,64,128)colorbar 2、imshow(X,MAP)显示索引图像，X为数据图像矩阵，MAP为调色板。【例】显示索引图像X,MAP=imread(canoe.tif);imshow(X,MAP);colorbar 【例】显示伪彩色图像I=imread(moo

11、n.tif);imshow(I,jet(256)colorbar3、imshow(RGB)显示真彩色图像。【例】rgb=imread(flowers.tif);imshow(rgb); 4、subimage多图显示多个调色板的图像。【例】load trees;x2,map2=imread(forest.tif);subplot(2,1,1),subimage(X,map);colorbarsubplot(2,1,2),subimage(x2,map2);colorbar 【例】不用专用函数显示多图load trees;x2,map2=imread(forest.tif);subplot(2,1

12、,1),imshow(x2,map2);colorbarsubplot(2,1,2),imshow(X,map);colorbar 图像的写入1、imwrite函数imwrite(I,文件名,文件格式)保存无调色板的图像imwrite(X,MAP,文件名,文件格式)保存有调色板的图像imwrite(.,文件名)文件名中含格式imwrite(.,参数,值)指定保存参数文件格式：bmp、jpg、pcx、tif。 【例】将tif图像保存为jpg图像x,map=imread(canoe.tif);imwrite(x,map,canoe.jpg,JPG,Quality,75) 主要内容遥感图像概述遥感图

13、像的读写与显示遥感图像辐射增强遥感图像几何变换遥感图像配准遥感图像滤波 遥感图像分割3 遥感图像辐射增强定义：将原来不清晰的图像变得清晰或突出某些特征，同时抑制一些不需要的信息的处理方法。目的：突出图像中的有用信息，扩大不同影像特征之间的差别，以便于进行判读和分析。基于直方图变换的增强方法直方图变换是一种通过直接改变图像中像元的亮度值来改变图像的对比度，从而改善图像质量的图像处理方法。灰度直方图是灰度级的函数，描述图像中具有该灰度级的像元的个数。其横坐标是灰度级，纵坐标是像元的个数（或该灰度出现的频率）。直方图反映灰度的总体结构，但是不反映空间的分布信息。直方图变换的类型线性变换非线性变换线性

14、变换：根据原图像直方图来确定好拉伸变换前的灰度值区间，然后把这一灰度值区间按某一直线方程关系拉伸或压缩而成为变换后灰度值区间。拉伸后的图像灰度值范围增大，对比度改善。线性变换按比例扩大原是灰度级的范围，将原始的相对集中的灰度值分布在0255范围内展开。线性变换分段线性拉伸对不同范围的灰度值进行不同的拉伸。非线性变换变换函数是非线性的，如指数变换、对数变换等。指数变换指数变换可以对图像高值区域进行拉伸对数变换对数变换对图像低值区域进行拉伸直方图均衡将原图像的直方图通过变换函数变为均匀的直方图，然后按均匀直方图修改原图像，从而获得一幅灰度分布均匀的新图像。直方图均衡计算步骤1. 计算原图像f的直方

15、图h直方图均衡计算步骤2. 求出图像f的总体像素个数N，计算每个灰度级的像素个数在整个图像中所占的比例hs。直方图均衡计算步骤3. 计算图像各灰度级的累积分布hp。直方图均衡计算步骤4. 求出新图像g的灰度值。因为直方图是近似的概率密度函数，所以用离散灰度级作变换一般得不到完全平坦的结果。另外，从上例可以看出，变换后的灰度级减少了，这种现象叫做“简并”现象。由于简并现象的存在，处理后的灰度级总是要减少的，这是像素灰度有限的必然结果。由于上述原因，数字图像的直方图均衡只是近似的。主要内容遥感图像概述遥感图像的读写与显示遥感图像辐射增强遥感图像几何变换遥感图像配准遥感图像滤波 遥感图像分割遥感图像

16、几何变换几何运算可以看成是像素在图像内的移动过程，该移动过程可以改变图像中物体对象(像素)之间的空间关系。几何运算可以是不受任何限制的，但是通常都需要做出一些限制以保持图像的外观顺序。完整的几何运算需要由两个算法来实现：空间变换算法和灰度插值算法。空间变换主要用来保持图像中曲线的连续性和物体的连通性，一般都采用数学函数形式来描述输入、输出图像相应像素间的空间关系。空间变换的一般定义为 其中，f表示输入图像，g表示输出图像，坐标(x,y)指的是空间变换后的坐标，要注意这时的坐标已经不是原来的坐标(x,y)了。a(x,y)和b(x,y)分别是图像的x和y坐标的空间变换函数。 灰度级插值主要是对空间

17、变换后的像素赋予灰度值，使之恢复原位置处的灰度值。在几何运算中，灰度级插值是必不可少的组成部分，因为图像一般用整数位置处的像素来定义。而在几何变换中，g(x,y)的灰度值一般由处在非整数坐标上的f(x,y)的值来确定，即g中的一个像素一般对应于f中的几个像素之间的位置，反过来看也是一样，即f中的一个像素往往被映射到g中的几个像素之间的位置。 图像的几何常用函数图像的缩放 imresize图像的旋转 imrotate图像的剪裁 imcrop图像的一般几何变换 imtransform 图像的插值运算最近邻插值方法（Nearest neighbor）双线性插值方法（Bilinear）双三次插值方法（

18、Bicubic）(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) (2) 图像旋转使用imrotate函数来旋转一幅图像。调用格式如下：B=imrotate(A,

19、ANGLE,METHOD,BBOX)，其中：A需要旋转的图像；ANGLE表示旋转的角度，正值为逆时针；METHOD插值方法；BBOX 取值loose（缺省）、crop【例】图像旋转clfI,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) (3) 图像剪切使用imcrop函数可以从一幅图像中抽取一个矩形的部分。imcr

20、op函数的调用格式如下：X2=imcrop(X,MAP,RECT) 其中，X表示有待剪切的图像，不指定X时，imcrop将当前坐标轴中的图像作为待剪切的图像。MAP表示X为索引图像时的调色板，RECT定义剪切区的矩形坐标。如果调用imcrop时不指定矩形的坐标，那么当光标位于图像中时会变成十字形，可以通过拖曳鼠标的方式交互式地选择一个矩形。imcrop函数根据用户的选择绘制一个矩形，释放鼠标键后将产生一个新的图像。 仿射变换仿射变换，可以用以下函数来描述：其中，A是变形矩阵，b是平移矩阵。 (1) 尺度变换变换矩阵： ,S0 【例】clf;I=checkerboard(20,2);subplo

21、t(121);imshow(I);axis on;title(原图)s=1.5;T=s 0;0 s;0 0;tf=maketform(affine,T);I1=imtransform(I,tf,bicubic,FillValues,0.3);subplot(122);imshow(I1);axis on;title(尺度变换) (2) 伸缩变换变换矩阵： 【例】clf;I=checkerboard(20,2);subplot(121);imshow(I);axis on;title(原图)t=2;T=1 0;0 t;0 0;tf=maketform(affine,T);I1=imtransfo

22、rm(I,tf,bicubic,FillValues,0.3);subplot(122);imshow(I1);axis on;title(伸缩变换) (3) 扭曲变换变换矩阵： 【例】clf;I=checkerboard(20,2);subplot(121);imshow(I);axis on;title(原图)u=0.5;T=1 u;0 1;0 0;tf=maketform(affine,T);I1=imtransform(I,tf,bicubic,FillValues,0.3);subplot(122);imshow(I1);axis on;title(扭曲变换) (4) 旋转变换变换矩

23、阵： clf;I = checkerboard(20,2);subplot(1,2,1);imshow(I);title(原图)angle=20*pi/180;sc=cos(angle);ss=sin(angle);T=sc -ss; ss sc;0 0;tf=maketform(affine,T);I1=imtransform(I,tf,bicubic,FillValues,0.3);subplot(122);imshow(I1);title(旋转变换) (5) 综合变换变换矩阵： 【例】clf;I = checkerboard(20,2);subplot(1,2,1);imshow(I);

24、title(原图)s=2;As=s 0;0 s; % 尺度t=2;At=1 0;0 t; % 伸缩u=1.5;Au=1 u;0 1; % 扭曲st=30*pi/180;sc=cos(st);ss=sin(st);Ast=sc -ss; ss sc; % 旋转T=As*At*Au*Ast;3 5;tf=maketform(affine,T);I1=imtransform(I,tf,bicubic,FillValues,0.3);subplot(122);imshow(I1);title(综合) 主要内容遥感图像概述遥感图像的读写与显示遥感图像辐射增强遥感图像几何变换遥感图像配准遥感图像滤波 遥感

25、图像分割图像配准是指依据一些相似性度量决定图像间的变换参数使得从不同传感器不同视角不同时间获取的同一场景的两幅或多幅图像变换到统一坐标系下在像素层上得到最佳匹配的过程。待配准图像相对于参考图像的配准可定义为两幅图像在空间和亮度上的映射两幅图像可定义为两个二维数组分别用I1(x,y)和I2(x,y)表示他们分别是两幅图像的亮度值则两幅图像间的映射可表示为：其中：f为二维空间坐标变换，g为一维亮度或其他亮度变换遥感图像配准步骤第一步：读取图像orthophoto = imread(westconcordorthophoto.png);figure, imshow(orthophoto)unregi

26、stered = imread(westconcordaerial.png);figure, imshow(unregistered)第2步：选取控制点cpselect(unregistered, orthophoto)第3步：保存控制点到MATLAB工作空间中File menu - choose the Export Points to Workspace第4步(可选)：控制点微调以提高配准精度cpcorr(?)第5步：指定变换类型及其参数mytform = cp2tform(input_points, base_points, projective);第6步：转换待配准图像register

27、ed = imtransform(unregistered, mytform);imshow(orthophoto)hold onh = imshow(registered, gray(256);set(h, AlphaData, 0.6)主要内容遥感图像概述遥感图像的读写与显示遥感图像辐射增强遥感图像几何变换遥感图像配准遥感图像滤波 遥感图像分割遥感图像滤波方法分类空间域滤波频率域滤波图像的空间域滤波滤波是一种图像修正或增强技术。可以突出图像的某些特征，也可以删除另一些特征。图像滤波的本质是一种邻域操作，输出图像的任一个像素值都是通过输入图像对应的像素邻域内的像素值利用一定的算法得到的。图像

28、的空间域滤波就是输出图像任一个像素值都是通过输入图像对应的像素邻域内的像素值线性组合得到的。图像的空间域滤波二维卷积函数 conv2滤波函数 filter2图像滤波 imfilter预定义滤波器 fspecial 设当前待处理像素为f (m,n) ，给出一个大小为33的处理模板。以模块运算系数表示即：常用的模板：【例】x=imread(cameraman.tif);A=1/25*ones(5,5);x2=conv2(double(x),A);imshow(x)figureimshow(int8(x2)B=0.5*0 0.25 0; 0.25 1 0.25;0 0.25 0;x3=conv2(d

29、ouble(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)卷积步骤：关于卷积核的中心，旋转卷积核180度。滑动卷积核，将卷积核的中心位于图像矩阵的每一个元素。将旋转后的卷积核作为权重，乘以对应的矩阵元素求加权和计算卷积输

30、出(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 default),返回全部二维卷积结果shape=same, 返回与A同样大小的卷积中心部分shape=valid, 不考虑边界补零，返回C大小为（ma-mb+1）（na-nb+1）相关运算实现滤波：filter2相关运算(correlation)滑动相关核，将相关核的中心位于图像矩阵的每一个元素。将相关核作为权重，乘以对应的矩阵元素

31、求加权和用相关运算实现滤波的函数是filter2B = filter2(h,A)计算相关输出(2,4)图像的线性滤波imfilterB=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), title(Original Image);figure, imshow

32、(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既可以用相关，也可以用卷积实现滤波操作，缺省是相关。A = magic(5);h = -1 0 1imfilter(A,h) % filter using correlationimfilter(A,h,conv) % filter using con

33、volutionoption: 边界补零(zero-padding)和边界复制(border replication)边界补零(zero-padding): 缺省option: 边界补零(zero-padding)和边界复制(border replication)边界补零(zero-padding): 缺省I = imread(eight.tif);h = ones(5,5) / 25;I2 = imfilter(I,h);imshow(I), title(Original Image);figure, imshow(I2), title(Filtered Image with Black B

34、order)缺点：滤波得到的图像边缘有一个dark bandoption: 边界补零(zero-padding)和边界复制(border replication)边界复制(border replication)option: 边界补零(zero-padding)和边界复制(border replication)边界复制(border replication) I3 = imfilter(I,h,replicate);figure, imshow(I3);title(Filtered Image with Border Replication)imfilter还有其他的边界补充选项，参考imfi

35、lter的帮助多维滤波imfilter既可以处理多维图像，也可以处理多维滤波器。用一个二维滤波器对一个三维图像滤波，相当于对三维图像的每个平面进行二维滤波。例 ：用同样的滤波器对一个真彩色图像的每个颜色平面进行滤波。% 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的

36、关系：filter2、conv2、 convn ：将输入转换为double类型，输出也是double的，输入总是补零（zero padded), 不支持其他的边界补充选项。imfilter：不将输入转换为double，输出只与输入同类型，有灵活的边界补充选项主要内容空间域滤波频率域滤波背景 Background法国数学家傅立叶(生于1768年)在1822年出版的热分析理论一书中指出:任何周期函数都可以表达为不同频率的正弦和或余弦和的形式,即傅立叶级数。20世纪50年代后期,快速傅立叶变换算法出现,得到了广泛的应用。傅立叶变换(Fourier Transform)傅里叶变换和频率域的介绍一维傅立

37、叶变换及其反变换二维DFT变换及其反变换二维DFT变换性质一维傅立叶变换及其反变换连续函数f(x)的傅立叶变换F(u):傅立叶变换F(u)的反变换:一维傅立叶变换及其反变换离散函数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。一维离散傅里叶变换离散傅里叶变换及其反变换总存在。用欧拉公式得每个F(u) 由f(x)与对应频率的正弦和余弦乘积和组成； u 值决定了变换的频率成份，因此，F(u) 覆盖的域(u值) 称为

38、频率域，其中每一项都被称为FT 的频率分量。与f(x) 的“时间域”和“时间成份”相对应。一维离散傅里叶变换傅里叶变换将信号分成不同频率成份。类似光学中的分色棱镜把白光按波长(频率)分成不同颜色，称数学棱镜。一维离散傅里叶变换傅立叶变换在极坐标下表示:频率谱相位谱功率谱一维离散傅里叶变换f(x)是一门函数，如图所示，它表示为：求其傅立叶变换F(u)一维离散傅里叶变换解：一维离散傅里叶变换对应的傅立叶谱为：一维离散傅里叶变换简单函数的傅里叶谱M 点离散函数及其傅里叶频谱(M=1024, A=1, K=8)； 对应的傅里叶频谱 曲线下面积：当x 域加倍时，频率谱的高度也加倍； 当函数长度加倍时，相

39、同间隔下频谱中零点的数量也加倍。二维DFT傅里叶变换一个图像尺寸为MN的函数f(x,y)的离散傅立叶变换F(u,v):F(u,v)的反变换的反变换:二维DFT傅里叶变换二维离散傅立叶变换在极坐标下表示:频率谱相位谱功率谱二维DFT傅里叶变换(u,v)=(0,0)位置的傅里叶变换值为即f(x,y) 的均值，原点(0,0) 的傅里叶变换是图像的平均灰度。F(0,0) 称为频率谱的直流分量(系数)，其它F(u,v) 值称为交流分量(交流系数)。离散傅里叶变换是对区间0,M-1 中的u 值表述的，变换结果是关于原点对称的两个半周期，要显示完全的周期，需要将变换的原点移到u=M/2，二维图像中心化亦是如

40、此简单二维函数的中心谱空间域和频率域抽样点之间的关系如下：简单二维函数的中心谱MATLAB中离散傅立叶(DFT)的计算：fft、fft2、fftn：分别实现一维、二维和n维离散快速傅立叶变换(DFT)ifft、ifft2、ifftn：分别实现一维、二维和n维离散快速傅立叶逆变换(IDFT)例：图像傅立叶变换和幅值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, imsho

41、w(F2,-1 5,InitialMagnification,fit);colormap(jet); colorbar上图：傅立叶变换的取样很粗，且零频率系数不像传统的那样显示在图像中心，而是显示在左上角。为获得傅立叶变换较精细的取样，采用补零的方法。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)

42、;figure,imshow(log(abs(F2),-1 5,InitialMagnification,fit);colormap(jet); colorbar频率域滤波频率域的基本性质 每个F(u,v)项包含了被指数项修正的f(x,y)的所有值：直观上将傅里叶变换和图像中的亮度变化联系起来并不困难： 直流分量F(0,0)对应一幅图像的平均灰度； 低频部分对应图像缓慢变化的分量； 高频部分对应图像边缘和灰度级突变的部分图为一幅集成电路的扫描电子显微镜(SEM)图像，放大将近2500倍。注意图中45的强边缘，和两个因热感应不足而产生的白色氧化突起图是上图的傅里叶频谱，沿着45方向对应上图边缘突

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

44、 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=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);ims

45、how(f)figureimshow(gf)一些基本的滤波器及其性质低通滤波器：使低频通过，高频衰减低频主要决定图像在平滑区域中总体灰度级的显示比原始图像少一些尖锐的细节部分高通滤波器：使高频通过，低频衰减高频决定图像细节部分，如边缘和噪声在平滑区域中减少灰度级变化，突出过渡（如边缘） 灰度级的细节部分，使图像更加锐化。基本的滤波器及其性质图像被模糊锐化F(0,0)=0:几乎没有平滑细节周期对称Lowpass filterHighpass filter 空间域滤波和频率域滤波之间的对应关系卷积定理： 空间域的乘法对应频域卷积空间域滤波和频率域滤波之间的对应关系重要性质:根据冲击函数和卷积定理的

46、性质，可知空间域和频率域的滤波器组成傅里叶变换对h(x,y)和H(u,v) 。给出频率域滤波器H(u,v) ，通过反傅里叶变换可以得到空间域相应的滤波器h(x,y) 。空间域滤波和频率域滤波之间的对应关系滤波器大小前述的所有函数均具有相同的尺寸MN。在实际中，指定一个频率域滤波器，进行反变换后会得到一个相同尺寸的空间域滤波器。如果两个域中滤波器尺寸相同，那么通常频域中进行滤波计算更为有效，更为直观，但空域中适用更小尺寸的滤波器，更为有效。空间域滤波和频率域滤波之间的对应关系方法：在频率域指定滤波器；做反变换；使用结果滤波器作为在空间域构建更小空间滤波模板的指导；基于高斯函数的滤波高斯滤波器函数

47、(低通):对应的空间域滤波器 为高斯曲线的标准差。基于高斯函数的滤波高斯滤波器的重要特性频域和空域高斯滤波器构成傅里叶变换对，且都是实高斯函数。处理时不用考虑复数，而且高斯曲线直观，易于操作。高斯滤波器傅里叶变换对之间有相互作用： 当H(u)有很宽轮廓时(大值)，h(x)很窄轮廓，反之亦然； 当趋于无穷时，H(u) 趋于常函数，h(x) 趋于冲击函数。基于高斯函数的滤波器频域滤波器越窄，滤除的低频部分越多，图像越模糊。在空域中意味着滤波器越宽，模板就越大(阶数越大)。空间域滤波和频域滤波之间的对应关系一些在空间域直接表述非常困难，甚至是不可能的增强任务，在频率域中变的非常简单；通过频率域实验选

48、择合适的滤波器，进行反变换获得空间滤波器，实际实施通常都是在空间域进行的。主要内容遥感图像概述遥感图像的读写与显示遥感图像辐射增强遥感图像几何变换遥感图像配准遥感图像滤波 遥感图像分割对于给定的一幅含有多个物体的数字图像，模式识别的过程由三个阶段组成，如图所示图像分割特征抽取输入图像物体图像特征矢量分类物体类型“Bar”检测出各种物体，并把他们的图像和其余景物分离对物体进行度量，即对物体进行定量分析估计输出仅仅是一种决策，确定每个物体应该归属的类别图像分割介绍定义将图像分成各具特性的区域，并提取出感兴趣目标的技术和过程灰度、颜色、纹理对应单个区域和多个区域图像处理过渡到图像分析的关键步骤，也是

49、一种基本的计算机视觉技术借助集合概念进行正式的定义：令R代表整个图像区域，对R的分割可看做将R分成若干个满足以下条件的非空子集(子区域)R1,R2,R3Rn(1) i=1nRi=R分割所得到的全部子区域的总和(并集)应能包括图像中所有象素，或者说分割应将图像中的每个象素都分进某1个子区域中(2) 对所有的i和j，ij，有RiRj=各个子区域是互不重叠的，或者说1个象素不能同时属于2各区域(3) i=1,2n,有P(Ri)=TRUE在分割后得到的属于同1个区域中的象素应该具有某些相同特性(4) 对ij,有P(RiRj)=FALSE在分割后得到的属于不同区域中的象素应该具有一些不同的特性(5) 对

50、i=1,2n,Ri是连通的区域要求同1个子区域内的象素应当是连通的分割准则应可适用于所有区域和象素分割准则应能帮助确定各区域象素有代表性的特性基于阈值的分割通过阈值对不同物体进行分割基于边缘的分割先确定边缘象素，并把它们连接在一起，以构成所需的边界基于区域的分割把各象素划归到各个物体或区域中图像分割方法分类基于阈值的分割图像阈值分割是一种广泛应用的分割技术，利用图像中要提取的目标物与其背景在灰度特性上的差异，把图像视为具有不同灰度级的两类区域(目标和背景)的组合，选取一个合适的阈值，以确定图像中每个象素点应该属于目标还是背景区域，从而产生相应的二值图像可以大量压缩数据，减少存储容量，而且能大大

51、简化其后的分析和处理步骤由此可见，阈值化分割算法主要有两个步骤：1、确定需要的分割阈值2、将分割阈值与象素值比较以划分象素在利用阈值方法来分割灰度图像时一般都对图像有一定的假设。基于一定的图像模型的。最常用的模型：假设图像由具有单峰灰度分布的目标和背景组成，处于目标或背景内部相邻象素间的灰度值是高度相关的，但处于目标和背景交界处两边的象素在灰度值上有很大的差别。如果一幅图像满足这些条件，它的灰度直方图基本上可看作是由分别对应目标和背景的两个单峰直方图混合构成的。阈值选取方法直方图阈值分割法类间方差阈值分割法二维最大熵值分割法直方图阈值分割简单直方图分割法最佳阈值简单直方图分割法图像的灰度级范围

52、为0,1,l-1，设灰度级i的象素数为ni，则一幅图像的总象素N为N=i=0l-1ni灰度级i出现的概率定义为：pi=ni/N灰度图像的直方图反映一幅图像上灰度分布的统计特性，成为利用象素灰度作属性的分割方法的基础Z1ZiZt Zj Zk暗亮PB1B2背景目标60年代中期，Prewitt提出了直方图双峰法，即如果灰度级直方图呈明显的双峰状，则选取两峰之间的谷底所对应的灰度级作为阈值。注意：应用灰度直方图双峰法来分割图像，也需要一定的图像先验知识，因为同一个直方图可以对应若干个不同的图像，直方图只表明图像中各个灰度级上有多少个象素，并不描述这些象素的任何位置信息。该方法不适合直方图中双峰差别很大

53、或双峰间的谷比较宽广而平坦的图像，以及单峰直方图的情况。70年代初，研究工作集中在直方图变换，但无论是直方图还是直方图变换法都仅仅考虑了直方图灰度信息而忽略了图像的空间信息最佳阈值所谓最佳阈值是指图像中目标物与背景的分割错误最小的阈值设一幅图像只由目标物和背景组成，已知其灰度级分布概率密度分布为P1(Z)和P2(Z)，且已知目标物象素占全图象素数比为，因此，该图像总的灰度级概率密度分布P(Z)可用下式表示： P(Z)= P1(Z)+(1-)P2(Z)假定阈值为Z，认为图像由亮背景上的暗物体所组成，即灰度小于Z的位目标物，大于Z的为背景P1(Z)P2(Z)Zt目标物背景如图所示，如选定Zt为分割

54、阈值，则将背景象素错认为是目标物象素的概率为：E1(Zt)=-ZtP2(Z)dZE2(Zt)=zt P1(Z)dZ将目标物象素错认为是背景象素的概率为：因此，总的错误概率E(Z)为：E(Zt)=(1-)E1(Zt)+ E2(Zt)最佳阈值就是使E(Zt)为最小值时的Zt，将E(Zt)对Zt求导，并令其等于0，解出其结果为：P1(Zt)=(1- )P2(Zt)设P1(Zt)和P2(Zt)均为正态分布函数，其灰度均值分别为1和2，对灰度均值得标准偏差分别为1和2，即将上两式代入，且对两边求对数，得到：简化为：AZt2+BZt+C=0上式是Zt的一个二次方程式，有两个解，因此，要使分割误差最小，需要

55、设置两个阈值，即上式的两个解。如果设2= 12 = 22，即方差相等，则上式方程存在唯一解，即：如果设=1- ，即1/2时，E1(Zt)E2(Zt)P1(Z)P2(Z)ZtZP从前面可以看出，假如图像的目标物和背景象素灰度级概率呈正态分布，且偏差相等(12 = 22)，背景和目标物象素总数也相等(1/2)，则这个图像的最佳分割阈值就是目标物和背景象素灰度级两个均值得平均基于边缘的分割先检测不连续的点，然后将点连接成边界点检测线检测边缘检测无论哪种方法，其处理过程为 R = w1z1+w2z2+w9z9 = S wkz其中 zk 是 与模板系数wk相联系的灰度级象素 R代表模板中心象素的值w1w

56、2w3w6w9w8w7w4w5边缘检测两个具有不同灰度值的相邻区域之间总存在边缘边缘是灰度值不连续的结果，这种不连续常可利用求导数方便的检测到，一般常用一阶和二阶导数来检测边缘常见的边缘剖面图有三种，如图所示图像水平方向剖面一阶导数二阶导数阶梯状处于图像中2个具有不同灰度值的相邻区域之间脉冲状对应细条状的灰度值突变区域屋顶状上升下降沿都比较缓慢一阶导数在图像由暗变明的位置处有1个向上的阶跃，而其它位置都为0，这表明可用一阶导数的幅度值来检测边缘的存在，幅度峰值一般对应边缘位置二阶导数在一阶导数的阶跃上升区有1个向上的脉冲，而在一阶导数的阶跃下降区有1个向下的脉冲，在这两个脉冲之间有1个过0点，

57、它的位置正对应原图像中边缘的位置，所以可用二阶导数的过0点检测边缘位置，而用二阶导数在过0点附近的符号确定边缘象素在图像边缘的暗区或明区对(a、b)而言对(c)而言，脉冲状的剖面边缘与(a)的一阶导数形状相同，所以(c)的一阶导数形状与(a)的二阶导数形状相同，而它的2个二阶导数过0点正好分别对应脉冲的上升沿和下降沿，通过检测脉冲剖面的2个二阶导数过0点就可确定脉冲的范围对(d)而言，屋顶状边缘的剖面可看作是将脉冲边缘底部展开得到，所以它的一阶导数是将(c)脉冲剖面的一阶导数的上升沿和下降沿展开得到的，而它的二阶导数是将脉冲剖面二阶导数的上升沿和下降沿拉开得到的，通过检测屋顶状边缘剖面的一阶导

58、数过0点，可以确定屋顶位置主要介绍以下几种边缘检测算子（1）梯度算子（2）方向算子（3）拉普拉斯算子（4）马尔算子（5）综合正交算子（6）坎尼算子梯度算子梯度对应一阶导数，梯度算子是一阶导数算子，在边缘灰度值过渡比较尖锐且图像中噪声比较小时，梯度算子效果好。对1个连续图像函数f(x,y)，它在位置(x,y)的梯度可表示为1个矢量：这个矢量的幅度(也常直接简称为梯度)和方向角分别为：有时用其它方式计算幅度，如：上面各式中的偏导数都需对一个象素位置计算，在实际中，常用小区域模板进行卷积近似计算梯度运算比较复杂一点，在数字图像梯度运算过程中，可以按图像内容试用一些近似运算，以获得既能满足要求又能使运

59、算简单的方法。对于数字图像，可用一阶差分代替一阶微分在数字图像中，还经常使用Robert和Sobel等算子检测边缘Robert梯度采用的是对角方向相邻两象素之差，即1-111Roberts模板Sobel梯度算子先做加权平均，然后再微分，即-1111-1-1111-1-1-1Prewitt模板-1121-1-2121-1-2-1Sobel模板水平水平垂直垂直Sobel算子是常用的，而且效果较其它两种算子好实例原图水平方向垂直方向梯度f |Gx| + |Gy|Sobel算子Roberts算子Prewitt算子原图拉普拉斯算子一阶微分是一种矢量，不但有其大小，还有方向，和标量相比较，它数据存储量大，

60、在具有相等斜率的宽区域上，有可能将全部区域都当作边缘提取出来Lapplacian 算子是一种二阶导数算子，是不依赖于边缘方向的二阶微分算子，是一个标量而不是矢量，具有旋转不变即各向同性的性质，对一个连续函数f(x,y)，它在图像中位置(x,y)的拉普拉斯值定义为：拉普拉斯算子是无方向性的算子，它比前述计算多个方向导数算子的计算量要小，因为只需用一个模板，且不必综合各模板的值。在数字图像中，计算函数的拉普拉斯也可以借助各种模板卷积实现。这里对模板的基本要求是对应中心象素的系数应是正的，而对应中心象素邻近象素的系数应是负的，且所有系数的和应为0，这样就不会产生灰度偏移-1-1-1-14-1-1-1