基于Retinex增强和拉普拉斯金字塔分解的图像融合 SSR-Laplacian Image Fusion [转载]

AN IMPROVED FUSION METHOD FOR INFRARED AND LOW-LIGHT LEVEL VISIBLE IMAGE - 2017

​   该论文是一篇SCI二区的论文，发表在TIP上，本文提出了一种融合低照度可见光图像和红外图像的方法。关于本文的创新点，并不是很强，主要是增加了图像的预处理步骤，后续过程多是一些传统处理方法的组合和延伸。由于创新点并不是很强，因此并没有细读，在这里做一个简单的记录。

算法框架

  由图可概括该算法的流程如下：

• 可见光图像的预处理：SSR算法，一种基于Retinex的人眼生理特征所发展起来的算法。这里用到了最基本的Single Scale Retinex算法，用于提升可见光图像的对比度。由于低照度图像本身的信息不多，因此增强可见光图像对于融合的质量理论上会有较大的提升。
• 拉普拉斯金字塔分解：拉普拉斯金字塔分解的图像融合本身就是一个很经典的算法，在这里并没有什么更多的创新之处。
• 融合策略：拉普拉斯金字塔的各种融合算法差异就在于融合策略的选取，本文所用的是一种局部能量特征法，该法的起源也可以说是很早了，所以也谈不上创新。

算法实现

  该算法本身的借鉴意义有限，但是几个子算法都有值得推敲的地方。在此，以Matlab代码的形式，对照上述的算法框架，做一个简单的实现。
  SSR算法，我们看到的图像一般都是光线和物体共同作用的结果，因此光线的好坏，会很大程度上影响图像的质量。SSR算法则是从人眼特性的角度，剥离光照的影响，还原物体的内在属性。关于此算法，有很多变体算法和延伸内容，如有兴趣可以自行学习，在此给出百度的SSR。本文，只给出单尺度SSR的示例。

function out = SSR(im, sigma)
%Single Scale Retinex
%分通道处理

if ~exist('sigma', 'var')
    sigma = 80;
end


hsize = size(im(:,:,1));
out = zeros([hsize, 3]);

k = fspecial('gaussian', hsize, sigma);
G = imfilter(im, k, 'replicate');
for i=1:size(im, 3)
    L = G(:,:,i);
    r = log(im(:,:,i)+1)-log(L+1);

    R = exp(r)-1;
    Max = max(R(:));
    Min = min(R(:));
    out(:,:,i) = (R-Min)/(Max-Min);
end
end

  效果如图，下图是采用不同的尺度得到的结果。如何选取sigma取决于你需要更好的还原颜色，还是更好的提取特征。其实我不知道这一步在实际做低照度夜视图像时有什么用。

  拉普拉斯金字塔分解，相关资料较多，不多做解释。此处给出彩色图像和红外图像的拉普拉斯金字塔融合效果。

function pyr = LPD(I, nlev, sigma)
%Laplacian Pyramid Decomposition
%拉普拉斯金字塔分解，在灰度图进行，多通道图片建议分通道处理。
%		nlev
%				可选，用于调节金字塔层数
%		sigma
%				可选，用于调节高斯模糊程度
%

if size(I,3)==3
    I=rgb2gray(I);
end

if ~exist('nlev', 'var')
    nlev = 4;
end
 
if ~exist('sigma', 'var')
    sigma=0.4;
end

f = fspecial('gaussian', [5, 5], sigma);
 
%构建拉普拉斯金字塔
pyr = cell(nlev,1);
J = I;
for l = 1:nlev-1   
    J_gauss = imfilter(J,f,'replicate');
    J_gauss_down = J_gauss(1:2:size(J_gauss,1)-1,1:2:size(J_gauss,2)-1); %downsample 
    J_gauss_high = imresize(J_gauss_down,[size(J_gauss,1) size(J_gauss,2)],'bilinear');
    pyr{l} = J-J_gauss_high;
    J=J_gauss_down;
end
pyr{nlev}=J_gauss_down; %拉普拉斯金字塔顶层即为高斯金字塔顶层

end

  融合策略，基于局部能量特征（存疑）。由于文中并未解释此法，在网上搜了一下，发现早在0几年的时候就用于小波变换的图像融合中，但是找到的文献公式，感觉并不正确，以下是个人理解。

  现在开始对于两个图像金字塔的各层进行融合。对于顶层图像，采用取两个图中较大值的方法融合。对于下面的几层，采取以下的融合规则：

  如图， M M M指两幅图像某个像素点的匹配度， E E E指对于某个像素点的能量。 E = E= E=某个点的局部窗口的平方和， M = M= M=两个图像对应点的局部窗口乘积和的平方除以各自的能量，懒的敲公式，并且觉得公式也有问题，具体可以看代码理解。当匹配度 M M M小于阈值 e e e，即匹配度低的时候，说明两个点差异大，有一个像素点价值较低，舍去。具体的取舍方法则是选择 E E E较大的点，即能量高。同理，当匹配度大于阈值，说明两个点接近，都有价值，则进行加权融合，权重计算方法：

​    
 
依旧是能量大的用大权重，能量小的用小权重。代码如下：

function pyr = pyr_fusion(pyr1,pyr2, e)
%pyramid fusion
%局部能量特征的拉普拉斯金字塔融合

if ~exist('e', 'var')
    e = 0.6;
end

n = size(pyr1, 1);
pyr = cell(n, 1);
pyr{n} = max(pyr1{n}, pyr2{n});%顶层取较大

f = fspecial('average', [3 3]);

% 对每一层融合
for i=n-1:-1:1
    F11 = pyr1{i};
    F12 = pyr2{i};
    E1 = imfilter(F11.^2, f, 'replicate');%计算能量
    E2 = imfilter(F12.^2, f, 'replicate');

    M = imfilter(F11.*F12, f, 'replicate').^2./(E1.*E2);%计算匹配度
    F11(E1<E2) = F12(E1<E2);%匹配度低时
    
    W_min = 0.5*(1 - (1-M)./(1-e));%计算权重
    W_max = 1 - W_min;
    F21 = W_min.*pyr1{i}+W_max.*pyr2{i};
    F22 = W_min.*pyr2{i}+W_max.*pyr1{i};
    F21(E1>E2) = F22(E1>E2);%匹配度高时
    
    F11(M>e) = F21(M>e);%按匹配度计算
    pyr{i} = F11;
end

  这一步完成了金字塔的融合，最后从金字塔重建，还原图像。附上主程序代码。

clc
clear
close all
addpath ./func
path = 'im/';

rgb = double(imread([path 'vis.jpg']));
nir = im2double(imread([path 'inf.jpg']));

rgb_s = SSR(rgb);
% rgb_s = rgb/255;

for i=1:3
    rgb_pyr = LPD(rgb_s(:,:,i));%金字塔分解
    nir_pyr = LPD(nir);
    
    pyr = pyr_fusion(rgb_pyr, nir_pyr, 0.6);%金字塔融合
    
%     for j=1:length(pyr)
%         figure;imshow(mat2gray(pyr{j}));
%     end

    out = pyr{4};%金字塔重建
    for l = 3 : -1 : 1
        out = pyr{l} + imresize(out,[size(pyr{l},1) size(pyr{l},2)],'bilinear');
    end
    O(:,:,i) = out;
end

imshow([rgb/255 nir O])

  融合结果如下所示，从左到右依次是可见光图像，红外图，融合图。

图像融合专题
​GitHub仓库