图像特征提取

图像特征提取

摘要

特征提取是计算机视觉和图像处理中的一个概念。它指的是使用计算机提取图像信息，决定每个图像的点是否属于一个图像特征。特征提取的结果是把图像上的点分为不同的子集，这些子集往往属于孤立的点、连续的曲线或者连续的区域。

至今为止特征没有万能和精确的图像特征定义。特征的精确定义往往由问题或者应用类型决定。特征是一个数字图像中“有趣”的部分，它是许多计算机图像分析算法的起点。因此一个算法是否成功往往由它使用和定义的特征决定。因此特征提取最重要的一个特性是“可重复性”：同一场景的不同图像所提取的特征应该是相同的。

特征提取是图象处理中的一个初级运算，也就是说它是对一个图像进行的第一个运算处理。它检查每个像素来确定该像素是否代表一个特征。假如它是一个更大的算法的一部分，那么这个算法一般只检查图像的特征区域。作为特征提取的一个前提运算，输入图像一般通过高斯模糊核在尺度空间中被平滑。此后通过局部导数运算来计算图像的一个或多个特征。

常用的图像特征有颜色特征、纹理特征、形状特征、空间关系特征。当光差图像时，常

常看到的是连续的纹理与灰度级相似的区域，他们相结合形成物体。但如果物体的尺寸很小

或者对比度不高，通常要采用较高的分辨率观察：如果物体的尺寸很大或对比度很强，只需

要降低分辨率。如果物体尺寸有大有小，或对比有强有弱的情况下同事存在，这时提取图像

的特征对进行图像研究有优势。

常用的特征提取方法有：Fourier变换法、窗口Fourier变换（Gabor)、小波变换法、最

小二乘法、边界方向直方图法、基于Tamura纹理特征的纹理特征提取等。

设计内容

课程设计的内容与要求（包括原始数据、技术参数、条件、设计要求等）：一、课程设计的内容

本设计采用边界方向直方图法、基于PCA的图像数据特征提取、基于Tamura纹理特征的纹理特征提取、颜色直方图提取颜色特征等等四种方法设计。

（1）边界方向直方图法

由于单一特征不足以准确地描述图像特征,提出了一种结合颜色特征和边界方向特征的图像检索方法.针对传统颜色直方图中图像对所有像素具有相同重要性的问题进行了改进,提出了像素加权的改进颜色直方图方法;然后采用非分割图像的边界方向直方图方法提取图像的形状特征,该方法相对分割方法具有简单、有效等特点,并对图像的缩放、旋转以及视角具有不变性.为进一步提高图像检索的质量引入相关反馈机制,动态调整两幅图像相似度中颜色特征和方向特征的权值系数,并给出了相应的权值调整算法.实验结果表明,上述方法明显地优于其它方法.小波理论和几个其他课题相关。所有小波变换可以视为时域频域的形式，所以和调和分析相关。所有实际有用的离散小波变换使用包含有限脉冲响应滤波器的滤波器段(filterbank)。构成CWT的小波受海森堡的测不准原理制约，或者说，离散小波基可以在测不准原理的其他形式的上下文中考虑。

通过边缘检测，把图像分为边缘区域和非边缘区域，然后在边缘区域内进行边缘定位．根据局部区域内边缘的直线特性，求得小邻域内直线段的高精度位置；再根据边缘区域内边缘的全局直线特性，用线段的中点来拟合整个直线边缘，得到亚像素精度的图像边缘．在拟合的过程中，根据直线段转角的变化剔除了噪声点，提高了定位精度．并且，根据角度和距离区分出不同直线和它们的交点，给出了图像精确的矢量化结果

图像的边界是指其周围像素灰度有阶跃变化或屋顶变化的那些像素的集合，边界广泛的存在于物体和背

景之间、物体和物体之间，它是图像分割所依赖的重要特征．边界方向直方图具有尺度不变性，能够比较好的

描述图像的大体形状．边界直方图一般是通过边界算子提取边界，得到边界信息后，需要表征这些图像的边

界，对于每一个边界点，根据图像中该点的梯度方向计算出该边界点处法向量的方向角，将空间量化为M级，

计算每个边界点处法向量的方向角落在M级中的频率，这样便得到了边界方向直方图．

图像中像素的梯度向量可以表示为[ ( ，)，)，( ，)，)] ，其中Gx( ，)，)，G ( ，)，)可以用下面的

Sobel算子[8]：

( ，Y)=，( 一1，Y+1)+2×，( ，Y+1)+，( +1，Y+1)一，( 一1，Y一1)一2

×，( ，Y一1)一，( +1，Y一1)

( ，Y)=，( +1，Y一1)+2×，( +1，Y)+，( +1，Y+1)一，( 一1，Y一1)一2

×，( 一1，Y)一，( 一1，Y+1) (4)

I(x，y)表示在图像的(x，y)点像素的亮度。为了减少由于数字化过程中产生的噪声的影响，边界方向直方图是

基于局部像素梯度向量的平均值．因为相反方向的梯度可能会相互抵消，所以并不能直接对局部邻居像素点的梯度向量平均．解决这个问题的一个方法就是在计算平均值前，先对用复数表示的向量进行平方运算，等价于把梯度向量的角度增加一倍．角度增加一倍的相反方向的两个梯度向量分别指向它们的对立梯度向量，从而互相得到增强；而处于正交的两个梯度将会互相抵消，经过平均运算后变回到它们的单一角度的表示。平方的向量可以表示为：

(G +_，G ) =G 一G +2_，G G (5)

梯度平方的平均值可以通过对局部邻居(利用非均匀窗口w)的平均计算得到，

DFx：Σ(G +G；)，DFy：Σ(2GxGr) (6)

W W

现在，平均梯度方向∈[0，180]可以由下面的公式得到：

r l80，DFx=0 and DFx=0

o(x，Y)={[l80 DF ．(7) ——丽，m№M

因此，边界方向直方图可以定义为：

H(k)：H(k)+1，／f ( ，Y)：k and k∈[0，180] (8)

= 180是一个特例．它并不意味着像素的边界方向是水平方向且与=0的相同．它表示在这个像素的周。

围没有亮度变化．边界方向直方图之间距离D 的计算采用与上面颜色直方图距离之间的相同的计算公式(3)而得到。

(2)基于PCA的图像数据特征提取

主成分分析（Principal Component Analysis ，PCA ）是一种掌握事物主要矛盾的统计分析方法，它可以从多元事物中解析出主要影响因素，揭示事物的本质，简化复杂的问题。计算主成分的目的是将高维数据投影到较低维空间。给定n 个变量的m 个观察值，形成一个n ′ m 的数据矩阵，n 通常比较大。对于一个由多个变量描述的复杂事物，人们难以认识，那么是否可以抓住事物主要方面进行重点分析呢？如果事物的主要方面刚好体现在几个主要变量上，我们只需要将这几个变量分离出来，进行详细分析。但是，在一般情况下，并不能直接找出这样的关键变量。这时我们可以用原有变量的线性组合来表示事物的主要方面，PCA 就是这样一种分析方法。

图像的有效检索手段逐渐引起人们的重视,传统的采用基于关键字或描述性文本的数据库检索方式,已远远不能满足人们的需要,近年来,国际上广泛开展了基于内容的图像检索(CBIR) 的研究,CBIR 是指根据图像的颜色、形状、纹理等特征以及这些特征的组合来查询图像,是计算机图像处理和数据库技术的有效结合。颜色是彩色图像的最显著特征之一。因此,基于颜色的查询是基于内容的图像检索中最基本的方法。这种查询可针对任何类型的彩色图像。目前,人们已提出了许多种借助于颜色特征对图像进行检索的方法。这些方法常用的彩色空间有RGB 和HSV ,提取的主要特征是颜色直方图。常用的计算直方图之间距离的

方法有直方图相交法、距离法等，以直方图作为图像的颜色特征进行图像检索,需存储大量的冗余特征信息,不能精确而简洁地描述图像。本文提出了一种基于HSV 颜色模型,利用主元分析法(PCA) 对图像的颜色进行特征提取的方法,大大降低了特征向量的维数,减小了特征数据库的规模。

图像的颜色有多种表示方式,其中HSV 颜色模型是一种适合肉眼分辨的模型,它把彩色信息表示为三种属性: 色调h 、饱和度s 和亮度v 。HSV 模型的色调h 是由颜色名称来辨别的,如红、橙、绿,它用角度0°～360°度量;亮度v 是颜色的明暗程度, 通常用百分比度量, 从黑0 到白100 %;饱和度s 指颜色的深浅,用百分比来度量,为从0 到完全饱和的100 %·这种颜色模型用Munsell 三维空间坐标系统表示,因坐标之间的心理感知独立性,因此,可以独立感知各颜色分量的变化;且这种颜色模型具有线性伸缩性,可感知的颜色差是与颜色分量的相应样值上的欧氏距离成比例的·在CBIR 中应用这种模型更适合用户的

主元分析法(PCA) 对L 的降维处理由节1. 2 可得到L 矩阵, L 为一m ×n 矩阵, m 、n 由图像尺寸决定。主元分析法( PCA) 是一种线性降维技术,其基本思想是通过对数据协方差矩阵的分解,在其n 个特征值中取前a 个特征值, ( n - a) 个特征值被滤出。a 远小于n ,而a个特征值对应的特征向量构成负荷矩阵P ∈R n ×a , L 到低维空间的投影就包含在得分矩阵中T = L P (3)由T 返回到n 维空间的投影为L = TPT (4)由^L 张成的子空间称为得分空间,其所包含的信息量近似原空间L 中所包含的信息量[9 ]。由矩阵L 可求出其协方差矩阵S 的特征值分解S = 1/ ( n - 1) L T L = VΛV T (5)Λ包含幅值递减的非负实特征值(λ1 ≥λ2 ≥?≥λn ≥0) ·为了最优地获取数据的变化量,这里a =6 ,即取前6 个特征值。

利用主元分析法(PCA) 对彩色图像颜色特征进行提取,可将原图像从m ×n ×3 降低至m ×a ( a 远小于n) ,大大降低了特征向量的维数,减小了特征数据库的规模,为图像的检索建立了良好的基础·将提取的颜色特征与纹理特征相结合作为支持向量机(SVM) 的输入,同时引入相关反馈来对图像进行分类与检索,以便近一步提高检索效率是有待进一步研究的内容。

(3)基于Tamura纹理特征的纹理特征提取

随着多媒体信息的广泛应用,并在数据库系统和计算机视觉两大研究领域的共同推动下,图像检索技术己逐渐成为一个非常活跃的研究领域。如何有效地对这些图像进行分析、存储和检索是一个亟待解决的问题。基于内容的图像检索技术能有效的解决这一问题,成为研究的热点。基于内容的图像检索(CBIR)不同于传统的基于文本的图像检索,它实际上是一种模糊查询技术,通过对图像提取一定的特征,找出在特征空间中与查询要求接近的图像,从而实现在图像数据库中自动地、智能地检索、查询和管理图像。在CBIR系统中,特征的提取和匹配算法是决定图像检索结果的关键。在已经存在的颜色特征,纹理特征,形状特征等几种特征提取方法中,由于纹理特征能够描述图像的平滑,稀疏,规则性等特性,因此本文决定采用纹理特征作为图像检索的依据。首先,介绍了国内外基于内容的图像检索系统的研究现状,从基于内容的图像检索系统的基本原理入手,论述了基于内容的图像检索系统的通用框架、关键技术、检索特点以及应用领域。其次,描述了纹理和纹理特征,介绍了几种常用的纹理分析方法,有灰度直方图、边缘方向直方图、空间灰度共生矩阵、Tamura纹理特征、Gabor滤波器和Gabor小波以及图像纹理谱分析方法,本文在传统的纹理特征提取方法的基础上,利用Gabor小波变换在数字图像中局部区域的频率,方向信息的优异性能,和共生矩阵对图像的整体区域有着较好的处理效果的特点,尝试了结合灰度共生矩阵和Gabor滤波器来提取纹理特

征的方法,并对提取的特征进行高斯归一化处理。另外,在图像相似性度量方面,阐述了几何矩阵模型和集合理论模型,给出了两种多特征组合的相似性度量结构,以及相似性度量的方法。同时介绍了在基于内容的图像检索技术中相关性反馈技术。最后,本文设计了一个基于纹理特征的图像检索原型系统,并对系统的各个模块及其功能进行了介绍,对本文提出的检索技术进行了仿真实验,给出了实验结果,并用排序平价法和平均检索率对算法进行了评价。

Tamura纹理特征表示方法,并对其相似性度量方法进行了改进,使用直方图表示了最重

要的粗糙度、对比度、方向度特征。本文还将粗糙度直方图改进为能够适应图片尺寸变化的方式,另外使用对数非均匀量化方式将对比度特征形成直方图,该处理方法能够更适应人眼的认知特点。同时,对每种纹理特征都给出了相似性度量仿真实例,证实该方法具有较好效果。最后使用三种特征的加权距离对图像的纹理特征相似性进行度量,并给出了图像检索的结果和性能评价。3)兴趣点特征提取:本文综合了Harris兴趣点标定方法和区域矩不变特性对图像兴趣点特征进行表达。该特征能够适应光照、视角、缩放等场景的变化,文中给出图片实例对这一特点进行了验证。本文根据实验提出了图像的兴趣点相似性度量方法,并给出了图像检索的结果和性能评价。本文第二部分研究了区域相似性评价和图片共性特征归纳方法。区域分割是该过程的基础,本文介绍了一种简单的基于颜色和纹理的区域分割法,并提出了图像空间结构相似性度量方法。在对具有相似内容的图片进行共性特征提取过程中,根据不同的应用场景选取了不同的特征子集以降低复杂度。本文还通过实例,阐明了共性特征提取方法、步骤及概念表达的方法,最后在给出的共性特征标示图中验证了该方法的效果和可行性。

(4)颜色直方图提取颜色特征

颜色直方图是在许多图像检索系统中被广泛采用的颜色特征，具有特征提取和相似计算简便的优点。Swain 和Ballard 首先把直方图应用于图像检索[1],Smith J.R.提出了在HSV （Hue，Saturation，Value）空间的166 色量化方法[2].但在对彩色图像进行检索时，颜色直方图还存在两个问题：（1）计算颜色直方图需要将颜色空间划分成若干个小的颜色区间，当区间的数目很大的时候不但会增加计算负担，也不利于在图像库中建立索引；（2）颜色直方图丢弃了图像的空间信息，对于两幅颜色分布相同，内容并不一样的图像缺乏区分度。

首先，根据HSV 颜色空间的色彩（Hue）、饱和度（Saturation）和亮度（Value）的值将图像的颜色划分为9 个区间，计算出每一种颜色区间所占的比重作为直方图信息。然后计算出每一颜色区间所包含象素的空间分布信息。在图像检索时，综合考虑颜色直方图和空间信息的相似度。

颜色直方图特征，HSV 颜色空间有三个属性：色彩（Hue），又称为色调，是指色彩的相貌，与色彩明暗无关；饱和度（Saturation），是色彩纯与不纯的分别；亮度（Value），是指色彩的明暗程度，光度的高低，越接近白色亮度越高，越接近灰色或黑色亮度越低。根据对HSV 空间的研究，人们能够分辨出来的颜色大体上包括：红、黄、绿、青、蓝、紫，他们在色彩上的分布并非是均匀的，于是可以把色彩分为不等间隔的六份 {Red，Yellow，Green，Cyan，Blue，Magenta}，而亮度则可以划分为不等间隔的黑灰白三份{Black，Gray，White}。这样，HSV 颜色空间就划分为9 个不同的空间区域。

对于图像中的一个象素，可以根据其Saturation 值的判定函数阈值将其划分到相应颜

色区间。[3]中采用了一个线性函数 1-0.8*Value 作为阈值，但是HSV 是一个不均匀的颜色空间，线性判别函数不能准确地划分彩色和非彩色的子空间。可以采用一种更有效的非线性

函数：

作为阈值。

定义阈值：

对于图像中的任一象素X(H,S,V)，判断其属于哪一个颜色区间的公式由表1 所示：

统计出每种颜色在图像中占的比重就可以得到图像的颜色直方图信息。设图像的颜色直

方图特征为His, 可由式（1）表示：

颜色直方图没有表达出颜色空间分布的信息，也没有反映出图像中物体的形状。因此，可以结合颜色区间的空间分布信息。颜色量化的结果也把图像分为了9 个颜色区间，可以计算各个颜色区间所含象素的分布状态作为其空间信息。

设图像颜色区间C 的空间信息为h(c)，它可以由公式（2）计算：

颜色直方图相似性计算：

常用的基于颜色直方图的相似性度量方法包括直方图相交、绝对值距离、欧式距离、加权距离等。本文采用直方图相交作为图像间的相似性度量标准。直方图的相交是指两个直方图在每个颜色区间中共有的象素数量，度量方法如下：

给定查询图像Q 和图像库中的图像T，提取的图像颜色直方图特征分别为：[ q1 c , …, qN c c ]和[ t1 c , …,tN c c ]，它们之间的相交距离由式（4）表示：

空间信息的相似性计算：

UFM[4]模糊匹配方法可以用来计算图像的空间信息相似度，定义图像颜色区间的平均距

离为f d ，计算f d 的公式是（6）：

其中， c N 为颜色区域数目。

对于输入图像q P 和目标图像t P ，它们的空间信息相似向量设为(P q ,P t ) M ,其计算公式：

其中，i h 和j h 分别为输入图像和目标图像的空间信息。

设图像颜色区间在图像中的大小权重向量为w 。有多种选择权重的方法，一种方法是每一个区间都赋予相同的权重；另一种方法是以区间的面积在整个图像中的比重作为其权重，即面积越大的区间越重要；还有一种方法是对位于图像边缘或者中心部分的区间赋予较大的权重。本文采用第二种方法。

输入图像q P 和目标图像t P 的空间信息相似度为k δ，k δ由公式（8）计算。

综合颜色直方图和空间信息的图像相似度：

颜色直方图虽然有计算简单、对平移和旋转不敏感的优点，但无法捕捉颜色组成之间的

空间关系，这无疑是不够准确的。颜色直方图量化本身就是一种粗略的划分方法，单独采用基于颜色直方图量化划分的空间信息，检索结果也不理想。本文采用综合两种特征的方法，对颜色直方图和空间信息的相似度以不同的权重累加，检索效果更加理想。

设输入图像q P 和目标图像t P 的相似度为δ，δ由公式（9）计算：

二、程序源代码：

2.1、边界方向直方图法提取形状特征

originailimage=imread('1.jpg');

image_I=rgb2gray(originailimage);

% [image_x,image_y]=size(image_I);

% if image_x

% image_I=image_I';

figure(1),subplot(121),title('灰度图'),

imshow(image_I);

% 找到边缘

I=image_I;

image_bw=edge(I,'canny');

subplot(122),title('边界检测图')

imshow(image_bw);

[L,NUM]=bwlabeln(image_bw,8); % 连通性

%~~~~~~~将图像边缘量化为101份(1-100 101) 画出边缘直方图，量化份数因情况不同而定

num_line=zeros(1,101);

for i=1:NUM

number=length(find(L==i));

for j=1:100

if number==j

num_line(j)=num_line(j)+1;

end

end

if number>100

num_line(101)=num_line(101)+1;

end

end

EH_1=[];

for i=1:101

EH_1=[EH_1 num_line(i)];

end

x1=1:101;

y1=EH_1;

figure(2),title('plot曲线直方图'),

plot(x1,y1);

%~~~~亦可将图像量化为13份（1-10 11 12 13）画出bar投影图像num=zeros(1,13);

for i=1:NUM

number=length(find(L==i));

for j=1:10

if number==j

num(j)=num(j)+1;

end

end

if number>10 & number<=20

num(11)=num(11)+1;

end

if number>20 & number<=100

num(12)=num(12)+1;

end

if number>100

num(13)=num(13)+1;

end

end

EH=[];

for i=1:13

EH=[EH num(i)];

end

x=1:13;

y=EH;

figure(3),title('bar投影直方图'),

bar(x,y);

2.2、基于PCA的图像数据特征提取

clear all

% Define database

x=imread('3.jpg');

figure(1)

plot(x(:,1),x(:,2),'ko')

xlabel('Valve Position (%)')

ylabel('Fermenter Temperature (C)')

% Calculate the covariance of dataset x (relationship between Valve and

% Temperature)

covariance=cov(x);

% Singular Value Decomposition x=T*E*P

% E is the eignvalue of x; U is orthnormal; V=U'.

%[U E V]=svd(covariance);

%U=U'; V=V';

[U E] = eig(covariance);

E = diag([E(2,2),E(1,1)]); % arrange the Eignvalue in a decending order

U = [U(:,2),U(:,1)]; % arrange the Eignvector to make it compatable with Eignvalue.

% Mean center the raw data

Valvebar=x(:,1)-mean(x(:,1));

Tempbar=x(:,2)-mean(x(:,2));

xbar=[Valvebar,Tempbar];

% Principal axis rotation of the covariance matrix.

z=U'*xbar';

z=z';

%covz is the covariance matrix of principal components which is equal to E

%matrix. The variance of transformed variable (z1 & z2) will have variance

%0.2283 and 0.0123 respectively.

covz=U'*covariance*U;

% The first column of U factor is - and + with nearly equal value which % means the first principal component is related to variability which both

% measurements have difference. The second column is both positive

which

% means the 2nd PC is concentrate on the varaiability of the common between

% the two.

% scaling of PCS.

for i=1:2

Vs(:,i)=sqrt(E(i,i))*U(:,i);

Ws(:,i)=U(:,i)/sqrt(E(i,i));

end

% after rescale: Vs'*Vs=E; V'*COV*Vs=E^2; Ws'*Ws=inv(E); Ws'*COV*Ws=I y=Ws'*xbar';

%define T2

T2=diag(y'*y);

% find the control limit of T2

p=2; %PC number

n=10; % sample number

T2Upper=p*(n-1)/(n-p)*finv(0.95,p,n-p);

% draw the T2 control charts

figure(2)

plot(T2,'ko')

hold on

plot(T2Upper*ones(1,13),'k')

xlabel('Sample Number')

ylabel('T2 Score')

hold off

% Draw the control ellipse

%

s1^2*s2^2/(s1^2*s2^2-s12^2)*[(x1-x1mean)^2/s1^2+(x2-x2mean)^2/s2^

2-2*s12(

% x1-x1mean)(x2-x2mean)/s1^2/s2^2]=T2^2(upper control limit)

% S=covariance;

% a=S(1)*S(4)/(S(1)*S(4)-S(2)^2);

%

b=xbar(:,1).^2/S(1)+xbar(:,2).^2/S(4)-2*S(2)*xbar(:,1).*xbar(:,2) ./S(1)/S(4);

% error detection

% two new points is added to be detected which are [10.1,25.2];[9.9;24];

% plot two values with the original 10.

Valve2=[10.2,9.9];Temp2=[25.1,25.5];

figure(1)

hold on

plot(Valve2,Temp2,'k*')

hold off

% plot the control charts for valve and temperature separately. figure(3)

hold on

plot(1:10,Valve,'ko')

plot(11:12,Valve2,'k*')

plot(0:13, mean(Valve)*ones(1,14),'k')

plot(0:13, (mean(Valve)-std(Valve)*2)*ones(1,14),'k')

plot(0:13, (mean(Valve)+std(Valve)*2)*ones(1,14),'k')

xlabel('Sample Number')

ylabel('Valve Position (%)')

hold off

figure(4)

hold on

plot(1:10,Temp,'ko')

plot(11:12,Temp2,'k*')

plot(0:13, mean(Temp)*ones(1,14),'k')

plot(0:13, (mean(Temp)+std(Temp)*2)*ones(1,14),'k')

plot(0:13, (mean(Temp)-std(Temp)*2)*ones(1,14),'k')

xlabel('Sample Number')

ylabel('Temperature (C)')

hold off

% transform new observations to PC axis

y2=Ws'*([Valve2-mean(Valve);Temp2-mean(Temp)]);

% plot the control charts for principal components

figure(5)

hold on

plot(1:10,y(1,:),'ko')

plot(11:12,y2(1,:),'k*')

plot(0:13, mean(y(1,:))*ones(1,14),'k')

plot(0:13, -2*ones(1,14),'k')

plot(0:13, 2*ones(1,14),'k')

xlabel ('Sample Number')

ylabel ('Principal Component 1')

hold off

figure(6)

hold on

plot(1:10,y(2,:),'ko')

plot(11:12,y2(2,:),'k*')

plot(0:13,mean(y(2,:))*ones(1,14),'k')

plot(0:13, -2*ones(1,14),'k')

plot(0:13, 2*ones(1,14),'k')

xlabel ('Sample Number')

ylabel ('Principal Component 2')

hold off

% plot the T2 chart for the two new points figure(2)

T22=diag(y2'*y2);

hold on

plot([11,12],T22,'k*')

hold off

% draw the control ellipse

ra = sqrt(T2Upper*E(1,1));

rb = sqrt(T2Upper*E(2,2));

V=[sqrt(E(1,1))*U(:,1),sqrt(E(2,2))*-U(:,2)]; k1 = V(2,1)/V(1,1);

k2 = V(2,2)/V(1,2);

ang = atan(k1); %atan(k1);

x0 = mean(x(:,1));

y0 = mean(x(:,2));

% i = 9.5:0.01:10.5;

% s = 9.72:0.01:10.3;

% j = k1*i + 42.6303;

% l = k2*s + 19.3279;

figure (7)

h=ellipse(ra,rb,ang,x0,y0,'k',1000);

hold on

plot(x(:,1),x(:,2),'ko')

plot(Valve2, Temp2,'k*')

% plot(i,j)

% plot(s,l)

xlabel ('Valve Position')

ylabel ('Temperature')

hold off

2.3、基于Tamura纹理特征的纹理特征提取

%计算Tamura纹理特征

close all

clear all

clc

t0=cputime;

I = imread('F-18.bmp');

[Nx,Ny] = size(I);

Ng=256;

G=double(I);

%计算粗糙度(coarseness)

Sbest=zeros(Nx,Ny);

E0h=zeros(Nx,Ny);

E0v=zeros(Nx,Ny);

E1h=zeros(Nx,Ny);

E1v=zeros(Nx,Ny);

E2h=zeros(Nx,Ny);

E2v=zeros(Nx,Ny);

E3h=zeros(Nx,Ny);

E3v=zeros(Nx,Ny);

E4h=zeros(Nx,Ny);

E4v=zeros(Nx,Ny);

E5h=zeros(Nx,Ny);

E5v=zeros(Nx,Ny);

flag=0;

for i=1:Nx

for j=2:Ny

E0h(i,j)=G(i,j)-G(i,j-1); end

end

E0h=E0h/2;

for i=1:Nx-1

for j=1:Ny

E0v(i,j)=G(i,j)-G(i+1,j); end

end

E0v=E0v/2;

%图片大小必须大于4*4才能计算E1h、E1v if (Nx<4||Ny<4)

flag=1;

end

if(flag==0)

for j=3:Ny-1

E1h(i,j)=sum(sum(G(i:i+1,j:j+1)))-sum(sum(G(i:i+1,j-2:j-1)));

end

end

for i=2:Nx-2

for j=2:Ny

E1v(i,j)=sum(sum(G(i-1:i,j-1:j)))-sum(sum(G(i+1:i+2,j-1:j)));

end

end

E1h=E1h/4;

E1v=E1v/4;

end

%图片大小必须大于8*8才能计算E2h、E2v

if (Nx<8||Ny<8)

flag=1;

end

if(flag==0)

for i=2:Nx-2

for j=5:Ny-3

E2h(i,j)=sum(sum(G(i-1:i+2,j:j+3)))-sum(sum(G(i-1:i+2,j-4:j-1))); end

end

for i=4:Nx-4

for j=3:Ny-1

E2v(i,j)=sum(sum(G(i-3:i,j-2:j+1)))-sum(sum(G(i+1:i+4,j-2:j+1))); end

end

E2h=E2h/16;

E2v=E2v/16;

end

%图片大小必须大于16*16才能计算E3h、E3v

if (Nx<16||Ny<16)

flag=1

end

if(flag==0)

for j=9:Ny-7

E3h(i,j)=sum(sum(G(i-3:i+4,j:j+7)))-sum(sum(G(i-3:i+4,j-8:j-1))); end

end

for i=8:Nx-8

for j=5:Ny-3

E3v(i,j)=sum(sum(G(i-7:i,j-4:j+3)))-sum(sum(G(i+1:i+8,j-4:j+3))); end

end

E3h=E3h/64;

E3v=E3v/64;

end

%图片大小必须大于32*32才能计算E4h、E4v

if (Nx<32||Ny<32)

flag=1;

end

if(flag==0)

for i=8:Nx-8

for j=17:Ny-15

E4h(i,j)=sum(sum(G(i-7:i+8,j:j+15)))-sum(sum(G(i-7:i+8,j-16:j-1)) );

end

end

for i=16:Nx-16

for j=9:Ny-7

E4v(i,j)=sum(sum(G(i-15:i,j-8:j+7)))-sum(sum(G(i+1:i+16,j-8:j+7)) );

end

end

E4h=E4h/256;

E4v=E4v/256;

end

%图片大小必须大于64*64才能计算E5h、E5v

if (Nx<64||Ny<64)

flag=1;

end

if(flag==0)

for i=16:Nx-16

for j=33:Ny-31

E5h(i,j)=sum(sum(G(i-15:i+16,j:j+31)))-sum(sum(G(i-15:i+16,j-32:j -31)));

end

end

for i=32:Nx-32

for j=17:Ny-15

E5v(i,j)=sum(sum(G(i-31:i,j-16:j+15)))-sum(sum(G(i+1:i+32,j-16:j+ 15)));

end

end

E5h=E5h/1024;

E5v=E5v/1024;

end

for i=1:Nx

for j=1:Ny

[maxv,index]=max([E0h(i,j),E0v(i,j),E1h(i,j),E1v(i,j),E2h(i,j),E2 v(i,j),E3h(i,j),E3v(i,j),E4h(i,j),E4v(i,j),E5h(i,j),E5v(i,j)]); k=floor((index+1)/2);

Sbest(i,j)=2.^k;

end

end

Fcoarseness=sum(sum(Sbest))/(Nx*Ny);

%计算对比度

[counts,graylevels]=imhist(I);

PI=counts/(Nx*Ny);

averagevalue=sum(graylevels.*PI);

u4=sum((graylevels-repmat(averagevalue,[256,1])).^4.*PI); standarddeviation=sum((graylevels-repmat(averagevalue,[256,1])).^ 2.*PI);

alpha4=u4/standarddeviation^2;

Fcontrast=sqrt(standarddeviation)/alpha4.^(1/4);

%计算方向度

PrewittH=[-1 0 1;-1 0 1;-1 0 1];

PrewittV=[1 1 1;0 0 0;-1 -1 -1];

%计算横向梯度

deltaH=zeros(Nx,Ny);

for i=2:Nx-1

for j=2:Ny-1

deltaH(i,j)=sum(sum(G(i-1:i+1,j-1:j+1).*PrewittH)); end

end

for j=2:Ny-1

deltaH(1,j)=G(1,j+1)-G(1,j);

deltaH(Nx,j)=G(Nx,j+1)-G(Nx,j);

end

for i=1:Nx

deltaH(i,1)=G(i,2)-G(i,1);

deltaH(i,Ny)=G(i,Ny)-G(i,Ny-1);

end

%计算竖向梯度

deltaV=zeros(Nx,Ny);

for i=2:Nx-1

for j=2:Ny-1

deltaV(i,j)=sum(sum(G(i-1:i+1,j-1:j+1).*PrewittV)); end

end

for j=1:Ny

deltaV(1,j)=G(2,j)-G(1,j);

deltaV(Nx,j)=G(Nx,j)-G(Nx-1,j);

end

for i=2:Nx-1

deltaV(i,1)=G(i+1,1)-G(i,1);

deltaV(i,Ny)=G(i+1,Ny)-G(i,Ny);

end

%梯度向量模

deltaG=(abs(deltaH)+abs(deltaV))/2;

%梯度向量方向

theta=zeros(Nx,Ny);

for i=1:Nx

for j=1:Ny

if (deltaH(i,j)==0)&&(deltaV(i,j)==0)

elseif deltaH(i,j)==0

theta(i,j)=pi;

else

theta(i,j)=atan(deltaV(i,j)/deltaH(i,j))+pi/2; end

end

end

theta1=reshape(theta,1,[]);

phai=0:0.0001:pi;

HD1=hist(theta1,phai);

HD1=HD1/(Nx*Ny);

HD2=zeros(size(HD1));

%定义一个阈值THRESHOLD

THRESHOLD=0.01;

for m=1:length(HD2)

if HD1(m)>=THRESHOLD

HD2(m)=HD1(m);

end

end

[c,index]=max(HD2);

phaiP=index*0.0001;

Fdirection=0;

for m=1:length(HD2)

if HD2(m)~=0

Fdirection=Fdirection+(phai(m)-phaiP)^2*HD2(m);

end

end

disp('粗糙度：');display(Fcoarseness)

disp('对比度：');display(Fcontrast)

disp('方向度：');display(Fdirection)

deltaT=cputime-t0;

display(deltaT);

2.4、颜色直方图提取颜色特征

originailimage=imread('1.jpg');

image_I=rgb2gray(originailimage);

% [image_x,image_y]=size(image_I);

% if image_x

% image_I=image_I';

figure(1),subplot(121),title('灰度图'),

imshow(image_I);

% 找到边缘

I=image_I;

image_bw=edge(I,'canny');

subplot(122),title('边界检测图')

imshow(image_bw);

[L,NUM]=bwlabeln(image_bw,8); % 连通性

%~~~~~~~将图像边缘量化为101份(1-100 101) 画出边缘直方图，量化份数因情况

不同而定

num_line=zeros(1,101);

for i=1:NUM

number=length(find(L==i));

for j=1:100

if number==j

num_line(j)=num_line(j)+1;

end

end

if number>100

num_line(101)=num_line(101)+1;

end

end

EH_1=[];

for i=1:101

EH_1=[EH_1 num_line(i)];

end

x1=1:101;

y1=EH_1;

figure(2),title('plot曲线直方图'),

plot(x1,y1);

%~~~~亦可将图像量化为13份（1-10 11 12 13）画出bar投影图像 num=zeros(1,13);

for i=1:NUM

number=length(find(L==i));

for j=1:10

if number==j

num(j)=num(j)+1;

end

end

if number>10 & number<=20

num(11)=num(11)+1;

end

if number>20 & number<=100

num(12)=num(12)+1;

end

if number>100

num(13)=num(13)+1;

end

end

EH=[];

for i=1:13

EH=[EH num(i)];