本科毕业设计__基于matlab的小波分析在图像处理中的应用
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
基于Matlab 的小波分析在图像处理中的应用
摘要:本文先介绍了小波分析得基本理论,包括连续小波变换、离散小波变换和小波包分析。
小波变换具有时频局部化的特点,因此不但能对图像提供较精确的时域定位,也能提供较精确的频域定位。
经过小波变换的图像具有频谱划、方向选择、多分辨率分析和天然塔式数据结构特点。
基于小波变换这些特性,讨论了MATLAB 语言环境下图像压缩,图像去噪,图像融合,图像分解,图像增强的基本方法。
关键词:小波分析;图像压缩;图像去噪;图像融合;图像分解;图像增强
1 引言
小波分析诞生于20世纪80年代, 被认为是调和分析即现代Fourier 分析发展的一个崭新阶段。
众多高新技术以数学为基础,而小波分析被誉为“数学显微镜”,这就决定了它在高科技研究领域重要的地位。
目前, 它在模式识别、图像处理、语音处理、故障诊断、地球物理勘探、分形理论、空气动力学与流体力学上的应用都得到了广泛深入的研究,甚至在金融、证券、股票等社会科学方面都有小波分析的应用研究。
在传统的傅立叶分析中,信号完全是在频域展开的,不包含任何时频的信息,这对于某些应用来说是很恰当的,因为信号的频率的信息对其是非常重要的。
但其丢弃的时域信息可能对某些应用同样非常重要,所以人们对傅立叶分析进行了推广,提出了很多能表征时域和频域信息的信号分析方法,如短时傅立叶变换,Gabor 变换,时频分析,小波变换等。
其中短时傅立叶变换是在傅立叶分析基础上引入时域信息的最初尝试,其基本假定在于在一定的时间窗内信号是平稳的,那么通过分割时间窗,在每个时间窗内把信号展开到频域就可以获得局部的频域信息,但是它的时域区分度只能依赖于大小不变的时间窗,对某些瞬态信号来说还是粒度太大。
换言之,短时傅立叶分析只能在一个分辨率上进行。
所以对很多应用来说不够精确,存在很大的缺陷。
而小波分析则克服了短时傅立叶变换在单分辨率上的缺陷,具有多分辨率分析的特点,在时域和频域都有表征信号局部信息的能力,时间窗和频率窗都可以根据信号的具体形态动态调整,在一般情况下,在低频部分(信号较平稳)可以采用较低的时间分辨率,而提高频率的分辨率,在高频情况下(频率变化不大)可以用较低的频率分辨率来换取精确的时间定位。
本文介绍了小波变换的基本理论,并介绍了一些常用的小波函数,它们的主要性质包括紧支集长度、滤波器长度、对称性、消失矩等,都做了简要的说明。
然后研究了小波分析在图像处理中的应用,包括图像压缩,图像去噪,图像融合,图像分解,图像增强等。
2 小波分析的基本理论
2.1 连续小波变换
定义:设)()(2R L t ∈ψ,其傅立叶变换为)(ˆωψ
,当)(ˆωψ满足允许条件(完全重构条件或恒等分辨条件)
⎰=R
d C ωωωψ
ψ2
)(ˆ< ∞ (1)
时,我们称)(t ψ为一个基本小波或母小波。
将母函数)(t ψ经伸缩和平移后得 )(
1)(,a
b
t a
t b a -=
ψψ 0;,≠∈a R b a (2) 称其为一个小波序列。
其中a 为伸缩因子,b 为平移因子。
对于任意的函数)()(2R L t f ∈的连续小波变换为
dt a
b
t t f a f b a W R
b a f )(
)(,),(2
/1,->==<⎰
-ψψ (3) 其重构公式(逆变换)为
⎰⎰∞
∞-∞
∞--=
dadb a
b t b a W a C t f f
)(),(11
)(2ψψ
(4) 由于基小波)(t ψ生成的小波)(,t b a ψ在小波变换中对被分析的信号起着观测窗的作用,所以
)(t ψ还应该满足一般函数的约束条件
⎰
∞
∞
-dt t )(ψ〈∞ (5)
故)(ˆωψ
是一个连续函数。
这意味着,为了满足完全重构条件式,)(ˆωψ在原点必须等于0,即
0)()0(ˆ==⎰∞
∞-dt t ψψ
(6) 为了使信号重构的实现在数值上是稳定的,处理完全重构条件外,还要求小波)(t ψ的傅立叶变化满足下面的稳定性条件:
∑∞
∞--≤≤B A j 2
)2(ˆωψ
(7) 式中0〈A ≤B 〈∞。
2.2 离散小波变换
在实际运用中,尤其是在计算机上实现时,连续小波必须加以离散化。
因此,有必要讨论连续小波)(,t b a ψ和连续小波变换),(b a W f 的离散化。
需要强调指出的是,这一离散化都是针对连续的尺度参数a 和连续平移参数b 的,而不是针对时间变量t 的。
这一点与我们以前习惯的时间离散化不同。
在连续小波中,考虑函数:
)()(2
/1,a
b t a t b a -=-ψψ
这里R b ∈,+∈R a ,且0≠a ,ψ是容许的,为方便起见,在离散化中,总限制a 只
取正值,这样相容性条件就变为 ∞<=⎰
∞ωω
ωψ
ψd C 0
)(ˆ (8) 通常,把连续小波变换中尺度参数a 和平移参数b 的离散公式分别取作j a a 0=,
0b ka b j =,这里Z j ∈,扩展步长10≠a 是固定值,为方便起见,总是假定10>a (由于m
可取正也可取负,所以这个假定无关紧要)。
所以对应的离散小波函数)(,t k j ψ即可写作
)()()(002/00
002
/0
,kb t a a a b ka t a t j j j
j j k j -=-=---ψψψ (9) 而离散化小波变换系数则可表示为
>=<=⎰
∞∞
-k j k j k j f dt t t f C ,*,,,)()(ψψ (10)
其重构公式为
∑∑∞
∞-∞
∞-=)()(,,t C C t f k j k j ψ (11)
C 是一个与信号无关的常数。
然而,怎样选择0a 和0b ,才能够保证重构信号的精度呢?显然,网格点应尽可能密(即0a 和0b 尽可能小),因为如果网格点越稀疏,使用的小波函数)(,t k j ψ和离散小波系数k j C ,就越少,信号重构的精确度也就会越低。
2.3 小波包分析
短时傅立叶变换对信号的频带划分是线性等间隔的。
多分辨分析可以对信号进行有效
的时频分解,但由于其尺度是按二进制变化的,所以在高频频段其频率分辨率较差,而在低频频段其时间分辨率较差,即对信号的频带进行指数等间隔划分(具有等Q 结构)。
小波包分析能够为信号提供一种更精细的分析方法,它将频带进行多层次划分,对多分辨率分析没有细分的高频部分进一步分解,并能够根据被分析信号的特征,自适应地选择相应频带,使之与信号频谱相匹配,从而提高了时-频分辨率,因此小波包具有更广泛的应用价值。
关于小波包分析的理解,我们这里以一个三层的分解进行说明,其小波包分解树如图
图1 小波包分解树
图1中,A 表示低频,D 表示高频,末尾的序号数表示小波分解的层树(也即尺度数)。
分解具有关系:
S=AAA3+DAA3+ADA3+DDA3+AAD3+DAA3+ADD3+DDD3
3 常用小波基介绍
(1)Haar 小波
Haar 于1990年提出一种正交函数系,定义如下:
⎪⎩
⎪
⎨⎧-=011
H ψ 其它12/12/10<≤≤≤x x (12)
这是一种最简单的正交小波,即
0)()(=-⎰
∞
∞
-dx n x t ψψ ,2,1±±=n …
(2)Daubechies (dbN )小波系
该小波是Daubechies 从两尺度方程系数{}k h 出发设计出来的离散正交小波。
一般简写为dbN ,N 是小波的阶数。
小波ψ和尺度函数吁中的支撑区为2N-1。
ϕ的消失矩为N 。
除N =1外(Haar 小波),dbN 不具对称性〔即非线性相位〕;dbN 没有显式表达式(除N =1外)。
但{}k h 的传递函数的模的平方有显式表达式。
假设∑-=+-=1
01)(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 = (13)
其中 ∑-=-=
120
02
1
)(N k ik k
e
h m ω
ω
(3)Biorthogonal (biorNr.Nd )小波系
Biorthogonal 函数系的主要特征体现在具有线性相位性,它主要应用在信号与图像的重构中。
通常的用法是采用一个函数进行分解,用另外一个小波函数进行重构。
Biorthogonal 函数系通常表示为biorNr.Nd 的形式:
Nr=1 Nd=1,3,5 Nr=2 Nd=2,4,6,8 Nr=3 Nd=1,3,5,7,9 Nr=4 Nd=4 Nr=5 Nd=5 Nr=6 Nd=8
其中,r 表示重构,d 表示分解。
(4)Coiflet (coifN )小波系
coiflet 函数也是由Daubechies 构造的一个小波函数,它具有coifN (N=1,2,3,4,5)这一系列,coiflet 具有比dbN 更好的对称性。
从支撑长度的角度看,coifN 具有和db3N 及
sym3N 相同的支撑长度;从消失矩的数目来看,coifN 具有和db2N 及sym2N 相同的消失矩数目。
(5)SymletsA (symN )小波系
Symlets 函数系是由Daubechies 提出的近似对称的小波函数,它是对db 函数的一种改进。
Symlets 函数系通常表示为symN (N=2,3,…,8)的形式。
(6)Morlet (morl )小波
Morlet 函数定义为x Ce x x
5cos )(2
/2
-=ψ,它的尺度函数不存在,且不具有正交性。
(7)Mexican Hat (mexh )小波
Mexican Hat 函数为
2/24
/12)1(3
2)(x e x x ---=
ψπ (14) 它是Gauss 函数的二阶导数,因为它像墨西哥帽的截面,所以有时称这个函数为墨西哥帽函数。
墨西哥帽函数在时间域与频率域都有很好的局部化,并且满足
0)(=⎰
∞
∞
-dx x ψ
由于它的尺度函数不存在,所以不具有正交性。
(8)Meyer 小波
Meyer 小波函数ψ和尺度函数ϕ都是在频率域中进行定义的,是具有紧支撑的正交小波。
⎪⎪⎪
⎩
⎪⎪
⎪⎨⎧--=ψ--0))123(2cos()2())123(2sin()2()(ˆ2/2/12
/2/1ωπυππωπυππωωωj j e
e ]
3
8,32[3
8343432ππωπ
ωππωπ∉≤
≤≤≤ (15) 其中,)(a υ为构造Meyer 小波的辅助函数,且有
⎪⎪
⎩⎪⎪
⎨⎧-=--0))123(2cos()2()2()(ˆ2/12/1ωπυπππωφ 3
434323
2πωπωππω>≤≤≤
(16) 4 小波分析在图像处理中的应用
4.1 小波分析用于图像压缩
4.1.1 基于小波变换的图像局部压缩
基于离散余弦变换的图像压缩算法,其基本思想是在频域对信号进行分解,驱除信号
点之间的相关性,并找出重要系数,滤掉次要系数,以达到压缩的效果,但该方法在处理过程中并不能提供时域的信息,在我们比较关心时域特性的时候显得无能为力。
但是这种应用的需求是很广泛的,比如遥感测控图像,要求在整幅图像有很高压缩比的同时,对热点部分的图像要有较高的分辨率,例如医疗图像,需要对某个局部的细节部分有很高的分辨率,单纯的频域分析的方法显然不能达到这个要求,虽然可以通过对图像进行分快分解,然后对每块作用不同的阈值或掩码来达到这个要求,但分块大小相对固定,有失灵活。
在这个方面,小波分析就优越的多,由于小波分析固有的时频特性,我们可以在时频两个方向对系数进行处理,这样就可以对我们感兴趣的部分提供不同的压缩精度。
下面我们利用小波变化的时频局部化特性,举一个局部压缩的例子,可以通过这个例子看出小波变换在应用这类问题上的优越性。
load wbarb
%使用sym4小波对信号进行一层小波分解
[ca1,ch1,cv1,cd1]=dwt2(X,'sym4');
codca1=wcodemat(ca1,192);
codch1=wcodemat(ch1,192);
codcv1=wcodemat(cv1,192);
codcd1=wcodemat(cd1,192);
%将四个系数图像组合为一个图像
codx=[codca1,codch1,codcv1,codcd1]
%复制原图像的小波系数
rca1=ca1;
rch1=ch1;
rcv1=cv1;
rcd1=cd1;
%将三个细节系数的中部置零
rch1(33:97,33:97)=zeros(65,65);
rcv1(33:97,33:97)=zeros(65,65);
rcd1(33:97,33:97)=zeros(65,65);
codrca1=wcodemat(rca1,192);
codrch1=wcodemat(rch1,192);
codrcv1=wcodemat(rcv1,192);
codrcd1=wcodemat(rcd1,192);
%将处理后的系数图像组合为一个图像
codrx=[codrca1,codrch1,codrcv1,codrcd1]
%重建处理后的系数
rx=idwt2(rca1,rch1,rcv1,rcd1,'sym4');
subplot(221);image(wcodemat(X,192)),colormap(map);title('原始图像');
subplot(222);image(codx),colormap(map);title('一层分解后各层系数图像');
subplot(223);image(wcodemat(rx,192)),colormap(map);title('压缩图像');
subplot(224);image(codrx),colormap(map);title('处理后各层系数图像');
%求压缩信号的能量成分
per=norm(rx)/norm(X)
per =1.0000
%求压缩信号与原信号的标准差
err =586.4979
原始图像
50
100
150
200
250
50100150200250
一层分解后各层系数图像
100
200
300
400
500
20406080100120
压缩图
像
50
100
150
200
250
50100150200250
处理后各层系数图像
100
200
300
400
500
20406080100120
图2 利用小波变换的局部压缩图像
从图1可以看出,小波域的系数表示的是原图像各频率段的细节信息,并且给我们提供了一种位移相关的信息表述方式,我们可以通过对局部细节系数处理来达到局部压缩的效果。
在本例中,我们把图像中部的细节系数都置零,从压缩图像中可以很明显地看出只有中间部分变得模糊(比如在原图中很清晰的围巾的条纹不能分辨),而其他部分的细节信息仍然可以分辨的很清楚。
最后需要说明的是本例只是为了演示小波分析应用在图像局部压缩的方法,在实际的应用中,可能不会只做一层变换,而且作用阈值的方式可能也不会是将局部细节系数全部清除,更一般的情况是在N 层变换中通过选择零系数比例或能量保留成分作用不同的阈值,实现分片的局部压缩。
而且,作用的阈值可以是方向相关的,即在三个不同方向的细节系数上作用不同的阈值。
4.1.2 利用二维小波分析进行图像压缩
二维小波分析用于图像压缩是小波分析应用的一个重要方面。
它的特点是压缩比高,压缩速度快,压缩后能保持图像的特征基本不变,且在传递过程中可以抗干扰。
小波分析用于图像压缩具有明显的优点。
下面给出一个图像信号(即一个二维信号,文件名为wbarb.mat ),利用二维小波分析对图像进行压缩。
一个图像作小波分解后,可得到一系列不同分辨率的子图像,不同分辨率的子图像对应的频率是不相同的。
高分辨率(即高频)子图像上大部分点的数值都接近于0,越是高频这种现象越明显。
对一个图像来说,表现一个图像最主要的部分是低频部分,所以一个最简单的压缩方法是利用小波分解,去掉图像的高频部分而只保留低频部分。
图像压缩可按如下程序进行处理。
load wbarb;
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);
axis square
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);
axis square
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('第二次压缩图像的大小为:');
输出结果如下所示: 压缩前图像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
图像对比如图所示。
可以看出,第一次压缩提取的是原始图像中小波分解第一层的低频信息,此时压缩效果较好,压缩比较小(约为1/3):第二次压缩是提取第一层分解低频部分的低频部分(即小波分解第二层的低频部分),其压缩比较大(约为1/12),压缩效果在视觉上也基本过的去。
这是一种最简单的压缩方法,只保留原始图像中低频信息,不经过其他处理即可获得较好的压缩效果。
在上面的例子中,我们还可以只提取小波分解第3、4、…层的低频信息。
从理论上说,我们可以获得任意压缩比的压缩图像。
原始图像
5010015020025050100150200250
分解后低频和高频信息
100200300400500100200300400500
第一次压缩
20406080100120
20406080100
120
第二次压缩
20
40
60
204060
图3 利用二维小波分析进行图像压缩
下面再给出用wdenemp 函数对一个图像(文件名tire.mat )进行压缩的程序。
%装入一个二维信号
load tire; %显示图像
subplot(221);image(X);colormap(map) title('原始图像');
%下面进行图像压缩
%对图像用db3小波进行2层小波分解 [c,s]=wavedec2(X,2,'db3');
%使用wavedec2函数来实现图像的压缩 [thr,sorh,keepapp]=ddencmp('cmp','wv',X); %输入参数中选择了全局阈值选项‘gbl ’,用来对所有高频系数进行相同的阈值量化处理 [Xcomp,cxc,lxc,perf0,perfl2]=wdencmp('gbl',c,s,'db3',2,thr,sorh,keepapp); %将压缩后的图像与原始图像相比较,并显示出来 subplot(222);image(Xcomp);colormap(map) title('压缩图像'); axis square
disp('小波分解系数中置0的系数个数百分比:'); perf0
disp('压缩后图像剩余能量百分比:'); perfl2
输出结果如下所示:
小波分解系数中置0的系数个数百分比: perf0 =49.1935
压缩后图像剩余能量百分比: perfl2 =99.9928 图像对比如图所示:
原始图像
50100150200
50100150200
压缩图像
50100150200
50100150200
图4 利用二维小波分析对图像进行压缩
利用二维小波变换进行图像压缩时,小波变换将图像从空间域变换到时间域,它的作用与以前在图像压缩中所用到的离散余弦(DCT )、傅立叶变换(FFT )等的作用类似。
但是要很好的进行图像的压缩,需要综合的利用多种其他技术,特别是数据的编码与解码算法等,所以利用小波分析进行图像压缩通常需要利用小波分析和许多其他相关技术共同完成。
4.1.3 基于小波包变换的图像压缩
小波分析之所以在信号处理中有着强大的功能,是基于其分离信息的思想,分离到各个小波域的信息除了与其他小波域的关联,使得处理的时候更为灵活。
全局阈值化方法作用的信息粒度太大,不够精细,所以很难同时获得高的压缩比和能量保留成分,在作用的
分层阈值以后,性能明显提高,因为分层阈值更能体现信号固有的时频局部特性。
但是小波分解仍然不够灵活,分解出来的小波树只有一种模式,不能完全地体现时频局部化信息。
而压缩的核心思想既是尽可能去处各小波域系数之间的信息关联,最大限度体现时频局部化的信息,因此,实际的压缩算法多采用小波包算法,而小波树的确定则是根据不同的信息论准则,以达到分解系数表达的信息密度最高。
下面我通过一个例子来说明小波包分析在图像压缩中的应用,并给出性能参数以便于同基于小波分析的压缩进行比较。
load julia
%求颜色索引表长度 nbc=size(map,1);
%得到信号的阈值,保留层数,小波树优化标准 [thr,sorh,keepapp,crit]=ddencmp('cmp','wp',X) %通过以上得到的参数对信号进行压缩
[xd,treed,perf0,perfl2]=wpdencmp(X,sorh,4,'sym4',crit,thr*2,keepapp); %更改索引表为pink 索引表 colormap(pink(nbc));
subplot(121);image(wcodemat(X,nbc));title('原始图像');
subplot(122);image(wcodemat(xd,nbc));title('全局阈值化压缩图像'); xlabel(['能量成分',num2str(perfl2),'%','零系数成分',num2str(perf0),'%']); plot(treed);
得到的压缩结果如图所示
原始图像
100200300
50100150200250全局阈值化压缩图像
能量成分99.7346%零系数成分85.81%
100200300
50
100
150
200
250
图5 基于小波包分析的图像压缩
压缩过程中使用的最优小波树如图6所示
图6 最优小波树
这两个命令是Matlab小波工具箱提供的自动获取阈值和自动使用小波包压缩的命令,后者将分解阈值化和重建综合起来。
在将小波包用于信号压缩的过程中,ddencmp命令返回的最优小波树标准都是阈值化标准。
根据这个标准确定的最优小波树可以使得压缩过程的零系数成分最高,并且自动降低计算量。
最后需要说明的一点,对高频成分很多的图像,小波包的分解细节信息的特点尤其能发挥其优势。
图像压缩是应用非常广泛的一类问题,所以其机器实现效率是至关重要的,在实际的应用中,如JPEG2000,一般不采用通常的mallat算法做小波分解,而是应用特定的双正交小波,利用其滤波器分布规则的特性,用移位操作来实现滤波操作。
4.2 小波分析用于图像去噪
对二维图像信号的去噪方法同样适用于一维信号,尤其是对于几何图像更适合。
二维模型可以表述为
s(i,j)=f( i,j)+δ·e(i,j) i,j=0,1,…,m-1
其中,e是标准偏差不变的高斯白噪声。
二维信号用二维小波分析的去噪步骤有3步:(1)二维信号的小波分解。
选择一个小波和小波分解的层次N,然后计算信号s到第N层的分解。
(2)对高频系数进行阈值量化。
对于从1到N的每一层,选择一个阈值,并对这一层的高频系数进行软阈值量化处理。
(3)二维小波的重构。
根据小波分解的第N层的低频系数和经过修改的从第一层到第N层的各层高频系数计算二维信号的小波重构。
在这3个步骤中,重点是如何选取阈值和阈值的量化
下面给出一个二维信号(文件名为detfinger.mat),并利用小波分析对信号进行去噪处理。
Matlab的去噪函数有ddencmp,wdencmp等,其去噪过程可以按照如下程序进行。
load tire
%下面进行早声的产生
init=3718025452; rand('seed',init);
Xnoise=X+18*(rand(size(X)));
%显示原始图像及它的含噪声的图像 colormap(map);
subplot(2,2,1);image(wcodemat(X,192)); title('原始图像X') axis square
subplot(2,2,2);image(wcodemat(X,192)); title('含噪声的图像Xnoise'); axis square
%用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('去噪后的图像'); axis square
输出结果从图中3个图像的比较可以看出,Matlab 中的ddencmp 和wdencmp 函数可以有效地进行去噪处理。
原始图像X
50100150200
50100150
200
含噪声的图像X noise 50100150200
50100150
200
去噪后的图像
5010015020050100150
200
图7 去噪一
再给定一个有较大白噪声的图像。
由于图像所含的噪声主要是白噪声,而且主要集中在图像的高频部分,所以我们可以通过全部滤掉图像中的高频部分实现图像的去噪。
具体去噪过程可按照如下程序进行。
load wmandril; %画出原始图像
subplot(221);image(X);colormap(map); title('原始图像'); axis square
%产生含噪图像
init=2055615866;randn('seed',init) x=X+38*randn(size(X)); %画出含噪图像
subplot(222);image(x);colormap(map); title('含噪声图像'); axis square;
%下面进行图像的去噪处理
%用小波函数sym4对x 进行2层小波分解 [c,s]=wavedec2(x,2,'sym4');
%提取小波分解中第一层的低频图像,即实现了低通滤波去噪 a1=wrcoef2('a',c,s,'sym4'); %画出去噪后的图像 subplot(223);image(a1); title('第一次去噪图像'); axis square;
%提取小波分解中第二层的低频图像,即实现了低通滤波去噪 %相当于把第一层的低频图像经过再一次的低频滤波处理 a2=wrcoef2('a',c,s,'sym4',2); %画出去噪后的图像
subplot(224);image(a2);title('第二次去噪图像'); axis square;
输出结果如图:
原始图像
20406080100120
20406080100120
含噪声图像
20406080100120
20406080100120
第一次去噪图像20406080100120
20406080100120
第二次去噪图像20406080100120
20406080100120
图8 去噪二
从上面的输出结果可以看出,第一次去噪已经滤去了大部分的高频噪声,但从去噪图像与原始图像相比可以看书,第一次去噪后的图像中还是含有不少的高频噪声;第二次去噪是在第一次去噪的基础上,再次滤去其中的高频噪声。
从去噪的结果可以看出,它具有
较好的去噪效果。
下面再给出定一个喊有较少噪声的facets.mat 图像。
由于原始图像中只喊有较少的高频噪声,如果按照上一个例子把高频噪声全部滤掉的方法将损坏图像中固有的高频有用信号。
因此这幅图像适合采用小波分解系数阈值量化方法进行去噪处理。
load facets;
%画出原始图像
subplot(221);image(X);colormap(map); title('原始图像'); axis square
%产生含噪声图像
init=2055615866;randn('seed',init) x=X+10*randn(size(X)); %画出含噪声图像
subplot(222);image(X);colormap(map); title('含噪声图像'); axis square
%下面进行图像的去噪处理
%用小波画数coif3对x 进行2层小波分解 [c,s]=wavedec2(x,2,'coif3');
%提取小波分解中第一层的低频图像,即实现了低通滤波去噪 %设置尺度向量n n=[1,2]
设置阈值向量p p=[10.12,23.28];
%对三个方向高频系数进行阈值处理 nc=wthcoef2('h',c,s,n,p,'s'); nc=wthcoef2('v',c,s,n,p,'s'); nc=wthcoef2('d',c,s,n,p,'s');
%对新的小波分解结构[nc ,s]进行重构 xx=waverec2(nc,s,'coif3'); %画出重构后图像的波形
subplot(223);image(X);colormap(map); title('去噪后的图像'); axis square
输出结果如图
原始图像
50100150200250
50100150200250
含噪声图像
50100150200250
50100150200250
去噪后的图像
50100150200250
50100150200250
图9 去噪三
二维信号在应用中一般表现为图像信号,二维信号在小波域中的降噪方法的基本思想与一维情况一样,在阈值选择上,可以使用统一的全局阈值,有可以分作三个方向,分别是水平方向、竖直方向和对角方向,这样就可以把在所有方向的噪声分离出来,通过作用阈值抑制其成分。
4.3 小波分析用于图像增强
4.3.1 图像增强问题描述
小波分析在二维信号(图像)处理方面的优点主要体现在其时频分析特性,前面介绍了一些基于这种特性的一些应用的实例,但对二维信号小波系数的处理方法只介绍了阈值化方法一种,下面我将介绍一下以前在一维信号中用到的抑制系数的方法,这种方法在图像处理领域主要应用于图像增强。
在图像处理领域,图像增强问题主要通过时域(沿用信号处理的说法,空域可能对图像更适合)和频域处理两种方法来解决。
时域方法通过直接在图像点上作用算子或掩码来解决,频域方法通过修改傅立叶变换系数来解决。
这两种方法的优劣很明显,时域方法方便快速但会丢失很多点之间的相关信息,频域方法可以很详细地分离出点之间的相关,但需要做两次数量级为nlogn的傅立叶变换和逆变换的操作,计算量大得多。
小波分析是以上两种方法的权衡结果,建立在如下的认识基础上,傅立叶分析的在所有点的分辨率都是原始图像的尺度,但对于问题本身的要求,我们可能不需要这么大的分辨率,而单纯的时域分析又显得太粗糙,小波分析的多尺度分析特性为用户提供了更灵活的处理方法。
可以选择任意的分解层数,用进可能少的计算量得到我们满意的结果。
给定一个wmandril.mat图像信号。
由于图像经二维小波分解后,图像的轮廓主要体现在低频部分,细节部分体现在高频部分,因此可以通过对低频分解系数进行增强处理,对高频分解系数进行衰减处理,从而达到图像增强的效果。
具体程序清单如下:load wmandril
subplot(121);image(X);colormap(map);title('原始图像');
axis square
%下面进行图像的增强处理
%用小波函数sym4对X进行2层小波分解
[c,s]=wavedec2(X,2,'sym4');
sizec=size(c);
%对分解系数进行处理以突出轮廓部分,弱化细节部分
for i=1:sizec(2)
if(c(i)>350)
c(i)=2*c(i);
else
c(i)=0.5*c(i);
end
end
%下面对处理后的系数进行重构
xx=waverec2(c,s,'sym4');
%画出重构后的图像
subplot(122);image(xx);colormap(map);title('增强图像');
axis square。