基于小波变换的图像处理

基于小波变换的数字图像处理
摘要：本文先介绍了小波分析的基本理论，为图像处理模型的构建奠定了基础，在此基础上提出了小波分析在图像压缩，图像去噪，图像融合，图像增强等图像处理方面的应用,最后在MATLAB环境下进行仿真，验证了小波变化在图像处理方面的优势。

关键词：小波分析；图像压缩；图像去噪；图像融合；图像增强
引言
数字图像处理是利用计算机对科学研究和生产中出现的数字化可视化图像
信息进行处理，作为信息技术的一个重要领域受到了高度广泛的重视。

数字化图像处理的今天，人们为图像建立数学模型并对图像特征给出各种描述，设计算子，优化处理等。

迄今为止，研究数字图像处理应用中数学问题的理论越来越多，包括概率统计、调和分析、线性系统和偏微分方程等。

小波分析，作为一种新的数学分析工具，是泛函分析、傅立叶分析、样条分析、调和分析以及数值分析理论的完美结合，所以小波分析具有良好性质和实际应用背景，被广泛应用于计算机视觉、图像处理以及目标检测等领域，并在理论和方法上取得了重大进展，小波分析在图像处理及其相关领域所发挥的作用也越来越大。

在传统的傅立叶分析中，信号完全是在频域展开的，不包含任何时频的信息，其丢弃的时域信息可能对某些应用同样非常重要，所以人们对傅立叶分析进行了推广，提出了很多能表征时域和频域信息的信号分析方法，如短时傅立叶变换，Gabor变换，时频分析，小波变换等。

但短时傅立叶分析只能在一个分辨率上进行，所以对很多应用来说不够精确，存在很大的缺陷。

而小波分析则克服了短时傅立叶变换在单分辨率上的缺陷，在时域和频域都有表征信号局部信息的能力，时间窗和频率窗都可以根据信号的具体形态动态调整。

本文介绍了小波变换的基本理论，并介绍了一些常用的小波函数，然后研究了小波分析在图像处理中的应用，包括图像压缩，图像去噪，图像融合，图像增强等,本文重点在图像去噪，最后用Matlab进行了仿真[1]。

1小波分析理论
小波分析的思想最早出现在1910年Haar 提出了小波规范正交基。

1981年，Stromberg 对Haar 系进行了改造，为小波分析奠定了基础。

1986年Meyer 和Lemarie 提出了多尺度分析的思想。

后来信号分析专家Mallat 提出了多分辨分析的概念，给出了构造正交小波基的一般方法，并以多分辨分析为基础提出了著名的快速小波算法——Mallat 算法。

Mallat 算法的提出标志着小波理论获得突破性进展，从此，小波分析从理论研究走向了应用研究。

通过小波分析，可以将各种交织在一起的由不同频率组成的混合信号分解成不同频率的块信号，能够有效地解决诸如数值分析、信号分析、图像处理、量子理论、地震勘探、语音识别、计算机视觉、CT 成像、机械故障诊断等问题。

1.1 小波及小波变换
小波的核心作用是用小波及其伸缩和平移来表示函数和信号，不但具有局部化时频分析能力，而且时间分辨率和频率分辨率均可以调整。

定义：设 )()(2R L t ∈ψ,其傅立叶变换为)(ˆωψ，当)(ˆωψ满足允许条件（完全重构条件或恒等分辨条件）
⎰=R
d C ωωωψ
ψ2
)(ˆ< ∞ 时，我们称)(t ψ为一个基本小波或母小波。

将母函数)(t ψ经伸缩和平移后得
)(
1)(,a
b
t a
t b a -=
ψψ 0;,≠∈a R b a 称其为一个小波序列。

其中a 为伸缩因子，b 为平移因子。

对于任意的函数
)()(2R L t f ∈的连续小波变换为
dt a
b
t t f a
f b a W R
b a f )(
)(,),(2
/1,->==<⎰
-ψψ 其重构公式（逆变换）为 ⎰⎰∞
∞-∞
∞--=
d a d b a
b t b a W a C t f f
)(),(11)(2ψψ
把连续小波变换中的尺度参数a 和平移参数b 进行离散化：j a a 0=，
00b ka b j =，其中Z j ∈,为了方便起见，总是假设a 0>0，则得到离散小波函数
)()()(002/00
002/0
,kb t a a a b ka t a t j j j
j j k j -=-=---ψψψ 相应的离散小波变换
dt kb t a t f a t f b a W j
R
a b a f )()()(,),(002/0,->==<--⎰
ψψ 其重构公式为
)(,)()()()(,),(,,,002/0,t f t dtf kb t a t f a t f b a W b a Z
b a b a a
R
a b a f ψψψψ∑⎰
∈--><
=->==<
由于基小波)(t ψ生成的小波)(,t b a ψ在小波变换中对被分析的信号起着观测窗的作用，所以)(t ψ还应该满足一般函数的约束条件
⎰
∞
∞
-dt t )(ψ〈∞
故)(ˆωψ
是一个连续函数。

这意味着，为了满足完全重构条件式, )(ˆωψ在原点必须等于0，即
0)()0(ˆ==⎰∞
∞
-dt t ψψ
为了使信号重构的实现在数值上是稳定的，除完全重构条件外，还要求小波
)(t ψ的傅立叶变化满足下面的稳定性条件： ∑∞
∞--≤≤B A j 2
)2(ˆωψ
式中0〈A ≤B 〈∞。

1.2常用小波基介绍[3]
(1)Haar 小波
Haar 于1990年提出一种正交函数系,定义如下：
⎪⎩
⎪
⎨⎧-=011
H
ψ 其它12/12/10<≤≤≤x x 这是一种最简单的正交小波，即
0)()(=-⎰
∞
∞
-dx n x t ψψ ,2,1±±=n …
(2)Daubechies （dbN ）小波系
该小波是Daubechies 从两尺度方程系数{}k h 出发设计出来的离散正交小波。

一般简写为dbN ，N 是小波的阶数。

小波ψ和尺度函数中的支撑区为2N-1。

ϕ的消失矩为N 。

除N ＝1外(Haar 小波)，dbN 不具对称性(即非线性相位)，没有显式表达式(除N ＝1外)。

但{}k h 的传递函数的模的平方有显式表达式。

假设
∑-=+-=1
1)(N k k k N k y C y P ，其中，k N k C +-1为二项式的系数，则有
)2
(sin )2
(cos )(2
2
2
0ω
ω
ωP m N =
其中 ∑-=-=120
021)(N k ik k e h m ω
ω
(3)SymletsA （symN ）小波系
Symlets 函数系是由Daubechies 提出的近似对称的小波函数，它是对db 函数的一种改进。

Symlets 函数系通常表示为symN （N=2，3，…，8）的形式。

1.3 双尺度关系与分解关系
)(t φ与)(t ψ的两尺度关系：
)(t φ与)(t ψ的分解关系：
其中
1.4 小波进行分解与重构
两尺度函数的两尺度关系是 由两尺度关系，得序列 再由 得 则有小波函数的两尺度关系是
进一步得分解关系
由分解算法 得
由重构算法 得
⎪⎪⎩
⎪⎪⎨⎧
-=-=∑∑∞
-∞=∞
-∞=n n n n n t q t n t p t )2()()2()(φψφφ ,2,1,0)}()({)2(22±±=-+-=-∑∞
-∞
=--l n t b n t a l t n n l n l ψφφn n n n h b g a --==2121::)12()2()(-+=t t t φφφ011,1p p ==,)1(1+--=n n n p q 011,1
q q ==-)12()2()(--=t t t φφψ[]1(2)()()2t t t φφψ=+[]1(21)()()2t t t φφψ-=-⎪⎩
⎪⎨
⎧==∑∑+-+-l l k n l n k l l
k n l n k c b d c a c ,12,,12,()(),001,011,1
1,01,1,001,011,11,01,112
12
k k k k k k k k k k c a c a c c c d b c b c c c ++++++++⎧
=+=+⎪⎪⎨⎪=+=-⎪⎩)
(2,2,,1l n l k l l n l k n k q d p c c --++=∑1,0,00,00,0,0k k k k k c c p d q c d +=+=+1,1,01,01,0,0
k k k k k c c p d q c d +=+=-
2 图像处理的模型分析
2.1 图像的数学模型[4]
物体反射或投射的物质能量在空间上的分布在数学上可以表述为一能量场 E(x,y,z,λ,t )，其中x ，y ，z 表示在几何空间中点的坐标，λ为辐射波长，t 为时刻。

适当选取坐标系使取图平面垂直于z 轴，设截距为Z 0，图像可看作只是记录在平面z=z 0上的能量分布[2]，实际中这种物质能量的记录值往往用亮度值表示。

则
)],([),(y x E g y x f =
式中x ，y 为像平面中点的坐标。

图像的数学模型是一个二元函数f(x,y)，它反映了图像上点坐标f(x,y)与该点上的能量值之间的对应关系。

对于一幅灰度照片，它上面点的明暗程度需用不同的数值代表。

图像在某点处的函数值称为图像在该点的灰度或亮度。

由于f 的值是能量的记录，故其是非负有界的实数，即
A y x f ≤≤),(0
一幅实际图像的尺寸是有限的，一般定义(x,y)在某一矩形域中，即有
⎩⎨
⎧≤≤-≤≤-y y
x
x L y L L x L 2.2 图像处理的小波模型
通常情况下，图像处理可以被抽象为一个输入—输出系统，即以各种形式的算符Q 来对图像F 进行处理，算符Q 的形式决定于图像处理的目的以及图像F 的数学模型
Q Q I −−−−−−−→−T
O perator Processing mage 0
其中，T 代表图像处理算子，如去噪、锐化、分割，压缩或图像修复。

输入数据Q 0可以是一幅图也可以是图像序列，输出数据Q 是所有希望得到的图像性质。

图像处理的两个基本问题是作为输入的Q 0和算子T 的建模，它们是相对独立但又紧密相连的：算子T 的性能很大程度上取决于输入数据的模型。

为了能够有效地处理图像，首先需要知道如何从数学的角度理解和表示图像。

图像模型和它的表示方法在很大程度决定了图像的处理模型。

把图像看作像素点集合的基础上，以像素点集合为基础的数学模型可表达三种不同的模型，即随机场模型，小波模型和正则空间模型，它们分别从概率统计、小波分析和偏微分程这三方面来
研究图像处理。

小波应用图像处理的本质是一个函数逼近问题，小波可以构成Hilbert 空间的规范正交基，给出了图像的多尺度表示：
∑∑∑∈-=∈><+><=Z
n J j Z n n J n J n J n J f f x 1
1,,,,,,)(f ψψφφ
对图像的稀疏逼近为图像的高效压缩等图像处理提供了可能。

2.3 小波分析在图像处理中的应用
小波分析在图像处理中的应用主要表现在以下几个方面： ⑴图像去噪
噪声的产生是一个随机过程，噪声分量灰度值是一个随即变量，其统计特性 由概率密度函数表征，白噪声、高斯噪声、泊松噪声是三种形式常见的重要噪声； 设长度为N 的信号f n 被噪声e n 所污染，所测得的含噪数据为：
n n e f X +=n
去噪的目标是从含噪数据X 得到信号x 的一个逼近信号x ’,使得在某种误差准则估计下x ’是x 的最佳逼近。

小波去噪的基本策略是将含噪信号进行多尺度小波变换，从时域变换到小波域，然后在各尺度下尽可能提取信号的小波系数，而去除属于噪声的变换系数，然后由小波逆变换重构信号。

阀值去噪的主要理论依据是，小波变换具有很强的数据去相关性，能够使信号的能量在小波域集中在少量的大的小波系数中，而噪声却分布在整个小波域，经小波分解后，信号的小波系数的幅值要大于噪声的小波系数的幅值，于是可以采用阀值的办法把信号的小波系数保留，而使大部分噪声的小波系数减少为零，阀值去噪主要思想是对小波分解后的除了最低频以外的各层系数模对大于和小于某阈值系数分别处理，然后对处理完的小波系数再反变换重构出一幅经去噪后的图像。

在小波阈值收缩去噪法中最为重要的就是如何选择阈值和阈值函数，大多数阀值选择过程是针对一组小波系数，即根据本组小波系数的统计特性，计算出一个阈值。

阈值函数主要可以分为如下三种： 硬阈值函数
它假定幅值大于门限T的小波系数都是由信号贡献的，而小于门限T的小波系数来至于噪声。

实验表明这种取门限的方式有时并不太令人满意。

软阈值函数
硬阈值方法可以很好保留图像边缘等局部特征，但图像会出现振铃、伪吉布斯效应等视觉失真，阈值方法处理结果相对平滑得多，但是软阈值方法会造成边缘模糊等失真现象。

半软阀值函数
该方法通过选择合适的阈值T
1和T
2
，可以在软阀值方法和硬阈值方法之间达到
很好的折中。

阀值去噪算法
阀值去噪具体的算法步骤如下：
①选择合适的小波函数，对于长度为N(设N=2J)的含噪信号(常采用周期延拓方法)进行小波变换，利用小波交换的快速算法获得低分辨率L(O≤0<J)下的尺度
系数{V
L,K K=1,2,…,2L}及各分辨率下的小波系数{W
j,k
，j=L,L+1,…J-1,k=
1,2,…,2L)其中尺度系数和小波系数共N个。

②为保持信号的整体形状不变，保留所有的低频系数{V
L,K
K=1,2,…,2L}，对小波系数进行非线性阈值处理。

对每个小波系数，采用软阈值和硬阈值方法进行处理。

③进行逆小波变换。

由所有低频尺度系数，以及经由阈值处理后的小波系数作逆小波变换进行重构，得到恢复的原始信号的估计值。

⑵图像压缩小波分析用于信号与图像压缩是小波分析应用的一个重要方面。

它的特点是压缩比高，压缩速度快，压缩后能保持信号与图像的特征不变，且在传递中可抗干扰。

由于小波分析固有的时频特性，我们可以在时频两个方向对系数进行处理，这样就可以对我们感兴趣的部分提供不同的压缩精度。

⑶图形增强图像增强问题主要通过时域和频域处理两种方法来解决。

时域方法通过直接在图像点上作用算子或掩码来解决，频域方法通过修改傅立叶变换系数来解决。

这两种方法的优劣很明显，时域方法方便快速但会丢失很多点之间的相关信息，频域方法可以很详细地分离出点之间的相关,但计算量大得多。

小波分析是以上两种方法的权衡结果，傅立叶分析在所有点的分辨率都是原始图像的尺度，但我们可能不需要这么大的分辨率，而单纯的时域分析又显得太粗糙，小波分析的多尺度分析特性为用户提供了更灵活的处理方法，可以选择任意的分解层数，用近可能少的计算量得到我们满意的结果。

⑷图像融合图像融合是将同一对象的两个或更多的图像合成在一幅图像中，以便它比原来的任何一幅图像更容易为人们所理解。

这一技术可应用于多频谱图像理解以及医学图像处理等领域。

在这些场合，同一物体部件的图像往往是采用不同的成像机理得到的。

3 MATLAB仿真
下面先利用小波分析对含噪的二维图像wbard进行去噪处理，对比去噪效果，然后对去噪后的图像进行锐化处理。

再对原二维图像wbarb进行压缩处理，之后再对wbarb和woman图像进行融合处理。

Matlab程序如下：
load wbarb
%下面进行噪声的产生
init=3718025452;
rand('seed',init);
Xnoise=X+50*(rand(size(X)));
%显示原始图像及它的含噪声的图像
colormap(map);
figure(1);
subplot(2,2,1);image(wcodemat(X,192));
title('原始图像X')
subplot(2,2,2);image(wcodemat(Xnoise,192));
title('含噪声的图像Xnoise');
%用sym5小波对图像信号进行二层的小波分解
[c,s]=wavedec2(X,2,'sym5');
%下面进行图像的去噪处理
%使用ddencmp函数来计算去噪的默认阈值
%使用wdencmp函数来实现图像的压缩
[thr,sorh,keepapp]=ddencmp('den','wv',Xnoise);
[Xdenoise,cxc,lxc,perf0,perfl2]=wdencmp('gbl',c,s,'sym5',2,thr,sorh,keepapp);
%显示去噪后的图像
subplot(223);image(Xdenoise);
title('去噪后的图像');
%对去噪后的图像做锐化处理
blur2=Xdenoise;
%对图像做2层的二维小波分解
[c,l]=wavedec2(Xdenoise,2,'db3');
csize=size(c);
%对低频系数进行放大处理，并抑制高频系数
for i=1:csize(2);
if(abs(c(i))<300)
c(i)=c(i)*2;
else
c(i)=c(i)/2;
end
end
%通过处理后的小波系数重构图像
blur2=waverec2(c,l,'db3');
subplot(224);image(wcodemat(blur2,192));colormap(gray(256));title('锐化去噪后的图像');
%图像压缩处理
figure(2);
subplot(221);image(X);colormap(map)
title('原始图像');
disp('压缩前图像X的大小：');
whos('X')
%对图像用bior3.7小波进行2层小波分解
[c,s]=wavedec2(X,2,'bior3.7');
%提取小波分解结构中第一层低频系数和高频系数
ca1=appcoef2(c,s,'bior3.7',1);
ch1=detcoef2('h',c,s,1);
cv1=detcoef2('v',c,s,1);
cd1=detcoef2('d',c,s,1);
%分别对各频率成分进行重构
a1=wrcoef2('a',c,s,'bior3.7',1);
h1=wrcoef2('h',c,s,'bior3.7',1);
v1=wrcoef2('v',c,s,'bior3.7',1);
d1=wrcoef2('d',c,s,'bior3.7',1);
c1=[a1,h1;v1,d1];
%显示分解后各频率成分的信息
subplot(222);image(c1);
title('分解后低频和高频信息');
%保留小波分解第一层低频信息，进行图像的压缩
%第一层的低频信息即为ca1，显示第一层的低频信息%首先对第一层信息进行量化编码
ca1=appcoef2(c,s,'bior3.7',1);
ca1=wcodemat(ca1,440,'mat',0);
%改变图像的高度
ca1=0.5*ca1;
subplot(223);image(ca1);colormap(map);
title('第一次压缩');
disp('第一次压缩图像的大小为：');
whos('ca1')
%保留小波分解第二层低频信息，进行图像的压缩，此时压缩比更大%第二层的低频信息即为ca2，显示第二层的低频信息
ca2=appcoef2(c,s,'bior3.7',2);
%首先对第二层信息进行量化编码
ca2=wcodemat(ca2,440,'mat',0);
%改变图像的高度
ca2=0.25*ca2;
subplot(224);image(ca2);colormap(map);
axis square
title('第二次压缩');
disp('第二次压缩图像的大小为：');
whos('ca2')
%图像融合
figure(3);
X2=X;map2=map;
subplot(221);image(X2);colormap(map2);
title('wbarb');
load woman;
X1=X;map1=map;
%画出原始图像
subplot(222);image(X1);colormap(map1);
title('woman');
%用小波函数sym4对X1进行2层小波分解
[c1,s1]=wavedec2(X1,2,'sym4');
%对分解系数进行处理以突出轮廓部分，弱化细节部分
sizec1=size(c1);
for i=1:sizec1(2)
c1(i)=1.2*c1(i); end
%用小波函数sym4对X2进行2层小波分解 [c2,s2]=wavedec2(X2,2,'sym4'); %下面进行小波变换域的图像融合 c=c1+c2; %减小图像亮度 c=0.5*c;
%对融合的系数进行重构 xx=waverec2(c,s1,'sym4'); %画出融合后的图像 subplot(223);image(xx);
title('wbarb 和woman 的融合图像'); 去噪和锐化仿真结果：
原始图像X
50
100
150
200
250
50100150200
250
含噪声的图像X noise
50
100
150
200
250
50100150200
250
去噪后的图像
50
100
150
200
250
50100150200
250
锐化去噪后的图像
50100150200250
50100150200
250
图像压缩仿真结果：
原始图像
50
100
150
200
250
50100150200
250
分解后低频和高频信息
100200300400500
100200300400
500
第一次压缩
2040
60
80100120
20406080100
120
第二次压缩
204060
2040
60
Matlab 输出：
压缩前图像X 的大小：
Name Size Bytes Class X 256x256 524288 double array Grand total is 65536 elements using 524288 bytes 第一次压缩图像的大小为：
Name Size Bytes Class ca1 135x135 145800 double array Grand total is 18225 elements using 145800 bytes 第二次压缩图像的大小为：
Name Size Bytes Class ca2 75x75 45000 double array Grand total is 5625 elements using 45000 bytes 图像融合仿真结果：
wbarb
50
100
150
200
250
50100150200
250
woman
50
100150200250
50100150200
250
wbarb 和woman 的融合图像50
100
150
200
250
50100150200
250
3结束语
本文详细讨论了小波在图像处理领域的应用，并用Matlab 进行了仿真，Matlab 语言对数字图像进行处理时具有编程简单、处理速度快的特点，结果显示用小波变换理论进行图像处理可以取得较好的效果。

在本文中只是实现了图像处理的基本功能,对于处理效果还有待改善,在后续进一步研究中可以选用更合适的小波基和算法进行处理,以使处理效果更佳,
参考文献
[1] 李朝晖、张弘编.数字图像处理及应用[M].北京:机械工业出版社，2004（06）. [2] 张亶、陈刚编.基于偏微分方程的图像处理[M].北京:高等教育出版社，
2004(01).
[3] 崔锦泰(著)，程正兴(译)．小波分析导论[M]．西安：西安交通大学出版社，
1995
[4] 龚声蓉、刘纯平、王强编.数字图像的处理与分析[M]．北京：清华大学出版
社2006．。