图像融合

图像融合是将两幅或多幅图像融合在一起，以获取对同一场景的更为精确、更为全面、更为可靠的图像描述。

融合算法应该充分利用各原图像的互补信息，使融合后的图像更适合人的视觉感受，适合进一步分析的需要；并且应该统一编码，压缩数据量，以便于传输。

图像融合可分为三个层次:1. 像素级融合2. 特征级融合3. 决策级融合其中像素级融合是最低层次的融合，也是后两级的基础。

它是将各原图像中对应的像素进行融合处理，保留了尽可能多的图像信息，精度比较高，因而倍受人们的重视。

像素级的图像融合方法大致可分为三大类:1. 简单的图像融合方法2. 基于塔形分解（如Laplace塔形分解、比率塔等）的图像融合方法3. 基于小波变换的图像融合方法小波变换是图像的多尺度、多分辨率分解，它可以聚焦到图像的任意细节，被称为数学上的显微镜。

近年来，随着小波理论及其应用的发展，已将小波多分辨率分解用于像素级图像融合。

小波变换的固有特性使其在图像处理中有如下优点：1. 完善的重构能力，保证信号在分解过程中没有信息损失和冗余信息；2. 把图像分解成平均图像和细节图像的组合，分别代表了图像的不同结构，因此容易提取原始图像的结构信息和细节信息；3. 具有快速算法，它在小波变换中的作用相当于FFT算法在傅立叶变换中的作用，为小波变换应用提供了必要的手段；4. 二维小波分析提供了与人类视觉系统方向相吻合的选择性图像。

——像素级图像融合的主要步骤以两幅图像的融合为例。

设A，B为两幅原始图像，F为融合后的图像。

若对二维图像进行N层的小波分解，最终将有(3N+1)个不同频带，其中包含3N 个高频子图像和1个低频子图像。

其融合处理的基本步骤如下：（1）对每一原图像分别进行小波变换，建立图像的小波塔型分解；（2）对各分解层分别进行融合处理。

各分解层上的不同频率分量可采用不同的融合算子进行融合处理，最终得到融合后的小波金字塔；（3）对融合后所得小波金字塔进行小波重构，所得到的重构图像即为融合图像。

图 1在图像融合过程中，小波基的种类和小波分解的层数对融合效果有很大的影响，对特定的图像来说，哪一种小波基的融合效果最好，分解到哪一层最合适，都是需要考虑的问题。

为此可以通过引入融合效果的评价来构成一个闭环系统。

如图2所示。

对图像而言，小波变换是将图像分解成频域上各个频率段的子图，以代表原图的各个特征分量。

这对后续的融合处理极为重要，使得融合处理可以根据不同的特征分量采用不同的融合方法以达到最佳融合效果。

图像的融合策略（方法）是图像融合的核心，方法与规则的优劣直接影响融合的速度与质量。

在一幅图像的小波分解中，绝对值较大的小波高频系数对应着亮度急剧变化的点，也就是图像中对比度变换较大的边缘特征，如边界、亮线及区域轮廓。

融合的效果就是对同样的目标，融合前在图像A中若比图像B中显著，融合后图像A中的目标就被保留，图像B中的目标就被忽略。

这样，图像A、B中目标的小波变换系数将在不同的分辨率水平上占统治地位，从而在最终的融合图像中，图像A 与图像B中的显著目标都被保留。

目前小波域的融合规则主要分为两种：一、基于单个像素的融合规则主要包括:（1）小波系数的直接替换或追加；（2）最大值选取；（3）加权平均。

二、基于区域特征的融合规则。

主要包括:（1）基于梯度的方法；（2）基于局域方差的方法；（3）基于局域能量的方法。

基于像素的融合规则在融合处理时表现出对边缘的高度敏感性，使得在预处理时要求图像是严格对准的，否则处理结果将不尽人意，这就加大了预处理的难度。

基于区域的融合规则由于考虑了与相邻像素间的相关性，降低了对边缘的敏感性，所以具有更加广泛的适用性。

图1针对不同类型的图像，下面介绍几种常用的融合方法：（1）取系数绝对值较大法适合高频成分较丰富，亮度、对比度较高的原图像，否则在融合图像中只保留一幅图像的特征，其他的特征被覆盖；融合图像中基本保留原图像的特征，图像对比度与原图像基本相同。

小波变换的实际作用是对信号解相关，并将信号的全部信息集中到一部分具有大幅值的小波系数中。

这些大的小波系数含有的能量远比小系数含有的能量大，从而在信号的重构中，大的系数比小的系数更重要。

（2）加权平均法权重系数可调，适用范围广，可消除部分噪声，原图像信息损失较少，但会造成图像对比度的下降，需要进行图像灰度增强。

（3）消除高频噪声法高频噪声基本消除，融合图像对比度较高，原图像特征可较好地保留在融合图像中，但在消除高频噪声的同时，损失了部分高频信息。

（4）双阈值法适于原图像中一幅图像的灰度分布均衡，高频成分较多；双阈值可选，增加了算法的实用性，但选择阈值时要考虑原图像灰度分布的特点，否则有可能出现边缘跳跃的现象。

前几天我们讨论了一种基于低频融合策略的小波图像融合算法（/chenyusiyuan/archive/2007/11/14/1883987.aspx），今天就来说说算法编程实现的几个细节问题。

（1）高频带的融合首先，高频系数融合均采用基于像素点绝对值取大的规则，这里要注意的是，比较的对象是两幅原始图像的像素点的绝对值，而融合后图像的像素点的取值是两幅原始图像像素点绝对值大的那一点的值，也就是说，要取绝对值大的那个点的实际值。

故相应的Matlab 程序为：% I1，I2 为两幅原始图像y1=mywavedec2(I1,dim);y2=mywavedec2(I2,dim);[r,c]=size(y1); % 根据低频融合算法进行图像融合for i=1:r % 首先取两幅源图像相应的小波分解系数绝对值最大者的值作为融合图像的分解系数for j=1:cif ( abs(y1(i,j)) >= abs(y2(i,j)) )y3(i,j)=y1(i,j);elseif ( abs(y1(i,j)) < abs(y2(i,j)) )y3(i,j)=y2(i,j);end% y3(i,j)=max( abs(y1(i,j)) , abs(y2(i,j)) ); % 这行代码是错误的，取了绝对值而非实际值endend（2）低频融合的对象低频融合的对象是N级分解后图像的LL-N子图像，相应的程序代码为：LLa=y1(1:r/(2^dim),1:c/(2^dim)); % 调用lowfrefus函数对低频部分的小波分解系数进行融合LLb=y2(1:r/(2^dim),1:c/(2^dim));y3(1:r/(2^dim),1:c/(2^dim))=lowfrefus(LLa,LLb);（3）低频融合策略的程序代码这一部分程序代码的实现难度不大，根据算法原理中列出的公式可以很快写出相应的代码。

不过有些细节问题还是要花心思处理，比如区域大小的确定、区域边界与图像边界的关系、区域中心与区域中各点的权值确定、区域中心在原始图像中的具体位置等等。

function y = lowfrefus(A,B);% 函数 Y=LOWFREFUS(A,B) 对输入的两个小波分解系数矩阵，根据低频融合算法，得出融合图像的低频小波分解系数[row,col]=size(A); % 求出分解系数矩阵的行列数alpha=0.5; % alpha是方差匹配度比较的阈值for i=1:row % 根据低频融合算法，先求出矩阵A,B中以点P为中心的区域方差和方差匹配度for j=1:col % 再根据方差匹配度与阈值的比较确定融合图像的小波分解系数[m2p(i,j),Ga(i,j),Gb(i,j)] = area_var_match(A,B,[i,j]);Wmin=0.5-0.5*((1-m2p(i,j))/(1-alpha));Wmax=1-Wmin;if m2p(i,j)<alpha % m2p表示方差匹配度if Ga(i,j)>=Gb(i,j) % 若匹配度小于阈值，则取区域方差大的相应点的分解系数作为融合图像的分解系数y(i,j)=A(i,j);elsey(i,j)=B(i,j);endelse % 若匹配度大于阈值，则采取加权平均方法得出相应的分解系数if Ga(i,j)>=Gb(i,j)y(i,j)=Wmax*A(i,j)+Wmin*B(i,j);elsey(i,j)=Wmin*A(i,j)+Wmax*B(i,j);endendendendfunction w = weivec(x,p);% 函数 W=WEIVEC(X,P) 对输入的r*c矩阵，计算出以点p为中心时矩阵各点的对应权值% 距离点p越近，权值就越大。

权值是通过行和列的高斯分布加权相加得到的。

[r,c]=size(x);p1=p(1); p2=p(2);sig=1;for i=1:rfor j=1:cw(i,j)=0.5*(gaussmf(i,[sig p1])+gaussmf(j,[sig p2]));endendfunction [m2p,Ga,Gb] = area_var_match(A,B,p);% 函数 area_var_match 计算两个输入矩阵以点p为中心的区域方差以及区域方差匹配度level=1; % 设置区域的大小[subA,mpa,npa]=submat(A,p,level); % submat 函数取输入矩阵中以点P为中心、阶数为（2*level+1）的方阵作为子矩阵[subB,mpb,npb]=submat(B,p,level);[r,c]=size(subA);w=weivec(subA,[mpanpa]); % 获取子矩阵的权值分布averA=sum(sum(subA))/(r*c); % 计算子矩阵的平均值averB=sum(sum(subB))/(r*c);Ga=sum(sum(w.*(subA-averA).^2)); % 计算子矩阵的区域方差Gb=sum(sum(w.*(subB-averB).^2));if (Ga==0)&(Gb==0) % 计算两个子矩阵的区域方差匹配度m2p=0;elsem2p=2*sum(sum(w.*abs(subA-averA).*abs(subB-averB)))/(Ga+Gb);endfunction [smat,mp,np] = submat(x,p,level);% 函数 submat 取输入矩阵中以点P为中心、阶数为（2*level+1）的方阵作为输出的子矩阵[row,col]=size(x);m=p(1); n=p(2);if (m>row)||(n>col)error('Point p is out of matrix X !');return;endif ((2*level+1)>row)||((2*level+1)>col)error('Too large sample area level !');return;end% 设置子矩阵的边界值up=m-level; down=m+level;left=n-level; right=n+level;% 若子矩阵的某一边界值超出输入矩阵的相应边界，就进行边界处理，% 即超出边界后往相反方向平移，使其恰好与边界重合if left<1right=right+1-left;left=1;endif right>colleft=left+col-right;right=col;endif up<1down=down+1-up;up=1;endif down>rowup=up+row-down;down=row;end% 获取作为输出的子矩阵，并计算点p在输出的子矩阵中的位置smat = x(up:down,left:right);mp=m-up+1;np=n-left+1;下面我们给出相应的图片示例，图1 是一幅远景清晰近景模糊的图片，图2 是一幅远景模糊近景清晰的图片，图3是融合效果图。