图像融合基于matlab拉普拉斯金字塔图像融合含Matlab 746期

Posted 2022-05-25 海神之光

一、获取代码方式

获取代码方式1：
通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。

获取代码方式2：
完整代码已上传我的资源：【图像融合】基于matlab拉普拉斯金字塔图像融合【含Matlab 746期】

备注：
订阅紫极神光博客付费专栏，可免费获得1份代码（有效期为订阅日起，三天内有效）；

二、拉普拉斯金字塔简介

1 拉普拉斯金字塔
在高斯金字塔的运算过程中，图像经过卷积和下采样操作会丢失部分高频细节信息。为描述这些高频信息，人们定义了拉普拉斯金字塔(Laplacian Pyramid， LP)。用高斯金字塔的每一层图像减去其上一层图像上采样并高斯卷积之后的预测图像，得到一系列的差值图像即为 LP 分解图像。
将Gl内插方法得到放大图像Gl，使Gl的尺寸与Gl-1的尺寸相同，即放大算子Expand。
该式子实现两个步骤：在偶数行和列插入0，然后使用下采样中的高斯核进行滤波处理，得到和l-1层一样大小的图像。
N为拉普拉斯金字塔顶层的层号LPl是拉普拉斯金字塔分解的第L层图像。由LP0，LP1、LP2…LPN构成的金字塔即为拉普拉斯金字塔。它的每一层L0图像是高斯金字塔本层G0图像与其高一层图像G1经内插放大后图像G1的差，此过程相当于带通滤波，因此拉普拉斯金字塔又称为带通金字塔分解。
内插方法：opencv中有实现的函数pyrup。可以得到*G1。然后在两个函数作差，相减就可以得到拉普拉斯金字塔。
求得每个图像的拉普拉斯金字塔后需要对相应层次的图像进行融合，具体的融合规则有，取大、取小，等等。

首先对原始图像Gi进行向下取样，然后向上采用，最后与原始图像相减，得到拉普拉斯金字塔图像。
下面这张图也是比较常见的：

使用原始图像 套入公式得到拉普拉斯金字塔第0层。
使用原始图像向下采样Od 代入公式 得到 拉普拉色金字塔第1层。
这张图在看一些资料的时候也经常看到：

2 融合应用
　　图像拉普拉斯金字塔分解的目的是将源图像分别分解到不同的空间频带上，融合过程是在各空间频率层上分别进行的，这样就可以针对不同分解层的不同频带上的特征与细节，采用不同的融合算子以达到突出特定频带上特征与细节的目的。即有可能将来自不同图像的特征与细节融合在一起。
（1）顶层处理
设LAl和LBl分别为源图像A,B经过拉普拉斯金字塔分解后得到的第l层图像，融合后的结果为LFl。当l=N时，LAN和LBN分别为源图像A，B经过拉普拉斯金字塔分解后得到的顶层图像。对于顶层图像的融合，首先计算以其各个像素为中心的区域大小为M*N(M、N取奇数且M >= 3、N >= 3)的区域平均梯度：
其中，Ix与Iy分别为像素f(x,y)在x与y方向上的一阶差分，定义如下：
因此对于顶层图像中的每一个像素LAN(i, j)和LBN(i, j)都可以得到与之相对应的区域平均梯度GA(i, j)和GB(i, j)。由于平均梯度反映了图像中的微小细节反差和纹理变化特征，同时也反映出图像的清晰度。一般来说平均梯度越大，图像层次也丰富，则图像越清晰。因此顶层图像的融合结果为：
（2）各层次处理
当0<l<N时，则对于经过拉普拉斯金字塔分解的第l层图像，首先计算其区域能量：
则其他层次图像的融合结果为：
在得到金字塔各个层次的融合图像LF1、LF2、LFN后。通过前面的重构，便可得到最终的融合图像。
第二种融合规则：
采用最高层系数取平均，其余各层系数绝对值取大的融合策略进行融合。融合后图像的系数（灰度值）越接近较清晰图像的灰度值就说明融合效果好。

三、部分源代码

clear all;
close all;
clc;

im1 = (imread('红外.jpg'));
im2 = imread('可见光.jpg');

figure(1);
imshow(im1);
title('红外图像');
figure(2);
imshow(im2);
title('可见光图像');

if length(size(im1)) ==3%判断输入图像是否是真彩色图像
    im1 = rgb2gray(im1);
end
if length(size(im2)) ==3
    im2 = rgb2gray(im2);
end
im1 = double(im1);
im2 = double(im2);
%%
%拉普拉斯滤波器
w = [1 4 6 4 1; 4 16 24 16 4; 6 24 36 24 6; 4 16 24 16 4; 1 4 6 4 1]/256;
G = cell(1,5);
H = cell(1,5);
I = cell(1,5);
G1 = im1;%第一层为原图像
H1 = im2;

for k=2:5
    Gk = imfilter(Gk-1,w,'replicate');%滤波
    Hk = imfilter(Hk-1,w,'replicate');%滤波
    Gk = Gk(1:2:size(Gk,1),1:2:size(Gk,2)); %下采样
    Hk = Hk(1:2:size(Hk,1),1:2:size(Hk,2)); %下采样
end

for k=5:-1:2
    Gk-1 = Gk-1(1:2*size(Gk,1),1:2*size(Gk,2));%调整图像尺寸
    Hk-1 = Hk-1(1:2*size(Hk,1),1:2*size(Hk,2));%调整图像尺寸
end

for k=1:4
    
    Gk = Gk - us(Gk+1,w);%获得拉普拉斯金字塔
    Hk = Hk - us(Hk+1,w);
    a = abs(Gk) > abs(Hk);%两幅图像的1-4级拉普拉斯金字塔取绝对值大的
    Ik = a.*Gk + ~a.*Hk;
end
I5 = H5;%定义融合后顶层图像尺寸，并为图像边缘赋值，图像中间的值会被比较平均梯度后所取的值替换
gg = g(G5);
gh = g(H5);

for i = 2:size(G5,1)-2
    for j = 2:size(G5,2)-2
        if gg(i,j) > gh(i,j)%比较某点邻域平均梯度大小，取大的
            I5(i,j) = G5(i,j);
        else
            I5(i,j) = H5(i,j);
        end
    end
end

for k=4:-1:1
    Ik = Ik+us(Ik+1,w);%重构图像
end

function f = g(b)%区域平均梯度函数（3*3区域）
[Y,X] = size(b);
g = zeros(Y,X);%保存梯度
f = g;%保存平均梯度

for i = 1:Y-1%求梯度，最后一行、列求不到，后面会
    for j = 1:X-1
        g(i,j)  = sqrt((((b(i+1,j)-b(i,j))^2+(b(i,j+1)-b(i,j))^2))/2);%某点梯度
    end
end
function f = us(im,w)
w=w*4;
[i,j] = size(im);
[m,n] = size(w);

w11 = w(1:2:m,1:2:n);%插入滤波器
w12 = w(1:2:m,2:2:n);
w21 = w(2:2:m,1:2:n);
w22 = w(2:2:m,2:2:n);

四、运行结果

五、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 蔡利梅.MATLAB图像处理——理论、算法与实例分析[M].清华大学出版社，2020.
[2]杨丹,赵海滨,龙哲.MATLAB图像处理实例详解[M].清华大学出版社，2013.
[3]周品.MATLAB图像处理与图形用户界面设计[M].清华大学出版社，2013.
[4]刘成龙.精通MATLAB图像处理[M].清华大学出版社，2015.