求在matlab上实现gabor滤波器对图像处理的程序

这是别人写的一个gabor滤波器的程序，你看看吧！

% GABORFILTER Bi-dimensional Gabor filter with DC component compensation.

% [G,GABOUT]=GABORFILTER(I,S,F,W,P) filters the input image I with the 2D

% Gabor filter G described by the parameters S, F, W and P to create the

% output filtered image GABOUT.

% This version of the 2D Gabor filter is basically a bi-dimensional

% Gaussian function centered at origin (0,0) with variance S modulated by

% a complex sinusoid with polar frequency (F,W) and phase P described by

% the following equation:

%

% G(x,y,S,F,W,P)=k*Gaussian(x,y,S)*(Sinusoid(x,y,F,W,P)-DC(F,S,P)),

% where:

% Gaussian(x,y,S)=exp(-pi*S^2*(x^2+y^2))

% Sinusoid(x,y,F,W,P)=exp(j*(2*pi*F*(x*cos(W)+y*sin(W))+P)))

% DC(F,S,P)=exp(-pi*(F/S)^2+j*P)

%

% PS: The term DC(F,S,P) compensates the inherent DC component produced

% by the Gaussian envelop as shown by Movellan in [1].

%

% Tips:

% 1) To get the real part and the imaginary part of the complex

% filter output use real(gabout) and imag(gabout), respectively

%

% 2) To get the magnitude and the phase of the complex filter output

% use abs(gabout) and angle(gabout), respectively.

% Author: Stiven Schwanz Dias e-mail: [email protected]

% Cognition Science Group, Informatic Department,

% University of Esp韗ito Santo, Brazil, January 2007.

%

% References:

% [1] Movellan, J. R. - Tutorial on Gabor Filters. Tech. rep., 2002.

function [G,GABOUT]=gaborfilter(I,S,F,W,P)

if isa(I,'double')~=1

I=double(I)

end

size=fix(1.5/S)% exp(-1.5^2*pi) <0.1%

%k=2*pi*S^2

%F=S^2/sqrt(2*pi)

k=1

for x=-size:size

for y=-size:size

G(size+x+1,size+y+1)=k*exp(-pi*S^2*(x*x+y*y))*...

(exp(j*(2*pi*F*(x*cos(W)+y*sin(W))+P))-exp(-pi*(F/S)^2+j*P))

end

end

GABOUT=conv2(I,double(G),'same')

数字图像的恢复、增强,乃至复合处理,归根到底只是改善图像的品质,提高图像的可解译性。但处理系统(计算机)并未对图像上地物的类别作出“判决”(解译)。由计算按一定的判别模式来自动完成这一“判决”,便是图像分类处理的过程。

图像分类处理的最终目标是智能化,使遥感图像处理发展成为一种人工智能系统。广义的分类处理,既包括波谱信息的分类,也包括空间信息的分类。后者一般包括图形识别、边缘和线条信息的检测与提取,以及纹理结构分析等,通常也称图像的空间信息分析。关于这一部分对于地质工作者显然感兴趣的内容,可参阅文献[3]等著作。限于篇幅,这里仅介绍按波谱信息分类的基本概念。

(一)图像分类的依据

一般来说,同一类地物有着相似的波谱,在多波段遥感的数字图像中,可以粗略地用它们在各个波段上的像元值的连线(亨利曲线)来表示(图4-29A)由于受光照条件、环境背景等因素的影响,在实际的多维波谱空间中,它们的像元值向量往往不是一个点,而是呈点群分布(集群),不同地物的点群处在不同的位置(图4-29B)不仅如此,在实际图像中,不同地物的波谱集群还存在有交叉过渡,受图像分辨力的限制,一个像汪此元中可能包括有若干个地物类别,即所谓“混合像元”。因此,对不同集群的区分一般要依据它们的统计特征(统计量)。例如,集群位置用均值向量表示、点群的中心及离散度常用标准差或协方差来量度等等数字图像常用的几种统计量见表4-4。

图4-29 索尔顿湖和因佩里亚谷地陆地卫星MSS数字图像上主要几种地物的光谱反射比曲线和集群分布

表4-4 数字图像常用的统计量

图像分类处理的实质就是按概率统计规律,选择适当的判别函数、建立合理的判别模型把这些离散的“集群”分离开来敏粗,并作出判决和归类。通常的做法是,将多维波谱空间划分为若干区域(子空间),位于同一区域内的点归于同一类。子空间划分的标准可以概括为两类:①根据点群的统计特征,确定它所应占据的区域范围。例如,以每一类的均值向量为中心,规定在几个标准差的范围内的点归为一类②确定类别之间的边界,建立边界函数或判别函数。不论采取哪种标准,关键在于确定同一类别在多维波谱空间中的位置(类的均值向量)、范围(协方差矩陈)及类与类边界(判别函数)的确切数值。按确定这些数据是否有已知训练样本(样区)为准,通常把分类技术分为监督和非监督两类。

(二)非监督分类

非监督分类是在没有已知类别的训练数据及分类数的情况下,依据图像数据本身的结构(统计特征)和自然点群分布,按照待分样本在多维波谱空间中亮度值向量的相似程度,由计算机程序自动总结出分类参数,进而逐一对像元作归类,通常也称聚类(集群)分析。使用的方法有图形识别、系统聚类、分裂法和动态聚类等。

其中,比较实用的是动态聚类。它是首先根据经验和分类数,选定若干个均值向量,作为“种子”,建立一批初始中心,进行初步概略的分类,然后根据规定的参数(阈值)检验分类结果,逐步修改调整分类中心,再重新分类,并根据各类离散性统计量(如均方差等)和不同类别之间可分离性统计量(如类间标准化距离等),进行类的合并或分裂此后再修改中心,直至分类结果合理为止。动态聚类中,聚类中心和分类数可以按客观的波谱特征自动调整,分类效果一般比较好,但分类结果的确切含义(类别的属性)需另作分析,从实况调查或已有的地面资料中去确定它们的地物类型。

非监督分类由于事先不需训练样本,故处理速度较快,较客观,并能为监督分类的训练样区选择提供参照,一般在有目的的监督分类之前进行。

(三)监督分类

监督分类一般是先在图像中选取已知样本(训练区)的统计数据,从中找出分类的参数、条件,建立判别函数,然后对整个图像或待分类桥陵镇像元作出判别归类。遥感图像处理中常用的监督分类方法有最小距离法、费歇尔线性判别法、贝叶斯线性和非线性判别法(最大似然法)等。

其中,最小距离法在算法上比较简单:首先在图像显示屏上选出训练样区,并且从图像数据中求出训练样区各个波段的均值和标准差尔后再去计算其它各像元的亮度值向量到训练样区波谱均值向量之间的距离。如果距离小于指定的阈值(一般取标准差的倍数),且与某一类的距离最近,遂将该像元归为某类。该分类法的精度取决于训练样区(地物类别)的多少和样本区的统计精度。由于计算简便,并可按像元顺序逐一扫描归类,一般分类效果也较好,因而是较常用的监督分类方法。

最大似然法也是常用的监督分类方法之一。它是用贝叶斯判别原则进行分析的一种非线性监督分类。简单地说,它可以假定已知的或确定的训练样区典型标准的先验概率,然后把某些特征归纳到某些类型的函数中,根据损失函数的情况,在损失最小时获得最佳判别。该法分类效果较好,但运算量较大。

监督分类的结果明确,分类精度相对较高,但对训练样本的要求较高,因此,使用时须注意应用条件,某一地区建立的判别式对别的地区不一定完全适用。此外,有时训练区并不能完全包括所有的波谱样式,会造成一部分像元找不到归属。故实际工作中,监督分类和非监督分类常常是配合使用,互相补充的。

图像分类处理目前在农林、土地资源遥感调查中应用较广。对于地质体的分类,由于干扰因素较大,不容易取得十分理想的效果,故在地质应用上尚不很普遍。但最近已陆续出现了一批使用分类技术的遥感地质应用成果,较多的是用经变换(比值、K-L等)处理的图像再作分类处理,用于岩性填图或热液蚀变填图等,是值得重视的发展方向。

图像分类是与图像信息提取和增强不同的遥感图像处理中另一重要的埋首方面，与图像增强后仍需人为解译不同，它企图用计算机做出定量的决定来代替人为视觉判译步骤。因此，分类处理后输出的是一幅专题图像。在此图像中，原来图像中的每一个象元依据不同的统计决定准则被划归为不同的地表覆盖类，由于是一种统计决定，必然伴随着某种错误的概率。因此，在逻辑上的合理要求是，对每一个象元所做的决定，应是使整个被分类面积即对大量单个象元的分类的某个错误判据为最小。

以下是几种常用的遥感图像分类方法：

1.最大似然分类（maximum likelihood classification）

最大似然分类是一种基于贝叶斯判别准则的非线性监督分类方法，需要知道已知的或确定的训练样区典型标准的先验概率P（wi）和条件概率密度函数P（wi，x）。P（wi）通常根据各种先验知识给出或假定它们相等：P（wix）则是首先确定其分布形式，然后利用训练样本估计其参数。一般假设为正态分布，或通过数学方法化为正态分布。其判别函数集为：

Di（x）=P（wix），i=1，2，…，m （2-2）

如果Di（x）≥ Dj（x），则x属于wi类。其中，j≠i，j=1，2，…，m。m为类别数。

从上述最大似然分类的说明看，其关键就在于已知类别的定义，先验概率的确定，参与分类的变量的好坏和结果误差评价。直到现在，最大似然分类至少还有两个缺点：一是事先大量人力已知光谱类的选择和定义：二是需要长时间的计算机分类计算时间。实际上这也使得最大似然分类法遥感应用受到了限制，因此许多人专门研究改进算法以便解决和缩减图像分类的时间，提高分类的精度。Solst和Lillesand（1991）为了解决已知类别定义消耗大量人力的缺点，发展了半自动训练法进行已知光谱类的定义。Fabio Maselli等（1992）利用Skidmore和Tumer提出的非参数分类器计算出各已知类训练集的先验概率，然后将它们插入常规的最大似然分类过程中进行分类。该方法融合了非参数和参数分类过程的优点，提高了分类的精度。

通常情况下，地形会影响到训练集数据，这样训练集光谱数据就偏离了最大似然分类的假设条件正态分布，从而常规的最大似然分类法在地形起伏较大的地区效果并不太好。为了解决这一问题，C.Conese和G.Maracchi和F.Maselli（1993）提出了一种改进的最大似然分类算法，即去掉每一类数据集中与第一主成分相关的信息（地形信息）然后再进行分类。通过试验，这种方法是有效的，分类精度得到了提高。

K.Arai（1993）用光谱和空间信息进行分类改进了最大似然分类方法。该方法简单易行，大大提高了正确分类的概率。C.Conese和Fabio Maselli（1992）用误差矩阵提高最大似然分类面积估计的精度。Irina Kerl（1996）加最大似然分类精度的一种方法，即多概率比较法。他对同一遥感数据的原始波段、主成分和植被指数的22种组合进行了最大似然分类，发现没有一种波段组合的分类能给出图像中所有土地利用类型的精确分类，每一波段组合仅对图像中的一两类土地利用类型分类有效。因此他提出将能有效区分出所要决定的土地利用类型的几个波段组合的分类结果进行组合来进行图像分类，并称这种方法为多概率比较法，这种方法的基础就是图像数据不同波段组合的分类结果之间分类概率大小的比较。应用这种方法提高了分类的精度。

2.最小距离分类（minimum distance classification）

最小距离分类是一种线性判别监督分类方法，也需要对训练区模式样本进行统计分析，是大似然分类法中的一种极为重要弯卜数的特殊情况。最小距离分类在算法上比较简单，首先需选出要区分类别的训练样区，并且从图像数据中弊腔求出各类训练样区各个波段的均值和标准差，然后再计算图像中其他各个象元的灰度值向量到各已知类训练样区均值向量之间的距离。如果距离小于指定的阈值（一般取标准差的倍数），且与某一类的距离最近，就将该象元划归为某类。因此称为最小距离分类。该方法的精度主要取决于已知类训练样区的多少和样本区的统计精度。另外，距离度量的方法不同，分类的结果也不相同，常见的有：

（1）明氏距离（minkowski distance）

中亚地区高光谱遥感地物蚀变信息识别与提取

当q=2时，即为欧氏距离，而当q趋于无穷时，得到切比雪夫距离。明氏距离，特别是其中的欧氏距离，在实际中用得较多，但它存在着两方面的缺点：一是它与各指标的量纲有关；为克服这一缺点，常常采用先将数据规格化的方法。二是它没有考虑变量之间的相关性。一种改进的距离就是马氏距离。

（2）马氏距离（mahalanobis distance）

中亚地区高光谱遥感地物蚀变信息识别与提取

当 中各特征间完全不相关，这时的马氏距离即为欧氏距离。

总之，最小距离分类是一个能在程序上经济有效实现的简单方法，与最大似然方法不同，它在理论上并不使平均分类错误为最小，所得到的精度与最大似然分类法可相比拟，而计算时间却只有后者的一半。

3.平行管道分类（parallelepiped classification）

平行管道分类是一种最简单的分类方法，是通过研究训练样区数据的各个光谱成分的直方图来进行分类的图像直方图中灰度值的上下限描述了图像中每个波段中类别的灰度值范围。某一光谱类在所有波段的灰度值范围定义了一个多维的平行管道。通过分类计算，图像中的未知象元被划分到属于已知一光谱类的平行管道内，因此称该方法为平行管道分类。这种方法简单易行，但也有几个缺点：①各已知光谱类的平行管道之间必定具有一定的间隔，如果图像中的象元落在这些间隔内，则这些象元不被分类。②对于图像中相关性强的光谱类，它们所定义的平行管道之间必定具有某些重叠，结果落在重叠区的这些象元不能被很好地分离。③没有考虑已知光谱类的先验概率。④根据直方图定义的已知光谱类的平行管道仅仅是通常用来定义光谱类特征的椭圆平面的粗略表示。

T.W.Kellenberger和K.I.Iten（1996）提出改进方法。改进过的平行管道方法基于每个图像波段内两类之间累积百分比直方图的最大差值自动分离已知目标类，这种分离值和图像波段都是自动选择的。他利用改进过的平行管道分类对TM遥感数据进行了森林覆盖分类，结果取得了较好的效果，提高了分类的精度。

4.模糊分类（fuzzy classification）

模糊理论（fuzzy theory）是处理模糊性的理论的总称，它是以1965年由Zadeh提出的模糊集合论为基础的。模糊分类即是建立在模糊理论之上的分类方法。在常规遥感图像分类中，图像中的每个象元都被划归为某一种地物类型，象元和类别之间一一对应，并没有考虑遥感图像中实际存在的混合象元问题，即某一象元中同时存在有几种地物类型，所以这种分类方法的结果必然存在有误差，这也是常规遥感图像分类方法的局限性。模糊分类法正好克服了这一不足，它利用象元隶属度表示象元的归属问题。当一象元对某一类地物的隶属度为0时，表示该象元不属于该类；当一象元对某一类地物的隶属度为1时，表示该象元属于该类；当一象元对某一类地物的隶属度为0和1之间的数时，表示该象元以该隶属度属于该类，也即表明该象元为混合象元。模糊分类的方法有很多，大多是常规分类方法的扩展。F.Wang（1990）详细讨论了遥感图像模糊监督分类的意义，他认为模糊监督分类与常规分类方法相比，在3个方面有所改进：①用模糊集合表示地学信息更合理；②光谱空间的模糊划分；③分类参数的模糊估计。

Paolo Gamba等（1996）利用完全模糊分类链对多光谱遥感图像进行了分类，他们不仅考虑了图像中象元的光谱特征，而且还考虑了图像中邻近象元之间空间上的相互关系，结果图像的分类精度得到了提高。

5.神经元网络分类（neural network classification）

神经元网络的结构包含一个输入层、一个输出层以及一个或多个隐层。输入层结点数与参加分类的特征数相同，输出层结点数与最终类别数相同，而中间隐含层结点数则由实验来确定。神经元网络分类主要就是网络的训练迭代，要达到一定的精度，往往需要很多次的迭代，这是非常费时的，然而网络训练一经完成，则可较快地应用于分类识别。

F.Amar等（1995）基于分类误差百分比和用户CPU时间比较了4 种神经元网络的遥感应用，他发现后向传播神经元网络算法（back-propagation neural network algorithm）需要最多的迭代次数，是其中最慢获得88%分类精度的算法。动态学习神经元网络算法（dynamic-learning neural network algorithm）只需迭代两次，占用591 s的CPU时间就可获得85%的分类精度，达到86%的分类精度需要迭代5次，占用CPU时间1613.7 s；功能连接神经元网络算（functional link neural network algorithm）分类精度比快速学习神经元网络算（fast-learning neural network algorithm）算法低，但占用了较少的CPU时间；快速学习神经元网络算法提供了最高的分类精度，但两次迭代运算后比功能连接神经元网络算法和动态学习神经元网络算法需要更长的CPU时间。另外，许多研究者利用神经元网络算法对遥感图像数据进行了分类（Mcclellan et al.，1989；Benediktsson et al.，1990；Kanellopoulos et al.，1991；Heermann and Khazenie，1992；Bischof et al.，1992），这些研究表明神经元网络分类能够给出好的分类结果，该结果与常规最大似然分类算法取得的分类精度大致相当，如神经元网络分类算法应用恰当，则分类精度更高。K.Sanjo（1996）研究了神经元网络技术的遥感图像分类，发现监督分类中训练样本数据的排列影响着分类的精度，因此必须注意由于简单地神经元网络算法学习阶段训练样本数据的重新排列所造成的分类结果的改变。另外，他还研究了混合象元问题，结果表明通过学习混合数据，神经元网络算法能有效地减少分类的误差。

6.N维概率密度函数（N-dimensional probability density functions）

Haluk Cetin和Donald W.Levandowski（1991）利用N维概率密度函数对多维遥感数据进行了分类和填图。N维概率密度函数是一种用来显示、分析和分类数据的算法，源于常称的频率透视图，但又克服了早先方法的内在限制。利用N维概率密度函数算法进行的交互式分类过程是一种新的多维数据的分类方法，它提供了遥感数据分布的清晰透视和监督分类中被选择的训练区分布的清晰透视。经过多维数据和训练区分布的N维概率密度函数的制图后，N维概率密度函数空间根据训练区数据的分布被划分，然后将N维概率密度函数的这种划分当作查询表（look-up table）分类遥感图像数据。对非监督分类，N维概率密度函数图可以提供数据分布的一种有价值的代表，数据分布可直接用来选择类别数和数据初始聚类时类均值的位置。与传统的统计分类方法，例如最大似然分类和最小距离分类，需要大量计算机内存、计算缓慢、对显示重叠类能力有限相比，N维概率密度函数分类速度极快，可利用无限制的波段数，图形地显示数据和类别的分布。Haluk Cetin和Timothy A.Warner和Donald W.Levandowski（1993）利用N维概率密度函数对TM等多种遥感数据源数据进行了分类、可视化和增强，结果取得了良好的效果，取得了比传统分类方法较高的分类精度。

7.其他分类方法

N.B.Venkateswarlu和P.S.V.S.K.Raju（1993）提出了一种新的遥感图像快速分类器该分类器是一种非参数分类器，叫做整体平均分类（ensemble-average classifier），利用了最小距离的概念，算法步骤如下：

①计算每一组的整体平均值（一般为均值）Mi，i=1，2，...，C

②如果满足下列两式，那么任一随机象元X将被划归wi组。

XT（Mi-Mj）＜ Tij，j≠ I （2-5）

中亚地区高光谱遥感地物蚀变信息识别与提取

式中Tij=-Tij。

③经过①②步后，随机象元X被划归为正确的类。

另外，通过对参与计算变量的排序和部分一总和逻辑的考虑，可大大降低该算法计算的时间。与最小距离（欧氏距离）和最大似然分类器相比，整体平均分类器所用时间最少，分类精度与最小距离大致相同，对像农田面积和森林这样的名义类型的分类十分有效。

Haluk Cetin（1996）提出了一种分类方法：类间距离频率分布法（interclass distance frequency dis-tribution），这是多光谱数据非参数分类方法的一种。类间距离频率分布过程简单，是一种有力的可视化技术，它图形地显示多光谱数据和类分布。首先选择感兴趣的类，这些类的统计信息从典型的训练样区可获得。利用类的平均测量矢量计算多光谱数据中每个象元的距离，并存放在一个两维数据分布数组中。选择其他类的训练区，训练区数据的分布通过距离计算可获得。通过可视化地检查结果，建立分类查询表（look-up table），然后利用分类查询表进行多光谱图像数据的分类，具体细节请参见原文。

H.N.Srikanta Prakash等（1996）改进了遥感数据凝聚聚类分析，这是一种基于相互近邻概念，用来进行多光谱数据分类的非参数、层次、凝聚聚类分析算法。该方法定义了围绕象元的感兴趣区域（area of interest around each pixel），然后在它内部寻找分类时初始合并 *** 作需要的k最近邻，将象元的特征值、波段值和象元的相对位置值一起考虑，提出了改进的距离量度，这样，大大减少了计算的时间和内存的需求，降低了分类的误差概率。

Steven E.Franklin和Bradley A.Wilson（1992）设计了3阶段分类器进行遥感图像的分类，它由一个基于四叉树的分割算子、一个高斯最小距离均值测试和一个包括辅助地理网数据和光谱曲线测量的最终测试构成。与最大似然分类技术相比，3阶段分类器的总体分类精度得到了提高，减少计算时间，另外仅需最少的训练样区数据（它们在复杂地形区很难获得）。

---