数字图像基础之图像变换

数字图像基础之图像变换

来源http://software.intel.com/de-de/node/256949

2.3 图像变换

 

2.3.1 图像变换的目的

有时候，直接对图像进行处理会遇到一些困难，如图像中周期干扰信号的滤除、图像的压缩编码等。而将图像数据转换为其他某些具有明显特征的形式时，可以找到简单而且有效的多的处理方法。因而在很多时候需要将图像的空间信息变换为其它形式的信息进行处理，将信息进行处理后，再通过对应的反变换将图像信息变换回空间信息的形式，我们将这种图像处理的过程通常称为图像的变换域处理（频域处理）。图像的频率处理在后续的多个章节都有应用，如第三章的图像编解码、第四章的图像频域增强处理等。

为了让大家对变换域处理的这种方法的目的和意义有一个感性认识，我们举一个一维信号变换域处理的例子：我们的生活中，有各种各样的信号，如声音信号。如图2.13所示为我们在某一海域的轮船上通过麦克风接收到的海里的声音信号（进行了适当的放大处理），从时间域（时间轴为自变量，横轴）上来看，我们不能从这些杂乱无章的信号里得出任何信息，但是，如果我们将这个信号通过适当的处理（乘上一个窗口信号，然后进行傅里叶变换（频域变换）、并通过减噪处理后），我们会发现这个语音信号在13Hz有一个尖峰，在26和39Hz也有较低的波峰，通过分析我们发现，这是一种周期信号（非正弦信号）的频谱，13hz称为基频，26和39Hz分别称为二次和三次谐频。基于此分析，我们可以推断这个声音的信号是海底潜水艇的三叶片推进器所发出的（推进器以每秒4.33转的速度旋转，每个叶片旋转速度为每秒13转）。这种频谱分析方法是声纳技术的基础。

图2.13 对麦克风接收的海里信号进行傅里叶变换处理过程

我们再以一个图像处理的实例来说明变换域处理的目的和意义。这个实例是用来验证是否在一个图像中存在一个预先定义好的模式，比如识别一个一块钱美元的硬币（如自动售货机），如图2.14所示为一个100*100像素的硬币图像。中间为乔治华盛顿的肖像，我们需要确认硬币中是否存在29*29像素大小的华盛顿的脸部肖像。因此，我们的问题定义是，给定一副图像（如一块钱硬币的图像）和一个模式（如图2.14中29*29像素大小的华盛顿脸部的肖像），怎样确认该副图像中是否存在该模式并且定位该模式的位置。

图2.14 识别货币中是否存在某一图像的处理过程

我们可以看出，这个问题是一个相关性的问题（correlation），（相关性是衡量两个对象的相似程度的一种方式，如果两个对象是相同的或相似的，那么他们的相关性最大，在数学中，我们可以用内积的（协方差）方式来判断两个向量的相关性，如果内积值为正值，说明两者正相关，如果为负值，则为负相关，如果为0，则表示不相关），而相关性的计算方式和卷积是相似的（请读者参考卷积运算定义）。因此可以用卷积来计算相关性，但是，卷积计算式一种非常耗时的运算，并没有多少优势，而在数学中，已经证明在时域或空间域中的卷积等于在频域中的乘积，而乘积的运算非常快，因此我们将图像进行傅里叶变换，再将两个变换的图像相乘，最后，用反傅里叶变换求的其相关性的值。我们从图中可以看出，最终在图像中有一个亮点，我们成功找出该模式并且将其准确定位。这只是变换的一种目的：方便计算，其实变换还有很多很重要的目的，比如图像压缩等，在图像的压缩存储中，我们认为图像中单个像素的信息时很少的，很大程度上图像中某一个像素和周围像素的相关性很大，因此，通常会将这种空间上相关性很大的信息通过变换转换成不相关的信息，从而通过降低信息的冗余度来实现压缩的目的。

信号的正交变换 设f(x,y)为任意一个空间域上的函数，有n个函数f₁(x,y),f₂(x,y),...,f_n(x,y)在区间[(x₁,y₁),(x₁,y₂)]构成一个正交函数空间，用这n个正交函数的线性组合近似f(x,y)得：

把这种将一个函数分解成正交函数的线性组合称为信号的正交分解，将一个函数通过正交分解映射到正交函数空间的数学变换称为正交变换。

2.3.2 图像变换的种类（DFT/FFT、DCT）

傅里叶变换 傅里叶（1768-1830）是法国的数学家和物理学家，他提出连续的周期信号都可以表示为不同频率的正弦和余弦的和形式，而每个正弦和余弦量前需要乘以不同的系数（称为傅里叶级数）。无论这个信号多么复杂，只要它是周期的（需要满足某些软的数学条件），都可以用这样的和来表示。如图2.15所示，左图中我们可以看出，其信号比较复杂，难以看出规律，通过傅里叶分析，我们可以得出，左边的信号是四个余弦信号的叠加形成的，这四个余弦函数的频率分别为10，25，50，100HZ（频率f=2π/T）。并且其系数均为1。因此，这个看似复杂的信号的可以表示成：

图2.15 一维信号的傅里叶变换

上述对周期信号的分解方式称为傅里叶级数。但是，现实生活中，我们遇见的信号通常都是非周期的，有一定定义域限制的，但仍可以用加权的正弦和余弦的和来表示，的我们在这种情况下的分析称为傅里叶变换。变换前的信信息是在时间域（声音信号）或空间域（图像信号），变换后的信息则在频率域（一维信号的频率域信息为频幅，频率大小和相位，二维信号的频率域信息为频幅，频率大小（含方向）和相位）。经过变换后，我们可以再频率于对信号进行处理，再转换回函数的原始域时不会丢失任何信息。

傅里叶变换是一种线性积分变换，在计算机中只能将积分进行离散近似。因此计算机实现的傅里叶变换通常称为离散傅里叶变换（DFT，Discrete Fourier Transformation），由于离散傅里叶变换的速度比较慢，在上世纪60年代，发明了离散傅里叶变换的快速实现方式，我们称之为快速傅里叶变换（FFT，Fast Fourier Transformation）。傅里叶变换在物理学、信号处理、通讯等众多领域有着广泛的应用，是数字信号处理领域一种很重要的算法。

实数一维离散傅里叶变换（DFT，Real Discrete Fourier Transformation） 为了让读者较好理解图像（二维）的傅里变换，我们先介绍一维傅里叶变换。二维的傅里叶变换在一维的基础上进行的延伸。理解一维傅里叶变换是二维傅里叶变换的基础。

在一维空间的信号通常又称为时域信号（time domain signals），因为信号的采集是每个一定时间间隔进行的（regular intervals of time），一维离散傅里叶变换通常将时域信号分解（decompose，analysis, forward DFT）成不同频率的正弦和余弦的系数。一维信号通常用小写的数组表示，如x[]，信号的长度通常由变量N来表示，而信号的频域信息（frequency domain）通常用大写的数组表示，如X[]，频域信息包含两个部分：ReX[]（Real part of x[]）和ImX[]（Imaginary part of x[]）。ReX[]对应于不同的余弦函数的幅度（amplitudes），而ImX[]对应于不同的正弦函数的幅度。假设有一个时长为1秒、采样频率为128Hz的信号（N=128，），由采样定理可知，该信号的带宽（bandwidth）应为0至截止频率，该截止频率为采样频率的一半（64Hz，该频率通常又称为Nyquist frequency），如果该信号中的高频成份的频率大于64Hz，则用128Hz的采样频率进行采样，就会产生著名的“混叠效应（aliasing）”。这也就是说，在对信号进行傅里叶变换时，最大的频率应为N/2。所以在对该信号进行傅里叶变换后，ReX[]从ReX[0]至ReX[64]，ImX[]从ImX[0]至ImX[64]。而在傅里叶变换中用到的正弦和余弦函数通常也称为基函数（basis functions） ，这些正余弦函数的幅度是都是基本单位1，设信号的角频率为w₀=2π/N，则进行傅里叶变换的基频为w₀，而最高频率为w=(2π/N)*(N/2)。因此，这些基函数可以描述为：

DFT的计算可以通过“相关性计算（Correlation）”来实现，公式如下：

上述两个公式的数学意义是：计算每个公式中两个函数（向量）的相关性，相关性越大，说明其越相似，这两个函数一个是信号本身x[i]，而另一个是某一个频率的正弦和余弦函数，通过相关性计算，就能得出信号中是否含有该频率成份，成份的多少表现在相关性的计算结果ReX[k]和 ImX[k]，即每个频率对应的幅值。

采用上述的计算方式对式所形成的信号进行分解，采用采样频率为300Hz后，所形成的Re[k] (k=0,1,...,150)中，有ReX[10]、ReX[25]、Re X[50] 、Re X[100]为1，其它均为0 ，同时ImX[k] (k=0,1,...,150)均为0。

复数一维离散傅里叶变换（DFT，Complex Discrete Fourier Transformation），在实际的应用中，并不是应用上述的实数傅里叶变换，而是采用其扩展形式：复数的傅里叶变换。

复数的傅里叶变换公式如下:

其反变换的公式表示如下：

在上式中，为简化表示，我们通常将e^-j2π/N用W_N代替，W_N^-1=e^j2π/N，则上述的公式也可以用矩阵表示如下：

其中，满足条件。

复数形式的傅里叶变换和实数形式的傅里叶变换的一个最大的不同在于其频率可以去负值。由于实数形式的傅里叶变换有着明显的物理意义，频率均为正值（变现为K的取值为从0到N/2）。而根据复数形式的傅里叶我们可以看出，其K的取值范围为从0至N，其前0- N/2的频率为正值，且N/2为最大频率。而由公式可得傅里叶变换的频谱亦是周期的（周期为N），因此，N/2-N-1的频率相当于- N/2至0的频率。这也就是说，在对信号进行复数的傅里叶分析后，其高频信息部分对应于X[k]的中间位置。对应的频率由中间向两边递减。

二维离散傅里叶变换 离散傅里叶变换很容易扩展到二维，一个尺寸为M*N大小图像的二维傅里叶变换的公式为（我们将二维函数表示为f(x,y)，其傅里叶变换为F(u,v)）：

同样，给出F(u,v)，可以通过傅里叶反变换获得f(x,y)，由下列表达式给出：

对图像进行傅里叶分析后所得的数据如下：

图像的实部数据（Real Part）为: R(u,v)=Real[F(u,v)]；

图像的虚部数据（Imaginary Part）为：I(u,v)=Imag[F(u,v)]；

图像频谱的幅值为：；

图像频谱的相位为：；

图像频谱的能量：；

图像频谱的总能量：。

上述公式具有“可分离性”，二维DFT变换公式还可以写成：

而中，y为常量。 因此，根据分离性，可以将二维离散傅里叶变换写成矩阵的形式如下：

计算二维图像傅里叶变换的过程是比较简单的：由于傅里叶变换的分离特性，我们很容易从一维傅里叶变换得到二维傅里叶变换。我们可以首先计算每一行的一维傅里叶变换，然后进行每一列的一维傅里叶变换，变换的结果即为二维的傅里叶变换（也可以先计算列，后计算行，傅里叶变换的可分离性）。其傅里叶反变换的过程是一样的（首先对行进行反傅里叶变换，然后对列进行反傅里叶变换，或者次序颠倒）。

和一维傅里叶变换类似，二维离散傅里叶变换的结果中，F(u,v)在F(M/2,N/2)处对应最高频率成份的分量，而离该位置越远，所对应的频率越低，在同多傅里叶变换对图像进行频域的滤波时，理解这些物理意义是至关重要的。我们会在后面的实验中有所体现。

离散余弦变换 离散余弦变换（Discrete Cosine Transform，简称DCT变换）是一种与傅立叶变换紧密相关的数学运算。在傅立叶级数展开式中，如果被展开的函数是实偶函数，那么其傅立叶级数中只包含余弦项，再将其离散化可导出余弦变换，因此称之为离散余弦变换。

离散余弦变换（DCT）是N.Ahmed等人在1974年提出的正交变换方法。它常被认为是对语音和图像信号进行变换的最佳方法。离散余弦变换（DCT）在目前图像编码中享有重要地位，成为H.261、JPEG、MPEG 等国际上公用的编码标准的重要环节。在视频压缩中，最常用的变换方法是DCT。

一维离散傅里叶变换（The One-Dimensional DCT） 一维离散余弦变换的公式表示如下：

将离散余弦变换和离散傅里叶变换进行比较可以看出，其基函数是不同的：首先，其没有正弦的基函数，因为它假设信号是偶函数，对偶函数进行傅里叶变换后正弦函数的分量都会为0；其次，其余弦基函数的频率不同，虽然频率都是从0至截止频率，但频率间隔也不同（间隔要小一半），DFT从余弦基函数实现的是从N各数到N/2+1个数的转换，而DCT从余弦基函数实现的是刚好从N个实数到另N个实数的转换。

二维离散傅里叶变换（The Two-Dimensional DCT）二维离散余弦变换的公式表示如下：

对图像进行二维DCT变换后，图像的频谱数据（频幅）出现在图像的左上角上，而高频数据出现在图像的右下角上，离左上角距离越远，所对应的频谱越大。在图2.16图像的分析中，我们根据对图像进行DCT变换得出一种有用的统计规律，即图像的能量谱集中在低频部分，这也是图像数据压缩的重要理论依据。

对于离散的数字信号，余弦变换为离散形式的DCT（discrete cosine transformation）。它是图像压缩中非常常用的一种变换编码方法。

图2.16 一副图像及其离散余弦变换后频谱图

2.3.3 图像变换的应用

正交变换在信号处理、通讯等领域有着广泛的应用。在图像处理领域，图像变换时其不可缺少的部分，也是数字图像处理的重要工具。在图像处理领域，应用最多的图像变换有傅里叶变换，离散余弦变换等。

二维图像的傅里叶变换和一维傅里叶变换类似。但是不同的是，将一维傅里叶变换有着更加直接的意义。举例来说，将声音信号进行傅里叶变换会使在时域中杂乱无章的信号在频域中变得更加容易理解（如上述2.3.2节所阐述的声纳技术）。但是对已经很直观的二维图像信息进行傅里叶变换，很难观察频域中信息的意义。在图像中，边缘信息是图像的主要信息，边缘信息会有很多不同的频率成份，因此也很难通过改变频谱的方式来处理图像。因此，图像滤波器通常在直观的空间域进行设计的，图像的平滑和增强都针对的是空间域的效果，而不是高通和低通滤波器。

尽管这样，二维图像的傅里叶变换仍然有一些很有用的属性。比如，在空间域中的卷积运算可以转化为对应频域中的乘积，由于乘积的操作比卷积快很多。这种特性可以用于加入图像处理的计算速度。

离散余弦变换是傅里叶变换的近亲。傅里叶变换用一系列的正弦和余弦来的和来表示信号，而余弦变换只用余弦信号来代表信号。余弦变换有好多个版本，之间有一些小的差异。离散余弦变换是数码率压缩需要常用的一个变换编码方法。DCT变换在数字图像的编码中得到很好的应用，大家都熟悉的JPEG格式的图像文件，其压缩算法中就使用了DCT变换。JEPG压缩算法先将整体图像分成N*N像素块，然后对N*N像素块逐一进行DCT变换。由于大多数图像的高频分量较小，相应于图像高频分量的系数经常为零，加上人眼对高频成分的失真不太敏感，所以可用更粗的量化。因此，传送变换系数的数码率要大大小于传送图像像素所用的数码率。到达接收端后可以通过反离散余弦变换回到样值，虽然会有一定的失真，但人眼是可以接受的。

2.3.4 IPP库中的二维图像线性变换函数

在本节，我们介绍怎样调用IPP库函数在图像内存数据上实现线性变换。这些变换包括快速傅里叶变换（Fast Fourier Transform ，FFT）、离散傅里叶变换（Discrete Fourier Transform，DFT）和离散余弦变换（Discrete Cosine Transform，DCT）。

傅里叶变换 我们首先来介绍一下进行傅里叶变换所需的数据结构。我们首先分析一下傅里叶变换的过程。

根据公式：

及其矩阵形式：

，我们以一个4*4大小的图像为例，图像矩阵X、矩阵F₄分别为

因此其傅里叶变换为：

其实部和虚部分别为：

我们可以看出具有某些对称规律，我们再用一个示例来详细说明。对如下列矩阵的8*8图像数据经过上述的傅里叶变换。

我们求得其实部Xr和虚部Xi的矩阵值分别为：

图2.17 图像傅里叶变换实部对称图

我们对求得的实部和虚部进行分析，得出其如下的对称规律（如图2.17所示，以实部为例）：在用标号对应的每个个部分，都是对称的。虚部有一点不同，其相反数是对称的。

在将一个二维实数的图像数据进行傅里叶变换后， Intel® IPP采用包数据格式“RCPack2D”来存储和索引这种复数矩阵，在傅里叶反变换中，会将这种包数据转化为实数数据。

实数-复数包格式（Real - Complex Packed (RCPack2D) Format） RCPack2D格式利用了共轭对称的特性，因此只需存储矩阵的一半数据（另一半数据为其堆成的数据的共轭），对应于一个 的变换，傅里叶变换将所求得的复数的傅里叶系数存放在一个 的实数数组中。针对 是偶尔还是奇数，所存储的值有一点不同，如下表所示（M为偶数的时候）：

表2.5 M为偶数时RCPack2D格式存储的图像实部和虚部数据布局

因此针对上述Xr和Xi矩阵，针对上表存放的规律，ipp中RCPack2D的存储值如下（用实验验证了其值正为调用IPP中的傅里叶变换后所得的结果值）：

针对M为奇数的时候，RCPack2D格式存储的元素布局如表2.6所示：

表2.6 M为奇数时RCPack2D格式存储的图像实部和虚部数据布局

上述的格式针对的是单一通道，如果是针对RGB三个通道，其存储的结构将是每个通道的数据排放在一起，采取如下的格式：通道1-Re A(0,0)；通道2-ReA(0,0)； 通道3-Re A(0,0); 通道1-Re A(0,1)；通道12-Re A(0,1)；通道3-Re A(0,1)；通道1-Im A(0,1)；通道2-Im A(0,1)；...（详细介绍请参考IPP使用手册。）实现傅里叶变换所需的一些函数如表2.7所示。

调用傅里叶变换函数的一段代码如下图所示：我们来分析一下调用傅里叶变换所需要的一些步骤，如表2.8所示。

表2.7 IPP库图像傅里叶变换相关函数

函数名	函数解释
FFTInitAlloc	分配和初始化图像傅里叶变换所需要的内存空间
FFTFree	释放在傅里叶变换中所需的内存的空间
FFTGetBufSize	计算傅里叶变换所需额外空间的大小
FFTFwd	对一幅图像实现傅里叶变换
FFTInv	对复数数据进行傅里叶反变换转化为图像数据
MulPack	将两个以packed格式的图像相乘结果放入另一个packed格式的图像中

表2.8 IPP库图像傅里叶变换函数调用步骤

步骤	调用的函数
1、初始化图像数据，定义目标数据
2、申请IppiFFTSpec_R_32f变量，初始化内存分配，分配在计算傅里叶变换中所需的额外内存	ippiFFTInitAlloc
3、调用傅里叶变换函数	ippiFFTFwd_RToPack
4、释放申请的额外内存	ippiFFTFree_R_32f

基于以上步骤，我们实现一个简单的傅里叶变换的代码如下：

离散余弦变换 由于离散余弦变换将实数的二维图像数据转化为实数数据（虚部全为0），

不需要刚才所需的包格式结构来存储结果，因此要简单一些，离散余弦变换在图像编码中经常有应用，在编码领域，经常需要将8*8的图像进行离散余弦变换，IPP库能够直接支持这样的转换，其中有一个函数ippiDCT8x8Fwd实现了这种转换，上面是这种转换的简单代码。