第16章 小波分析用于图像处理其他领域

在前面的章节中，我们重点介绍了小波分析在图像压缩、图像去噪、图像增强(包括图像钝化和锐化、图像平滑等)方面的应用，这些应用在图像处理中非常常见。
本章将介绍小波分析用于图像处理的其他方面，比如图像融合、图像分解及图像特征提取等。随着小波理论与技术的日渐成熟，这些领域也渐渐地引进了小波分析方法，从而拓展了图像处理的技术道路，也改进了传统图像处理的效果。然而，任何一项理论与技术在现今都不能够称为“完美”，只能是相对而言。
本章的结构安排为：首先介绍图像融合所涉及的小波技术，如小波变换法、小波包法、多小波及小波框架等；随后介绍了图像分解及图像特征提取等。在每节都例举了这些小波应用实例，希望对此领域感兴趣的广大读者和技术爱好者可以从中找到快乐，愉悦心情。
学习目标：
(1)熟练掌握小波图像增强包含的内容和方法
(2)熟练掌握小波图像分解
(3)熟练掌握小波图像特征提取
16.1 图像融合
图像融合是将同一场景的多幅图像的互补信息合并成一幅新图像，以便更好地对场景进行监视和侦察。
16.1.1 小波分析用于图像融合的方法
小波分析具有多分辨等特点，可以有效地将特征明显、分辨率高的图像融合在一起，得到比任何一幅源图像效果都好的图像。
图像融合是信息融合的一个重要分支，在自动目标识别、遥感、机器人视觉、智能系统、医学图像处理以及军事领域有着广泛的应用。对融合后的图像的要求为：
•充分利用各源图像互补信息；
•更适合人的视觉感受；
•适合进一步分析的需要；
•统一编码，压缩数据量，便于通信等。
图像融合分几个方面：信号级融合、像素级融合、特征级融合、符号级融合。像素级融合是将各源图像中对应的像素进行融合处理，它是后两级的基础。近年来，随着小波理论及其应用的发展，已将小波多尺度分解用于像素级图像融合。小波变换的固有特性使其在图像处理中有如下优点：
•完善的重构能力，保证信号在分解过程中没有信息损失和冗余信息；
•把图像分解成逼近图像和细节图像之和，分别代表了图像的不同结构，因此容易提取原始图像的结构信息和细节信息；
•具有快速算法(Mallat算法)，它在小波变换中的作用相当于FFT 在傅立叶变换中的作用，为小波变换应用提供了必要的手段；
•二维小波分析提供了与人类视觉系统方向相吻合的选择性图像。
因此，基于小波分解的图像融合越来越受到重视。小波变换的目的是将原始图像分别分解到一系列频率通道中，利用分解后的金字塔或树结构，对不同分解层、不同频带进行融合处理，可有效地将来自不同图像的细节融合在一起。Campbell和Robson的实验表明，人的视网膜图像是在不同频带上分别以不同算子进行融合的，基于小波分解的图像融合也是以同样方式进行的，因而可获得与人的视觉特性更为接近的融合效果。
若对二维图像进行N层的小波分解，最终将有(3N+1)个不同频带，其中包含3N个高频带和1个低频带。
融合基本步骤如下：
(1)对每一源图像分别进行小波分解，建立图像的小波金字塔分解；
(2)对各分解层分别进行融合处理，采用不同的融合算子对各分解层的不同频率分量进行融合处理，最终得到融合后的小波金字塔；
(3)对融合后所得的小波金字塔进行小波逆变换，得到重构图像。
16.1.2 融合规则和融合算子
针对不同类型的图像有以下融合规则。
(1)取系数绝对值较大法。该法适合高频成分较丰富，亮度、对比度较高的源图像，否则在融合图像中只保留一幅源图像的特征，其他的特征被覆盖。融合图像中基本保留源图像的特征，图像对比度与源图像基本相同。小波变换的实际作用是对信号解相关，并将信号的全部信息集中到一部分具有大幅值的小波系数中。这些大的小波系数含有的能量远比小系数含有的能量大，从而在信号的重构中大的系数比小的系数更重要。
(2)加权平均法。权重系数可调适用范围广，可消除部分噪声，源图像信息损失较少，但会造成图像对比度的下降，需要进行图像灰度增强。
(3)消除高频噪声法。高频噪声基本消除，融合图像对比度较高，源图像特征可较好地保留在融合图像中，但在消除高频噪声的同时损失了部分高频信息。
(4)双阈值法。该法适于源图像中一幅图像的灰度分布均衡，高频成分不明显，另一幅图像高频成分丰富，高灰度级成分较多。双阈值可选增加了算法的实用性，但选择阈值时要考虑源图像灰度分布的特点，否则有可能出现边缘跳跃的现象。
16.1.3 小波包图像融合
一种自然的做法是将尺度子空间Vj
 和小波子空间Wj
 用一个新的子空间
 统一起来表征：
则Hilben空间的正交分解
 的分解统一起来：
定义子空间
 的闭包空间，由上式构造的序列{wn
 (t)}(n∈Z+
 )称为由基函数w0
 (t)=φ(t)确定的小波包。
小波分解的各层具有不同的分辨率，但组合成滤波器组时每一层所取的子带是固定的。小波包分解的每一层滤波器子带都覆盖信号所占的频率，只是各层的频率分辨率不同。
在滤波器组的构造上，小波包分解提供了层间多种组合选择的可能。小波包分解为多少层，以及在各层选择哪些子带来使用十分灵活。对不同特性的信号可以构成相应的最佳滤波器组，小波包图像融合可保存更多的图像细节信息。
16.1.4 小波框架图像融合
一个基函数φ∈L2
 (R)称为小波框架{φj，k
 }的生成元。若
对于有限大的正常数A和B成立，其中A和B分别称为框架的下、上边界，小波系数cj，k
 构成一框架的条件是：
小波框架一般是冗余的，因为它的元素不是线性无关的。一言以蔽之，由框架产生的小波一般是非正交小波。利用小波框架可以得到离散小波和离散化参数各自应该满足的约束条件，并可获得信号重构的精度分析。
16.1.5 多小波图像融合
多小波是指由两个或两个以上的函数作为尺度函数。生成的小波可以同时拥有许多良好的性质，如对称性、短支撑性、正交性和高阶消失矩等。
多小波的基本思想是将单小波中由单个尺度函数生成的多分辨分析空间扩展为由多个尺度函数生成，以此获得更大的自由度。多小波基是由多个小波母函数经过伸缩平移生成的。在实际应用中，可以把光滑性、紧支性、对称性等完美地结合在一起。
图像处理中正交性能保持能量，对称性适合于人眼的视觉系统，使信号边界易于处理。紧支撑的小波对应的滤波器是有限脉冲响应滤波器，使得快速小波变换的和是有限的和。
1994年，Good man 等人基于多分辨分析建立了多小波的理论框架并给出了多小波的例子，同年Geronimo、Har din和Massopust 应用分形插值的方法成功地构造出具有短支撑、正交的、对称的和二阶消失矩的两个尺度函数。
1996年，Geronimo、Hardin和Massopust再次应用分形插值的方法构造了著名的GHM小波。Strela在1996年的博士论文中对多小波在时域和频域上的性质做了比较深刻的分析，提出了两尺度相似变换的概念，并利用这一变换改进和构造多小波。2001年，有学者指出多小波用于图像融合是一个有发展前途的方向。
小波变换法与金字塔图像融合法相比有下列优点。
(1)小波变换的大小与图像相同，而金字塔的大小是图像的4/3，即加大了数据量，小波变换更为紧凑。
(2)小波表达式提供了方向信息，而金字塔表达式未将空间方向选择性引入分解过程。
(3)由于可以选择正交小波核，因此不同分辨率包含的信息是唯一的，而金字塔分解在两个不同的尺度之间含有冗余。另外，金字塔不同级的数据相关。很难知道两级之间的相似性是由于冗余还是图像本身的性质引起的。
(4)金字塔的重构过程可能具有不稳定性，特别是两幅图像存在明显差异区域时，而小波变换图像融合没有类似的问题。
16.1.6 小波分析用于图像融合的实例
图像融合是将同一对象的两个或更多的图像合成在一幅图像中，以便它比原来的任何一幅图像更容易为人们所理解。这一技术可应用于多频谱图像理解以及医学图像处理等领域。在这些场合，同一物体部件的图像往往是采用不同的成像机理得到的。
下面用二维小波分析将两幅图像融合在一起。
【例16.1】利用小波变换法将两幅图像融合。MATLAB代码设置如下：
load laure;
X1=X;
%画出原始图像
image(X1);
title('laure');
axis square
结果如图16-1所示。
load wbarb;
X2=X;
map2=map;
for i=1:256
for j=1:256
if (X2(i，j)>100)
X2(i，j)=3.2*X2(i，j);
else
X2(i，j)=0.5*X2(i，j);
end
end
end
figure;
image(X2);
colormap(map2);
title('wbarb');
axis square
结果如图16-2所示。
图16-1 laure图
图16-2 wbarb图
%用小波函数db4对X1进行2层小波分解
[c1，s1]=wavedec2(X1，2，'db4');
%对分解系数进行处理以突出轮廓部分，弱化细节部分
sizec1=size(c1);
for i=1:sizec1(2)
c1(i)=3.2*c1(i);
end
%用小波函数db4对X2进行2层小波分解
[c2，s2]=wavedec2(X2，2，'db4');
%下面进行小波变换域的图像融合
c=c1+c2(1:size(c1));
%减小图像亮度
c=0.5*c;
%对融合的系数进行重构
xx=waverec2(c，s1，'db4');
%画出融合后的图像
figure;
image(xx);
title('融合图像');
axis square
结果如图16-3所示。
【例16.2】利用小波包法将两幅图像融合。MATLAB代码设置如下：
clear
load woman2;
X1=X;
map1=map;
%画出原始图像
image(X1);
colormap(map1);
title('woman2');
axis square
结果如图16-4所示。
图16-3 融合图
图16-4 图1
load wbarb;
X2=X;
map2=map;
figure;
image(X2);
colormap(map2);
title('wbarb');
axis square
结果如图16-5所示。
%用小波函数db4对X1进行2层小波分解
t1=wpdec2(X1，2，'db4');
%用小波函数db4对X2进行2层小波分解
t2=wpdec2(X2，2，'db4');
%下面进行小波变换域的图像融合
[Lo_D1，Hi_D1，Lo_R1，Hi_R1]=get(t1，'Lo_D'，'Hi_D'，'Lo_R'，'Hi_R');
[Lo_D2，Hi_D2，Lo_R2，Hi_R2]=get(t2，'Lo_D'，'Hi_D'，'Lo_R'，'Hi_R');
Lo_D=(Lo_D1+Lo_D2)*0.1;
Hi_D=(Hi_D1+Hi_D2)*0.5;
Lo_R=(Lo_R1+Lo_R2)*0.5;
Hi_R=(Hi_R1+Hi_R2)*0.5;
t2=set(t2，'Lo_D'，Lo_D，'Hi_D'，Hi_D，'Lo_R'，Lo_R，'Hi_R'，Hi_R);
%对融合的系数进行重构
xx=wprec2(t1);
%画出融合后的图像
figure;
image(xx);
colormap(map);
title('融合图像');
axis square
结果如图16-6所示。
图16-5 图2
图16-6 融合图
16.2 图像分解
回顾从一维离散小波变换到二维的扩展，二维静态小波变换采用相似的方式，对行和列分别采用高通和低通滤波器。这样分解的结果仍然是四组图像、近似图像、水平细节图像、竖直细节图像和对角图像，与离散小波变换不同的只是静态小波分解得到的四幅图像与原图像尺寸一致，原理与一维情况相同。
二维离散小波变换同样只提供了一个函数swt2，因为它不对分解系数进行下采样，所以单层分解和多层分解的结果是一样的，不需要另外提供多层分解的功能。
下面举一个用命令行使用swt命令的例子，大家可以对比它和dwt处理结果的区别。
【例16.3】利用小波变换进行图像分解。MATLAB代码设置如下：
load woman
[swa，swh，swv，swd]=swt2(X，3，'db4');
%使用db4小波对图像进行三层静态小波分解
whos
colormap(map);
kp=0;
for i=1:3
subplot(3，4，kp+1)，image(wcodemat(swa(:，:，i)，192));
title(['Approx，cfs，level'，num2str(i)])
%显示第i层近似系数图像，以192字节为单位编码
subplot(3，4，kp+2)，image(wcodemat(swh(:，:，i)，192));
title(['Horiz.Det.cfs level'，num2str(i)])
subplot(3，4，kp+3)，image(wcodemat(swv(:，:，i)，192));
title(['Vert.Det.cfs level'，num2str(i)])
subplot(3，4，kp+4)，image(wcodemat(swd(:，:，i)，192));
title(['Diag.Det.cfs level'，num2str(i)])
kp=kp+4;
end
结果如图16-7所示。
图16-7 分解图
同idwt2类似，MATLAB对二维静态小波重建提供了iswt2命令，同一维情况类似，对经过重建滤波后的信号不做上采样(因为近似系数和细节系数大小都与原信号一致)。
同一维的静态小波重建一样，将用例子说明如何将 iswt2 单纯用作滤波器来实现各层系数的重建，与一维的情况不同的只是为了重建第j层近似系数，需要4次用到iswt2作为重建滤波器对第j+1 层的系数进行滤波，在对某一个近似系数滤波的过程中同样需要把其他的3个系数指定为0。
为了便于比较，下例接上面的二维静态分解的例子，直接利用对图像的分解结果，从中重构各级系数。
【例16.4】重构系数。MATLAB代码设置如下：
load woman
[swa，swh，swv，swd]=swt2(X，3，'db4');
%使用db4小波对noiswom图像进行三层小波分解
mzero=zeros(size(swd));
A=mzero;
A(:，:，3)=iswt2(swa，mzero，mzero，mzero，'db4');
%使用iswt2的滤波器功能，重建第3层的近似系数，为了避免iswt的合成运算，注意在重建过程中应保证其他各项系数为零
H=mzero;V=mzero;D=mzero;
for i=1;3
swcfs=mzero;swcfs(:，:，i)=swh(:，:，i);
H(:，:，i)=iswt2(mzero，swcfs，mzero，mzero，'db4');
swcfs=mzero;swcfs(:，:，i)=swv(:，:，i);
V(:，:，i)=iswt2(mzero，mzero，swcfs，mzero，'db4');
swcfs=mzero;swcfs(:，:，i)=swh(:，:，i);
H(:，:，i)=iswt2(mzero，mzero，mzero，swcfs，'db4');
end
%分别重建1～3级的各个细节系数，同样在重建某一系数的时候，要令其他系数为0
A(:，:，2)=A(:，:，3)+H(:，:，3)+V(:，:，3)+D(:，:，3);
A(:，:，1)=A(:，:，2)+H(:，:，2)+V(:，:，2)+D(:，:，2);
%使用递推的方法建立第1层和第2层近似系数
colormap(map)
kp=0;
for i=1:3
subplot(3，4，kp+1)，image(wcodemat(A(:，:，i)，192));
title(['第'，num2str(i)，'层近似系数图像'])
subplot(3，4，kp+2)，image(wcodemat(H(:，:，i)，192));
title(['第'，num2str(i)，'层水平细节系数图像'])
subplot(3，4，kp+3)，image(wcodemat(V(:，:，i)，192));
title(['第'，num2str(i)，'层竖直细节系数图像'])
subplot(3，4，kp+4)，image(wcodemat(D(:，:，i)，192));
title(['第'，num2str(i)，'层对角细节系数图像'])
kp=kp+4;
end
%画出通过手工方法重建的各级小波系数图像
err=norm(A(:，:，2)-swa(:，:，2))
%求出用这种算法重建的第2层近似系数和分解系数之间的误差
结果如图16-8所示。
图16-8 结果图
16.3 图像特征提取
图像特征在图像识别和计算机分析上十分有用，在前面章节中已经介绍了信号特征提取。本节限于篇章，只介绍相关图像特征提取的算例。
【例16.5】利用小波变换提取其中一个子带的直方图。MATLAB代码设置如下：
close all
clear all
clc;
Original_image=imread('moon.tif');
Original_image=double(Original_image);
image((Original_image));
colormap('default');
title('原始图');
结果如图16-9所示。
%小波变换进行系数分解
[ll1，h1，v1，vd]=dwt2(Original_image，'db4');
HL1=v1;
HL1_1D=sort(HL1(:)).';
bin_size=2;
range_value=81;
step=[-range_value:bin_size:range_value];
figure
hist(HL1_1D);
title('小波垂直分解直方图');
结果如图16-10所示。
图16-9 原图
图16-10 直方图
Extract_histogram1=histc(HL1_1D，step);
% 提取的 histogram 步长为1
c_elements=cumsum(Extract_histogram1);
figure;
bar(c_elements，'BarWidth'，1)
title('提取出的子直方图')
结果如图16-11所示。
图16-11 提取子直方图
【例16.6】利用小波变换，对给定图像进行旋转变换后，提取其中一个子带图。MATLAB代码设置如下：
clear
Original_image=imread('saturn.png');
imshow(Original_image);
title('原图')
结果如图16-12所示。
Rescale_image=imrotate(Original_image，10);
figure;
imshow(Rescale_image);
title('旋转图')
结果如图16-13所示。
图16-12 原图
图16-13 旋转图
Rescale_image=double(Rescale_image);
[ll4，h4，v4，d4]=dwt2(Rescale_image，'db4');
HL4=v4;
HL4_1D=sort(HL4(:)).';　%将HL子带变为一维向量
bin_size=2;
range_value=81;
step=[-range_value:bin_size:range_value];
Extract_histogram4=histc(HL4_1D，step); %提取的 histogram 步长为1
[m4，n4]=size(Extract_histogram4);
n1=n4;
% hist(sort(HL(:))，n);
relation4=zeros(1，n1);
for i=2:3:(n1/2)
relation4(1，i)=2*Extract_histogram4(1，i)/(Extract_histogram4(1，i+1)+Extract_histogram4(1，i-1));
end
for i=((n1/2)+4):3:(n1-2)
relation4(1，i)=2*Extract_histogram4(1，i)/(Extract_histogram4(1，i+1)+Extract_histogram4(1，i-1));
end
%%% Extract watermark W
W_number4=24;　%水印的个数 W_number4
H4=zeros(1，3*W_number4);
H4(1:(3*W_number4/2))=Extract_histogram4(1:(3*W_number4/2)); % H的第一个直方图的起点是-range_value
H4((3*W_number4/2+1):3*W_number4)=Extract_histogram4((n4/2+3):(n4/2+2+3*W_number4/2));% H的第(2*W_number+1)个直方图的起点是-range_value+bin_size*n/2
a=zeros(1，W_number4);
b=zeros(1，W_number4);
c=zeros(1，W_number4);
H_matrix4=reshape(H4，3，W_number4);
a=H_matrix4(1，:);
b=H_matrix4(2，:);
c=H_matrix4(3，:);　% a(i)，b(i) and c(i)是嵌入第i个水印的3个bin的大小
extract_W4=zeros(1，W_number4);
for i=1:W_number4
if (2*b(1，i) >=a(1，i)+c(1，i))
extract_W4(1，i)=1;
else if (2*b(1，i) < (a(1，i)+c(1，i)))
extract_W4(1，i)=-1;
end
end
end
extract_watermark4=extract_W4;
X4=2:1:(n1-2);
Y4=relation4(1，X4);
figure;
plot(X4，Y4，'gd');
title('旋转后的二进提取');
结果如图16-14所示。
图16-14 结果图
【例16.7】提取小波分解概貌系数。MATLAB代码设置如下：
clear
Original_image=imread('saturn.png');
subplot(121);
imshow(Original_image);
title('原图')
Input=Original_image;
A=double(Input);
Xrgb=0.2990*A(:，:，1)+0.5870*A(:，:，2)+0.1140*A(:，:，3);
NbColor=255;
X1=wcodemat(Xrgb，NbColor);
%伪彩矩阵压缩
[c2，l2]=wavedec2(X1，4，'sym4');
%二维小波分解
ch11=detcoef2('h'，c2，l2，1);
%提取小波分解细节系数
ch12=detcoef2('h'，c2，l2，2);
ch13=detcoef2('h'，c2，l2，3);
ch14=detcoef2('h'，c2，l2，4);
cv11=detcoef2('v'，c2，l2，1);
cv12=detcoef2('v'，c2，l2，2);
cv13=detcoef2('v'，c2，l2，3);
cv14=detcoef2('v'，c2，l2，4);
cd11=detcoef2('d'，c2，l2，1);
cd12=detcoef2('d'，c2，l2，2);
cd13=detcoef2('d'，c2，l2，3);
cd14=detcoef2('d'，c2，l2，4);
ca14=appcoef2(c2，l2，'sym4'，1);
%提取小波分解概貌系数
subplot(122);
image(ca14);
colormap;
axis off
结果如图16-15所示。
图16-15 概貌图
16.4 本章小结
本章主要介绍了图像处理中小波方法的应用。小波分析的应用是与小波分析的理论研究紧密地结合在一起的。现在，它已经在科技信息领域取得了令人瞩目的成就。电子信息技术是六大高新技术中的一个重要领域，图像和信号处理又是电子信息技术领域的重要方面。现今，信号处理已经成为当代科学技术工作的重要组成部分。
现在，对性质随时间稳定不变的信号处理的理想工具仍然是傅立叶分析。但在实际应用中，绝大多数信号是非稳定的，小波分析正是适用于非稳定信号的处理工具。图像处理是针对性很强的技术，根据不同应用、不同要求需要采用不同的处理方法。
图像处理是“信息处理”的一个方面，这一观点现在已经为人所熟知。它可以进一步细分为多个研究方向：图片处理、图像处理、模式识别、景物分析、图像理解、光学处理等。小波分析在图像处理方面的应用主要是进行图像压缩、图像去噪、图像增强(包括图像钝化和图像锐化)、图像融合、图像分解及图像特征提取等。
本书还详细介绍了用MATLAB来实现这些图像处理的技术。MATLAB的语法规则比FORTRAN语言和C语言等高级语言更简单，更重要的是，它具有贴近人的思维方式的编程特点。使用MATLAB编写程序有如在便笺上列公式和求解，所以了解这些技术对读者选用工具箱中合适的函数来解决自己的问题很有帮助。