2.3 数字图像处理的数学工具

2.3.1 傅里叶变换图像处理

图像频域变换工具主要有傅里叶（Fourier）变换、余弦变换和小波变换等，其中傅里叶变换是最基本的一种，在数字图像处理中应用广泛。傅里叶变换对图像的分解可比喻成一个玻璃棱镜对光信号的分解。棱镜是可以将光分解为不同颜色的物理仪器，每个成分的颜色由波长（或频率）来决定。傅里叶变换可以看作数学上的棱镜，能将函数基于频率分解为不同的成分，从而可以更好地识别图像中的各种成分。

图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度。如：大面积的沙漠在图像中是一片灰度变化缓慢的区域，对应的频率值很低，而地表属性变换很大的边缘区域在图像中是一片灰度变化剧烈的区域，对应的频率值较高。傅里叶变换有非常明显的物理意义，设f是一个能量有限的模拟信号，则其傅里叶变换就表示f的谱。从纯粹的数学意义上看，傅里叶变换是将一个函数转换为一系列周期函数来处理的。从物理效果看，傅里叶变换是将图像从空间域转换到频率域，其逆变换是将图像从频率域转换到空间域。换句话说，傅里叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数，傅里叶逆变换是将图像的频率分布函数变换为灰度分布函数。

（1）图像傅里叶变换

傅里叶变换是数字图像处理中应用最广的一种变换，图像增强、图像复原和图像分析与描述等每一类处理方法都要用到图像变换，尤其是图像的傅里叶变换。傅里叶变换将时域信号分解为不同频率的正弦和余弦和的形式。它是数字图像处理技术的基础，其通过在时域和频域来回切换图像，来对图像的信息特征进行提取和分析。

离散傅里叶变换的定义如下。

二维离散傅里叶变换为

逆变换为

式中，u，x∈{0，1，…，M-1}，v，y∈{0，1，…，N-1}。在离散傅里叶变换对中，F（u，v）称为离散信号f（x，y）的频谱，而|F（u，v）|称为幅度谱，φ（u，v）为相位角，功率谱为频谱的平方，它们之间的关系为

图像傅里叶变换与反变换效果如图2-16所示。

图像傅里叶变换的matlab程序如下。

（2）傅里叶频谱性质以及图像旋转

对灰度图像进行二维傅里叶变换，对变换频谱进行观察，会发现一些性质。图2-17为图像旋转后频谱图的比较。从图中可以发现，随着图像的旋转，图像的纹理方向发生变换，同时图像的频谱也跟着旋转。实际上，图像纹理旋转了多少度，其对应的频谱图像也旋转相应的角度。

傅里叶频谱性质MATLAB程序如下：

（3）傅里叶变换滤波

傅里叶变换对信号进行高通、带通、低通滤波可实现混合信号的分离。如图2-18所示，原始图为低频信号与高频信号（条纹纹理）的组合，可使用二维傅里叶变换实现信号的分离。

图像傅里叶变换滤波实现信号分离的程序如下。

2.3.2 离散余弦变换

离散余弦变换（Discrete Cosine Transform，DCT）是与傅里叶变换相关的一种变换，它类似于离散傅里叶变换，但是只使用实数。离散余弦变换相当于一个长度大概是它两倍的离散傅里叶变换。离散余弦变换的类型包括Ⅰ型、Ⅱ型和Ⅲ型等，其中最常用的是第二种类型。它的逆被称为反离散余弦变换或逆离散余弦变换（IDCT）。

Ⅱ型和Ⅲ型离散余弦变换公式分别为

式中，F（u，v）为f（x，y）的DCT变换；（x，y）为图像空间坐标；（u，v）为变换域坐标；N为图像的大小。

离散余弦变换经常在信号处理和图像处理中使用，用于对信号和图像进行有损数据压缩。这是由于离散余弦变换具有很强的能量集中特性：大多数自然信号（包括声音和图像）的能量都集中在离散余弦变换后的低频部分。离散余弦变换常被用于JPEG图像压缩。此外，它也经常被用于求解偏微分方程。

图2-19所示为图像经过DCT变换，变换系数处理以及反变换重构后的结果。图2-20所示为采用8×8离散余弦变换压缩算法对图像进行压缩的结果。从图2-20可以看出采用不同的压缩参数（程序中的mask）时压缩效果具有差别。图2-21为基于离散余弦变换的相位展开算法处理结果，其中采用离散余弦变换求解偏微分方程进而实现了相位展开。

图2-19中的处理实现程序如下。

图2-20中的处理程序如下。

图2-21中的处理程序如下。

2.3.3 偏微分方程图像处理

偏微分方程图像处理是这几年兴起的一种图像处理方法，主要针对底层图像处理，在图像去噪、修复、分割等方向的应用中取得了不错的效果。偏微分方程具有各项异性的特点，在图像处理中可以在去噪的同时很好地保持边缘。比如，在图像修复中，利用偏微分方程对图像进行建模，可以使得待修复区域周围的有效信息沿等照度线自动向内扩散，在保持图像边缘的基础上平滑噪声。

偏微分方程在图像处理中的成功应用可以归功于两个主要因素。首先，许多变分问题或它们的正则化逼近通常可以通过其欧拉-拉格朗日方程进行有效计算。其次，像经典的数学物理模型一样，偏微分方程对于描述、建模，以及模拟许多动态和平衡的现象（如扩散、对流或输运、反应等）非常有效。

图像处理中的偏微分方程与变分模型存在一定关系。一方面，图像处理中的偏微分方程不一定总满足一个变分模型。这在物理中是常识，著名的例子有流体力学中的纳维-斯托克斯方程、量子力学中的薛定谔方程、电磁学中的麦克斯韦方程。另一方面，在图像处理中，许多偏微分方程都有对应的变分模型，比如热扩散方程、平均曲率运动方程和全变分方程。下面介绍几种典型的偏微分方程图像处理模型。

1. PM模型及其改进模型

Perona和Malik对最早的各项同性热传导方程滤波模型进行了重要改进，将控制扩散速度函数g（|▽u|）引入模型中，又于1990年提出了第一个具有方向性的偏微分方程滤波模型——各向异性扩散模型，即PM模型。该模型用一个非线性方程代替热传导方程

式中，函数g（|▽u|）用于控制扩散速度，是图像梯度▽u的非增函数，因此扩散速度在图像边缘处（相当于▽u大的地方）慢，使图像边缘能够得到保持，而在均匀区域扩散速度较大（即▽u小的区域），能够更好地平滑非边缘区域。该模型极大地改进了滤波效果。

常用的控制扩散速度函数为

式中，k为预先设定的常数。在PM模型中，扩散根据控制扩散速度函数进行，降低了边缘附近的平滑作用，但它对噪声极为敏感，许多学者仍在不断地对该函数进行改进。

复数扩散偏微分方程是PM模型的改进之一，其目的是避免或者减少PM模型的块相应。其模型为

式中，u（x，y，t）为随时间t演化的图像；Im（u）为图像u的虚部；▽·为散度操作符；D扩散速率定义为

式中，i为虚数单位；k为阈值参数；θ∈（-π/2，π/2），为相位角。其中不包含图像的导数项，这是比实数扩散方法优越的一点。Gilboal等研究结果表明，复数扩散偏微分滤波在光滑区域幅度加强，在边缘区域幅度减弱。

2. 全变分（Total Variation，TV）模型

基于变分方法的全变分模型最早起源于Osher，Rudin，Fatemi的工作，于1992年提出，简称为“ROF模型”。ROF模型中图像的光滑部分使用TV范数描述，纹理部分采用L2（平方可积）刻画，因此也叫作TV-L2模型。TV-L2模型的能量泛函表示为

式中，f为原始图像；u为带求解图像；λ为正则参数。

L2范数对于零均值的Gaussian信号具有很好的描述。

3. 热扩散方程的图像去噪

热扩散方程表示为

式中，Δ代表拉普拉斯算子；代表热扩散方程演化图像。

首先对时间项离散化

式中，n为离散化的迭代次数；Δt为离散化的时间步长，是一常数，比如取为0.1。

然后，对x和y离散化（采用前向差分，对应为i和j）

进一步可表示为

所以可以用简单的数组与循环来实现热扩散方程。图2-22为热扩散方程图像去噪结果。

图2-22所用的程序如下。

4. 全变分图像去噪

全变分图像去噪的能量泛函表示为

式中，u_x和u_y为u关于x和y的偏导数。

对应的欧拉-拉格朗日方程求解过程为

从而可以用简单的数组与循环实现全变分模型的演化方程。图2-23为全变分图像去噪结果。

全变分图像去噪的程序如下。

2.3.4 小波变换等时频分析方法

（1）小波变换

傅里叶变换将信号展开为无限个正弦和余弦的线性组合，这样虽可得到信号的频谱信息，但是不能获得信号时间方面的信息。换句话说，傅里叶变换提供了图像中所出现的所有频率，但是并不能反映它们源于何处。窗口傅里叶变换（即短时傅里叶变换）可以定位信号的变化。它将信号分解为小窗口并将其看作周期函数做局部处理。然而，窗口傅里叶变换仍然存在很大缺陷：窄窗口带来较差的频率分辨率，而宽窗口带来较差的定位效果。小波变换比窗口傅里叶变换更进一步，不仅能提供较精确的时域定位也能提供较精确的频域定位。

小波变换兴起于20世纪80年代中期，由于其具有良好的时频局域分析能力，所以成为信号处理最重要的工具之一。对含“点奇异”的一维信号，小波具有最优的逼近性能。小波变换具有如下优点。

● 稀疏性：小波系数稀疏分布，经小波变换后，只用很少的非零系数即可实现对对象的表达。

● 多分辨率：对不同分辨率根据信号和噪声分布消除噪声。

● 去相关性：小波变换可对信号去相关，因此相对于时域，更有利于去噪。

● 基函数多样性：可根据信号特点，结合滤波效果，来选择或构造不同的小波基函数。

如果函数ψ（ω）满足容许条件

则积分变换

为f（x）以ψ（x）为基的小波变换，引入符号ψ_a，b（x）

对f（x，y）∈L²（R²），其对应的小波变换为

式中，a为尺度因子；b、b₁、b₂为平移量。

以维纳滤波器为例

式中，y（t）=x（t）∗h（t）。

小波滤波的基本实现过程如下。

第一步：对图像f（x，y）进行小波变换，得到各级变换系数

第二步：对高频系数进行修正处理，可以得到系数的修正值

式中，η()为修正函数。

第三步：用a_-J（k）及进行重构，得到滤波后的图像。

由一维小波张成的二维小波只具有有限方向，即水平、竖直、对角，多方向的缺乏导致其不能最优地表达线或面奇异性，对线或面奇异函数的逼近性能差强人意。

设f（x）∈L²（R²），定义连续小波变换：

式中，，ψ()为母小波；a为尺度因子；b为平移因子。参数a和b均连续变化，所以称为连续小波变换。其内积表示形式为

参数a的变化决定着小波窗函数的形状和频谱结构。当a减小时，ψ_a，b（x）的频谱集中于高频部分，此时窗口尺寸较小，小波函数具有较好的空间分辨率；当a增大时，ψ_a，b（x）的频谱集中于低频部分，此时窗口尺寸变小，空间分辨率降低。

小波函数必须满足容许条件，小波变换才存在逆变换，其容许条件为

式中，。容许条件表明，可以用作基本小波的函数ψ_a，b（x）必须满足Ψ（0）=0的条件，即ψ_a，b（x）为均值为0的振荡波形，Ψ（ω）具有带通性。

（2）曲波变换

小波变换在图像边缘信息表达上存在严重不足，由一维小波张成的可分离小波只具有有限的方向性，不能“最优”地表示具有线性奇异性或面奇异性的高维函数。小波方向性的单一使得其不能很好地逼近具有高维奇异的信号。由于双树复小波的多尺度、近似平移不变性和多方向性，使得其比传统小波在图像去噪时更易保持边缘。双密度双树离散小波变换将信息提高到16个方向，但每个方向值由一个小波表示，使图像的分解与重构精度受到限制。

多尺度几何分析具有三个优点：即多分辨率分析特性、局域性和方向性。目前，多尺度几何分析工具主要有Brushlet变换、Ridgelet变换、Ridgelet变换、Curvelet变换、Wedgelet变换、Beamlet变换、Bandlet变换、Contourlet变换、Directionlet变换以及Shearlet变换等。由于它们主要以变换为核心，因此也称多尺度方向变换，为了能充分利用原函数的几何正则性，这些变换的基的支撑区间表现为“长条形”，用最少的系数来逼近奇异曲线。

定理2.1 设，令，其中曲线γ二阶可导，则函数f的曲波变换的m项非线性逼近误差为

边缘属于突变信号，是弯曲的，而不是直线，这些信号在数学上表现为奇异性。小波、脊波等不能有效表达出图像的这些特征。和小波相比，小波系数属于“过”边缘的表达，而曲波是“沿”边缘的表达，如图2-24所示，所以传统小波在处理图像上表现出很大的局限性。

在曲波变换尺度体系中，基函数本身对方向高度敏感。从图2-24中可以看出，在同一尺度下，两个不同位置的基元在方向上的表现只是略有差别，可见曲波具有很强的方向性。经过曲波分解，如果曲波基元与图像的边缘重合，那么在曲波系数上表现为大值，不重合则曲波系数很小，所以曲波系数能表达出该尺度、该方向上的边缘信息。更进一步，为了捕获图像各个方向的边缘，曲波可以进行不同尺度的分解，经旋转、平移后实现边缘特征的捕获，在矩阵上表现为稀疏性。

Candes等提出了两种曲波变换的快速实现方法，分别是USFFT（基于非均匀采样的快速傅里叶变换，Unequally Spaced FFT）和Wrap（基于频域特殊采样的卷绕规则Wrapping算法，Wrapping-based FFT）算法。这两种算法的主要区别在于不同尺度、不同方向上空间网格的选择方法不同，算法实现步骤简单介绍如下。

基于USFFT的离散Curvelet变换实现基本过程如下。

第一步：对输入的笛卡儿坐标下的二维函数（图像）进行快速傅里叶变换，得到二维频域表示

第二步：在频域，对每一个（j，k）重采样得到采样值

第三步：将内插后的与窗函数相乘便可得到

第四步：对进行二维快速傅里叶逆变换，可以得到离散Curvelet系数c^D（j，l，k）。

基于Wrapping的离散Curvelet变换实现方法核心思想是：围绕原点包裹（wrap），即在具体实现时对任意区域通过周期化技术一一映射到原点的仿射区域。该方法是在USFFT方法的基础上增加wrap步骤，具体过程如下。

第一步：对输入的一个笛卡儿坐标下的二维函数图像进行二维快速傅里叶变换，得到二维频域表示

第二步：在频域对每一个（j，k）重采样得到采样值

第三步：将内插后的与窗函数相乘便可得到

第四步：围绕原点wrapping局部化；

第五步：对进行二维快速傅里叶逆变换，可以得到离散曲波系数c^D（j，l，k）。

（3）剪切波变换

剪切波变换是一类新的多维函数逼近方法。剪切波通过一个基本函数的膨胀、剪切和平移变换来构造，体现了函数的几何和数学特性，如近年来许多领域的学者所强调的函数方向性、尺度和振荡等。剪切波变换是一种继承曲波和轮廓波优点的新型多尺度几何分析工具。

在L²（R²）空间中的二维函数f剪切波变换定义为

式中，Ψ_a，s，t为剪切波，它构成了在频域中连续尺度a＞0，位置t∈R²，曲率方向s∈R的局部化仿射系统。剪切波表示为

其中

它们分别为尺度矩阵和剪切矩阵。尺度矩阵A_a体现各向异性膨胀矩阵，剪切矩阵通过s参数化方向。

通过对尺度参数、剪切参数和平移参数离散化，得到离散的剪切波函数为

其中

相应的离散剪切波函数

对于ξ=（ξ₁，ξ₂），ξ₁≠0，令

式中，，，⊂[-1/2，1/16]∪[1/16，1/2]，并且⊂[-1，1]，能够得到以及。

假定

并且对∀j≥0

根据和的支集就能够得到函数的频域支集为

也就是说，是一个大小为2^2j×2^j的梯形对，方向为2^-jl。

和曲波相比，剪切波的构造有着根本的区别，剪切波是由一族算子操作在一个函数上产生的，而曲波并不是这样。图2-25为曲波和剪切波的频域剖分图。曲波和剪切波的不同之处在于：①曲波不是和一个固定的平移格相关联，如果考虑方向的敏感性，则剪切波的方向数目在每一个尺度加倍，而曲波则在每一个不同的尺度加倍；②剪切波定义在笛卡儿域，通过剪切变换获得各种方向，而曲波是在极坐标域构造的，通过旋转操作获得方向。

2.3.5 形态学处理

形态学本来指的是生物学中的一个分支，而此处的形态学为图像处理中的常用算法工具，它分析的基础是基于形态的数学运算。作为数学理论的一部分，它常用来调整分割后的图像形状以达到特定的目的，边缘检测、噪声抑制、特征提取甚至图像分割都可以将形态学运用于其中并且实际效果良好。

数学形态学是由一组形态学的集合运算组成的，其基本运算有膨胀、腐蚀、开运算/闭运算、击中/击不中、细化/粗化。这些基本运算在二值图像和灰度图像中各有特点。基于这些基本运算还可推导和组合出各种数学形态学实用算法。膨胀和腐蚀是数学形态学方法中最基本的运算。膨胀和腐蚀的原理是利用一个称作结构元素的“探针”收集图像的信息，当探针在图像中不断移动时，通过简单的逻辑运算便可考察图像各个部分之间的相互关系，从而了解图像的结构特征。

腐蚀操作的运算符号为“Θ”，如果对集合A用结构元素B来进行腐蚀操作，则可以表示为AΘB，并且定义

式（2-75）表明，结构元素B在集合A中平移x扫描后，结构元素B依旧在集合A中的所有参考点集合即为运算结果，具体操作例如图2-26所示过程。

图2-26a中的阴影区域为集合A，图2-26b中的阴影区域为结构元素B，设定红色“+”号标记点为参考点，图2-26c中的纯黑部分为腐蚀运算后集合A所留下的部分，即腐蚀的结果。

膨胀操作的运算符号为“⊕”，如果对集合A用结构元素B来进行膨胀操作，则可以表示为A⊕B，并且定义

膨胀运算与腐蚀运算略有不同，需要先对结构元素B进行关于原点的映射得到，结构元素在集合A中平移x扫描后，的位移与A至少有一个非零元素相交时B的参考点位置集合即膨胀结果，具体操作例如图2-27所示过程。

图2-27a中的阴影区域为集合A，图2-27b中的阴影区域为结构元素B，设定红色“+”号标记点为参考点，图2-27c中的阴影区域为结构元素B的映射，图2-27d中的纯黑区域为膨胀运算后集合A所增加的部分，整个阴影区域为膨胀的结果。

一般情况下，腐蚀运算和膨胀运算虽然不互为逆运算，但是可以通过恰当的组合变成新的形态运算，比如开运算和闭运算。使用相同的结构元素对图像先进行腐蚀运算再进行膨胀运算即为开运算，其运算符号为“°”，如果对集合A用结构元素B来进行开运算操作，则可以表示为A°B，并且定义

开运算的作用较多，比如可以削弱图像中狭长的区域，断开不同区域之间的连接等。虽然作用与腐蚀运算差不多，但是开运算可以基本保持待测物体尺寸不变。

使用相同的结构元素对图像先进行膨胀运算再进行腐蚀运算即为闭运算，其运算符号为“·”，如果对集合A用结构元素B来进行闭运算操作，则可以表示为A·B，并且定义

闭运算的作用也较多，也可以平滑图像轮廓，但与开运算不同的是，闭运算一般用来连接断开的临近区域。同样，它的作用虽然与膨胀运算差不多，但是可以基本保持待测物体尺寸不变。