基于快速傅里叶变换的形状检索

摘要：近年来，随着多媒体技术和数字设备的出现，如何有效地管理和访问图像信息已成为人们亟待解决的问题。因此，一种新的图像检索技术——基于内容的图像检索技术被提出来。本文在已有研究成果的基础上，提出了自己的基于形状特征的图像检索方法：基于二维快速傅里叶变换的形状检索算法。

Abstract： In recent years， with the emergence of multimedia technology and digital equipment， how to effectively manage and access information has become an urgent issue. Therefore， a new content-based image retrieval technology was raised. This article， on the basis of the research results， presented the approach to image retrieval based on shape feature： shape retrieval based on two-dimensional Fast Fourier transform algorithm.

关键词：图像检索；形状；傅里叶变换

Key words： image retrieval；shape；Fourier transform

中图分类号：TP39 文献标识码：A 文章编号：1006-4311（2012）32-0198-05

0 引言

本文旨在研究基于形状区域特征的图像检索，提出了一种基于二维快速傅里叶变换的的形状检索算法。第一节介绍了的相关知识，第二节介绍了基于快速傅里叶变换的频域特征提取的具体实现技术，第三节介绍了特征匹配的方法，最后是实验的结果与分析。

1 傅里叶变换

本文算法是在快速傅里叶变换的理论基础上提出的，所以本文首先对傅里叶变换进行了研究。

1.1 傅里叶变换的基本概念 傅里叶变换在数字图像处理中有着广泛的应用，其在数学中的定义非常严格，定义如下：设f（x）为x的函数，若f（x）满足狄里赫拉条件：①具有有限个间断点；②具有有限个极点；③绝对可积。

但图像必须在空间和灰度上都离散化才能被计算机处理，要对数字图像进行傅里叶变换，必须引入离散傅里叶变换（DTF，Discrete Fourier Transform）的概念。

1.1.1 一维离散傅里叶变换 对一个连续f（x）等间隔采样，可得到一个离散序列。设共采了N个样，则这个离散序列可表示为{f（0），f（1），f（2），…，f（N-1）}，则其离散傅里叶变换F（u）为：

F（u）=I{f（x）}=■■f（x）e■■u=0，1，…，N-1（1）

离散傅里叶反变换为：

f（x）=I■{F（u）}=■F（u）e■x=0，1，…，N-1（2）

根据欧拉公式e±ix=cosx±isinx，式（1）可写成：

exp（-j2？仔xu）=cos2？仔ux-jsin2？仔ux（3）

所以，一维离散傅里叶函数F（u）也可写成复数函数形式：F（u）=R（u）+jI（u）（4）

其中R（u）和I（u）分别为F（u）的实部和虚部。

其中F（u）=R■（u）+I■（u）■（5）

上式称为傅里叶变换的幅度或频率谱。

其中？准（u）=arctan[I（u）/R（u）]（6）

上式称为变换的相角或相位谱。

其中P（u）=F（u）■■=R■（u）+I■（u）（7）

上式称为功率谱。

1.1.2 二维离散傅里叶变换 在数字图像中，像素点是二维离散的，设一个图像尺寸为M×N（M为图像的高，N为图像的长）的函数为f（x，y），则图像可表示成为一个二维的离散点序列，如式（8）表示：

f（x，y）=■

（8）

对于图像的像素点f（x，y），其二维离散傅里叶变换为：

F（u，v）=I{f（x，y）}=■■■f（x，y）e■

（9）

其中：u=0，1，…，M-1；v=0，1，…，N-1

其傅里叶反变换为：

f（x，y）=I■{F（u，v）}=■■■F（u，v）e■

（10）

其中：x=0，1，…，M-1；y=0，1，…，N-1

考虑到在数字图像处理中，图像取样一般是方阵，即M=N，则二维离散傅里叶变换公式为：

F（u，v）=I{f（x，y）}=■■■F（x，y）e■（11）

其中：u，v=0，1，…，N-1

f（x，y）=I■{F（u，v）}=■■■F（u，v）e■（12）

其中：x，y=0，1，…，N-1

与在一维时的情况类似，二维傅里叶变换的频谱、相位角、功率谱如下：

F（u，v）=[R■（u，v）+I■（u，v）]■（13）

？准（u，v）=arctan[I（u，v）/R（u，v）]（14）

P（u，v）=F（u，v）■=R■（u，v）+I■（u，v）（15）

1.2 傅里叶变换的性质 由于我们主要研究的是二维数字图像，而数字化图像都是离散化了的图像，因此下面对二维离散傅里叶变换的一些重要性质进行简述。

①线性。离散傅里叶变换是一个线性变换。

②分离性。一个二维离散傅里叶变换也可以用二次一维离散傅里叶变换来实现。

③对称性。如果离散函数f（x，y）的离散傅里叶变换为F（u，v），那么：I[F（x，y）]=MN·f（-u，-v）（16）

④周期性。离散傅里叶变换和傅里叶反变换均是以N为周期。

⑤共轭性。如果离散函数f（x，y）的傅里叶变换为F（u，v），F*（-u，-v）为f（-x，-y）离散傅里叶变换的共轭函数，则它的傅里叶变换具有共轭对称性。其共轭性说明离散函数f（x，y）经过正交变换后得到的F（u，v）是以原点为中心对称的。利用该特性，只要求出半个周期内的值就可以得到整个周期的值。

⑥平移特性。离散傅里叶变换对的平移特性可写成：

I[f（x，y）e■]=F（u-u0，v-v0）（17）

I[f（x-x0，y-y0）]=F（u，v）e■（18）

式（17）表示将f（x，y）与一个指数相乘就相当于把其变换后的频域中心移动到新的位置。式（18）表明将F（u，v）与一个指数相乘就相当于把其反变换后的空域中心移动到新的位置。

⑦旋转特性。如果空间域离散函数旋转角度？兹0，则在变换域中该离散函数的离散傅里叶变换函数也将旋转同样的角度。

⑧尺度变换特性。如果函数f（x，y）的傅里叶变换为F（u，v），给定两个标量a和b，则：

I[af（x，y）]=aF（u，v）（19）

I[f（ax，by）]=■F（■，■）（20）

⑨平均值。一幅图像的灰度平均值可由DFT在原点处的值求得，即：

■（x，y）=■F（0，0）（21）

⑩卷积定理。设f（x，y）和g（x，y）分别用尺寸A×B和C×D，周期为M和N的离散数组表示。为了防止卷积后产生重迭误差，需要选择M？叟A+C-1和N？叟B+D-1。

则二维离散傅里叶卷积定理如下：

I[fe（x，y）*ge（x，y）]=Fe（u，v）·Ge（u，v）（22）

I[fe（x，y）·ge（x，y）]=Fe（u，v）？鄢Ge（u，v）（23）

其中，Fe（u，v）为离散函数fe（x，y）的离散傅里叶变换，Ge（u，v）为离散函数ge（x，y）的离散傅里叶变换。

{11}相关定理。对于离散函数f（x，y）和g（x，y）的相关运算？莓，同样必须先对f（x，y）和g（x，y）进行扩展处理，然后再定义相关运算？莓如下：

fe（x，y）？莓ge（x，y）=■■■fe（m，n）ge（x+m，y+n） ■（24）

则相关定理如下：

I[fe（x，y）？莓ge（x，y）]=Fe（u，v）·G■■（u，v）（25）

I[fe（x，y）·g■■（x，y）]=Fe（u，v）？莓G■（u，v）（26）

其中Fe（u，v）为函数fe（x，y）的傅里叶变换，Ge（u，v）为函数ge（x，y）的傅里叶变换；G■■（u，v）为G■■（u，v）的共轭，g■■（x，y）为ge（x，y）的共轭。

1.3 快速傅里叶变换 进行离散傅里叶变换需要的计算量太大，运算时间太长，在某种程度上限制了它的使用。按照式（1），计算一个长度为N的一维离散傅里叶变换，对u的每一个值需要做N次复数乘法和（N-1）次复数加法。那么对N个u，则需要N2次复数乘法和N（N-1）≈N2次复数加法。很显然，当N很大时，计算量是相当可观的。

库里（Cooley）和图基（Tukey）在1965年首先提出一种快速傅里叶变换（FFT）算法，采用该算法进行离散傅里叶变换，复数乘法和加法次数正比于Nlog2N，这在N很大时计算量会大大减少，如表1。

可见，采用FFT可以减少运算量，图像越大减少越多。对于长为1024的离散序列，用普通的离散傅里叶变换往往要计算几十分钟，而采用FFT，一般只要几十秒。快速傅里叶变换是傅里叶变换的一种改进算法，它分析了离散傅里叶变换中重复的计算量，并尽最大的可能使之减少，从而达到快速计算的目的。

在这里，只对一维情况下的快速傅里叶变换进行讨论，因为根据傅里叶变换的分离性可知，二维傅里叶变换可由连续两次一维傅里叶变换得到。

设WN=e■（27）

设N为2的整数幂，即：N=2n，并令M为正整数，满足N=2M，则式（1）可写成：

F（u）=■■f（x）W■■=■■■f（2x）W■■+

■■f（2x+1）W■■（28）

由式（27）可知W■■=W■■，所以式（28）可写成：

F（u）=■■■f（2x）W■■+■■f（2x+1）W■■W■■

（29）

现在定义：

G（u）=■■f（2x）W■■ u=0，1，2…，M-1 （30）H（u）=■■f（2x+1）W■■ u=0，1，2…，M-1 （31） 式（29）可转化为：F（u）=■[G（u）+H（u）W■■]（32）

同样，因为W■■=W■■和W■■=-W■■，式（30）通过式（32）得到：F（u+M）=■[G（u）-H（u）W■■]（33）

仔细分析式（30）至式（33）可知，一个N点变换可以通过把原始表达式分成两部分来计算，如式（32）和式（33）所示。计算F（u）的前半部分要对式（30）和式（31）给出的两个N/2点变换进行计算。G（u）和H（u）的计算结果被代入式（32）中得到F（u），u=1，2，…，（N/2-1），另外一半可直接从式（33）得到，而无需另外的变换计算。因此可以将求N个点的离散傅里叶变换转换成求N/2点的离散傅里叶变换，当N为2的整数幂时，式中的G（u）和H（u）可以再被细分成更短的序列，进一步加快运算速度。

1.4 傅里叶变换在图像检索领域的应用 通过上述研究可知，图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度。如：大面积的沙漠、海水、天空在图像中是一片灰度变化缓慢的区域，对应的频率值很低；而对于地表属性变换剧烈的边缘区域在图像中则是一片灰度变化剧烈的区域，对应的频率值较高。从物理效果看，傅立叶变换是将图像从空间域转换到频率域，其逆变换是将图像从频率域转换到空间域。换句话说，傅立叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数，傅立叶逆变换是将图像的频率分布函数变换为灰度分布函数。

作者认为傅里叶变换具有以下优点：第一，傅里叶变换是在图像的频域对图像进行特征分析，对噪声不敏感，而基于空域的图像分析方法是很难控制噪声对图像带来的影响的；第二，一幅图像由空域转换到频域，使高度相关的空间样值变为相关性较弱的变换系数，减少了空间样值之间的冗余度。正是这些优点，使得傅里叶变换被广范地应用于图像处理和图像分析领域。一维傅里叶变换应用于物体形状的轮廓描述上，它首先要确定图像中物体的轮廓，然后将轮廓映射到复平面上后，再对其进行傅里叶变换。而二维傅里叶变换不需要确定物体的边界，可直接对图像进行变换，对图像分别从x和y方向进行傅里叶变换，应用于对物体形状的区域描述上。

2 基于快速傅里叶变换的特征提取

本文在快速傅里叶变换的理论基础上提出了一种基于区域形状特征的提取方法。该方法首先对图像进行二维的快速傅里叶变换，然后对变换的结果进行分块统计，由每块的频域统计值构成图像的形状特征向量，提取步骤如图1所示。

2.1 图像归一化处理 考虑到图像资源尺寸的不统一，使得傅里叶变换的结果矩阵大小也不统一，将造成后期无法对图像统一提取特征这一问题。作者对图像进行了大小的归一化处理。又考虑到离散快速傅里叶变换仅适用于大小为2n的图像，可以将图像归一化后的尺寸选择为64×64、128×128、256×256、512×512。

归一化的过程实际是对图像进行缩放，缩会使多点合并成为一个点，放会使图像产生空白点。作者采用的具体方法是：对多点合并一个点时，先计算周围8领域的平均值，然后再利用这个平均值作为合并点的值；对产生的空白点需要进行填充处理时，先计算周围8邻域像素点的平均值，然后利用该平均值填充。通过归一化处理的图像，从视觉上看是比较接近的，为后期特征的提取带来了方便。但作者也发现它有可能对图像检索产生一些负面影响，一是如果图像较大，图像压缩则可能丢失一定的图像信息；如果图像较小，这样的拉伸会使图像信息产生一些冗余的附加信息。二是由于原有图像的宽度和高度如果不是1：1的比例，这样的缩放会使图像的形状产生变形和扭曲。

在对图像进行归一化处理时，若原图像的比例不是1：1，则造成了图像的扭曲变形，为图像检索带来了负面影响。同时图像归一化的过小，图像可能会丢失掉一些重要的细节，影响图像的检索效果，而归一化的过大，又增加了傅里叶变换的计算量，影响图像的检索效率。作者权衡以上两个因素，将图像归一化为256×256个像素点。

2.2 对图像进行频域转换 傅里叶变换是将空间域转换为频率域的有效方法，利用傅里叶变换的这种特性就可以将图像的空间信息转换为频率信息进行处理。要提取图像的形状区域信息，则需要对图像进行二维离散傅里叶变换，1.1节中已经详细地介绍了二维图像傅里叶变换的公式，这里就不再重复。作者考虑到二维傅里叶变换需要的计算量依赖于图像的大小，计算处理时间随着图像的增大迅速增长。因此，作者采用了傅里叶变换的快速算法对图像进行二维傅里叶变换，这一做法使傅里叶变换的计算量降至Nlog2N，大大缩减了计算时间，其方法已在1.3节做了详细分析，这里不再重复。快速傅里叶变换虽然可以使计算量下降，但其局限性也是很显然的，它的N要求必须是N=2n，所以作者在2.1节对图像做的归一化处理时就考虑到了这个问题。

由1节可知傅立叶变换是基于复数的，其结果为复数形式。因此，要直接表示一幅图像的二维FFT变换的结果就必须用到两幅图像：一幅表示实部，一幅表示虚部。在频谱图像中，灰色表示的频域值为0，黑色表示的频域值为负值，白色表示的频域值为正值。由图2可以看到4个角上的黑色更黑，白色更白，表示其幅度更大，其实4个角上的系数表示的是图像的低频组成部分，而中心则是图像的高频组成部分，低频部分反映了图像的整体轮廓，高频部分反映了图像的细节。尽管这样表示的频谱图已经能够反映出原始图像的能量分布情况，但是不够清晰、直观，因此，另一种表示方法是引入变换结果的模作为值在频谱图中表示出来，用灰度的明暗来代表模的大小，如图3所示。

研究表明，二维快速傅里叶变换是一种正交变换，变换后信息不会丢失，具有能量保持性，并且能够对能量进行重新分配与集中，使高度相关的空间样值转变为了相关性较弱的变换系数，从而减少空间样值之间的冗余度，降低了对噪声的敏感度。同时，对不同的图像进行二维快速傅里叶变换后得到的频谱图是不同的，且二维快速傅里叶变换的结果直接反映了图像的形状区域信息，因此可以通过对变换后的图像频谱结果进行分析来提取图像的区域形状特征。

2.3 频域特征的提取 对一个N×N的灰度图像进行二维快速傅里叶变换后，结果为一个N×N的复数矩阵，它反映了图像中物体的区域形状信息。通过分析发现，直接从傅里叶系数计算频谱（即模）作为描述符，是不实用的，因为得到的特征不具备旋转不变性。对一个形状描述符来说，是否具有旋转不变是非常重要的，因为近似的形状可能具有不同的方向。作者通过对图2和图3（a）所示结果进行分析后发现，虽然由傅里叶系数实部、傅里叶虚部和系数的模绘制出的频谱图外观上虽然不尽相同，但它们的能量分布却是完全一致的，体现了相同的图像特征；通过实验，作者分析发现对变换的结果仅提取实部作为描述符就可以有效地解决旋转不变性问题，因此作者选择傅里叶系数的实部来描述形状的区域特征，具体做法如下：

由于对图像进行傅里叶变换前已经将图像归一化为256×256，因此这里得到的变换结果为一个256×256的矩阵，如果利用每一个频谱特征构成图像的特征向量，则一幅图像的数据量就达到216，大大增加了图像检索征匹配过程的计算量，降低检索效率。在这里，作者把矩阵均匀分割为16×16的块，如图4所示，共提取出256个特征值组成图像的区域形状特征向量。

由于每个块中包含了256个傅里叶系数，我们需要对它们进行计算得到该块的特征值，作者选择对每块内的256个系数统计求和再开方，这样做可以获得性能更一致的全方位的能量响应。方法如下：

假设傅里叶系数为：f=A+Bi

Fi=■A■■ i=0，1，2，…，255（34）

其中，Fi为第i个块，Aj为块中第j个傅里叶系数的实部，j的计算公式如下：

PixlNum[r][c][m][n]=（16*r+m）*256+16*c+n（35）

其中，r为第i个块所在的行，c为第i个块所在列，m为第i个块中第j个系数所在行，n为第i个块中第j个系数所在列。

由于傅里叶系数实部的值比较大，为了特征匹配后的结果易于理解和分析，作者对特征Fi进行了归一化处理，即S■=■（36）

这样就得到了图像频谱的256个特征分量S■（i=1，2，3，…，N），由它们构成图像的频谱特征向量，即图像的区域形状特征向量，在图像中的特征分量从左到右分别计为1-N。

3 特征匹配

图像检索时，作者用查询图像的频谱特征向量直接和图像库中的每一幅图像的频谱特征向量运用欧式距离法进行相似性比较，得到相似图像，并按相似性进行排序。

S （q，d）=■（S■-S■）■■（37）

式（37）中S （q，d）表示总得相似性距离，S■表示图像库中被检索的图像的第i个特征分量，S■表示待检索图像的第i个特征分量。

4 实验结果与分析

为了更好地评价本章算法的检索性能，本文选择了MPEG-7 ShapeB标准形状测试集。测试集中有1400幅图像，被分成70类，每类20幅图像。图像库内的图像均为灰度图像，且每幅图像目标唯一、背景单一，其中目标为白色，背景为黑色。

对性能的评价，本文采用“查准率”对本章算法输出的结果与人们期望结果的一致性进行比较，查准率的公式见式（38）。实验中，作者将P（G）（即输出近似图像的数目）分别选择为10和20。具体实验步骤如下：以MPEG-7 ShapeB形状库为基础，从中选择10类作为本实验的测试形状库，其中每类包含了20幅图像，这些图像的分类是基于语义的分类，具体分类如图5所示，其中每一类中所包含的图像在视觉上是相似的，一些不同类之间也有相似性，如camel和elephant。对每一类图像随机抽取10幅分别作为关键图像，用查准率公式计算出P10和P20，然后对每类计算平均查准率■10和■20。

precision=P（R|G）=■=■（38）

表2为本章算法（记作算法A）和采用傅里叶系数的模作为特征的检索算法（记作算法B）的平均检索准确率。

从表2中的结果可以看出，针对上述各测试图像类，本章算法的平均检索准确率均高于采用傅里叶系数的模作为特征值的算法。由于本章算法在频域提取图像的区域形状特征，故对噪声不敏感，具有平移、尺度不变性；只选用傅里叶系数的实部建立图像的特征向量对图像进行检索，获得了较好的旋转不变性，检索效果如图6，第一副为示例图像。

参考文献：

[1]章毓晋.图像处理和分析（图像工程上册）.北京：清华大学出版社，2004.

[2]张洁.数字图像边缘检测技术的研究.硕士学位论文，合肥工业大学，2009.

[3]Bober M.. MPEG-7 visual shape descriptor.IEEE Trans. on Cirrcuits and Systems for Video Technology.2001，11（6）：716-719.

[4]冈萨雷斯.数字图像处理.北京：电子工业出版社，2005.

[5]Latecki L.J. Shape Data for the MPEG-7 Core Experiment CE-Shape-1， http：//www.cis.temple.edu/～latecki/TestData/mpeg7shapeB.tar.gz，2002.

[6]Castleman K R著，等译.数字图像处理.北京：电子工业出版社，2002.