使用小波变换进行灰度图像的融合

参考：
https://blog.youkuaiyun.com/liu_guanzhang/article/details/42120861
https://wenku.baidu.com/view/a4cfe4c36137ee06eff918c0.html
本文旨在利用小波变换实现一个很简单的融合算法，对小波变换的应用有个粗略的认识。当然本文的算法其实是很简陋的，仅是作为一个例子，不具有广泛应用的实际。
在下面的程序中，我们基于这样一个简单的融合规则，即对2幅图像的低频部分采取均值，而对2幅图像的高频部分，我们留下对应位置最大的那个频率信息。
算法的流程图如下图所示：

融合程序

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % 使用小波变换对2幅图像进行融合，仅适用灰度图
% % % % % %参数说明：
%    M1     - 输入图像 A
%    M2     - 输入图像 B 
%    wtype  - 指定使用的小波基
%    Y      - 融合结果图像
% % 
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

clc;
clear;
close all;                
% % ===============================================
%   读入要进行融合的图像
% % ===============================================
imA = imread('pan.tif');
% 这里转化成double类型，否则使用小波变换输出的会有大量大于1的存在，会导致图像显示有问题
M1 = double(imA) / 256;   % 【问题1】
imB = imread('LR.tif');
imB = imresize(imB,4);
imB = rgb2gray(imB);
M2 = double(imB) / 256;


zt = 4;           % 设置小波分解的层数
wtype = 'haar';  % 设置小波基

% % ===============================================================
% 这里的融合规则是：在高频系数中，取2幅图像绝对值大的小波系数，
% 在低频系数中，取2幅图像的平均值。
% % ================================================================

% % ==============================================
% 对2幅图像分别进行小波分解
% % ==============================================
[c0,s0] = wavedec2(M1, zt, wtype);  
[c1,s1] = wavedec2(M2, zt, wtype); 


KK = size(c1);
Coef_Fusion = zeros(1,KK(2));  % 生成同系数矩阵c一样大的矩阵作为结果存储单元

Coef_Fusion(1:s1(1,1)*s1(1,2)) = (c0(1:s1(1,1)*s1(1,2))+c1(1:s1(1,1)*s1(1,2)))/2;


% % =============================================
% % 处理高频系数
% % =============================================
MM1 = c0(s1(1,1)*s1(1,2)+1:KK(2));
MM2 = c1(s1(1,1)*s1(1,2)+1:KK(2));
mm = (abs(MM1)) > (abs(MM2));
Y  = (mm.*MM1) + ((~mm).*MM2);
Coef_Fusion(s1(1,1)*s1(1,2)+1:KK(2)) = Y;


% % =============================================
% 重构
% % 这里我们已经生成新的系数矩阵了，并利用新的系数矩
% % 阵，进行小波逆变换
% % =============================================
Y = waverec2(Coef_Fusion,s0,wtype);

% % =============================================
% 显示图像
% % =============================================
subplot(2,2,1);imshow(M1);title('pan');
subplot(2,2,2);imshow(M2);title('LR');
subplot(223);imshow(Y,[]);title('融合图像');

结果

为什么要有这句代码M1 = double(imA) / 256

需要注意的是，为什么要将图像放缩到0-1？我们看下面一段程序就知道原因了。即不加这句代码，会造成图像显示为“全白”的现象。

clc;
clear;
close all;

x=imread('lena_color_512.tif');  %读取原图像
x = double(x)/256; % 注释掉这段代码，再看现象
figure(1);imshow(x);title('原始图像');
x=double(x);
[c,s ]=wavedec2(x,1,'db8'); %作小波变换,形成小波系数矩阵
r = waverec2(c,s,'db8');%还原图像文件
figure;imshow(r);title('重构之后的图像');