LDR图像融合HDR: MertensTMO方法

参考文献：

"Exposure Fusion" by Tom Mertens, Jan Kautz, Frank Van Reeth in Proceedings of Pacific Graphics 2007

 

（1）高斯金字塔

高斯金字塔是最基本的图像塔。首先将原图像作为最底层图像G0（高斯金字塔的第0层），利用高斯核（5*5）对其进行卷积，然后对卷积后的图像进行下采样（去除偶数行和列）得到上一层图像G1，将此图像作为输入，重复卷积和下采样操作得到更上一层图像，反复迭代多次，形成一个金字塔形的图像数据结构，即高斯金字塔。

高斯金字塔的当前层图像就是对其前一层图像首先进行高斯低通滤波，然后再进行隔行和隔列的降2采样而生成的。前一层图像大小依次为当前层图像大小的4倍。

Opencv中使用pyrdown函数可获得高斯金字塔。

function p = pyrGaussGen(img)
[r, c, ~] = size(img);
levels = floor(log2(min(r, c)));
list = [];
for i=1:levels
    %Detail layer
    ts   = struct('detail', img);
    list = [list, ts];  
    %Next level
    img = pyrGaussGenAux(img);
end
%Base layer
p = struct('list', list, 'base', img);
end

function imgOut = pyrGaussGenAux(img)
%5x5 Gaussian Kernel
kernel = [1, 4, 6, 4, 1];
mtx = kernel' * kernel;
mtx = mtx / sum(mtx(:));
%Convolution
imgB = imfilter(img, mtx, 'replicate');
%Downsampling
[r, c] = size(img);
imgOut = imgB(1:2:r, 1:2:c); %imresize(imgB, 0.5, 'bilinear');
end

（2）拉普拉斯金字塔

在高斯金字塔的运算过程中，图像经过卷积和下采样操作会丢失部分高频细节信息。为描述这些高频信息，人们定义了拉普拉斯金字塔(Laplacian Pyramid， LP)。用高斯金字塔的每一层图像减去其上一层图像的上采样并高斯卷积之后的预测图像，得到一系列的差值图像即为 LP 分解图像。

function p = pyrLapGen(img)
[r, c, ~] = size(img);
levels = floor(log2(min(r, c)));
list = [];
for i=1:levels
    %Calculating detail and base layers
    [tL0, tB0] = pyrLapGenAux(img);
    img = tL0;
    %Detail layer
    ts   = struct('detail', tB0);
    list = [list, ts];  
end
%Base layer
p = struct('list', list, 'base', tL0);
end

function [L0, B0] = pyrLapGenAux(img)
%5x5 Gaussian kernel
kernel = [1, 4, 6, 4, 1];
mtx = kernel' * kernel;
mtx = mtx / sum(mtx(:));
%Convolution
imgB = imfilter(img, mtx, 'replicate');
%Downsampling
[r, c] = size(img);
L0 =  imgB(1:2:r, 1:2:c); %imresize(imgB, 0.5, 'bilinear');
%Upsampling
imgE = imresize(L0, [r, c], 'bilinear');
%Difference between the two levels
B0 = img - imgE;
end

function imgOut = MertensTMO(imageStack, weights)
%default parameters if they are not present
if(~exist('weights', 'var'))
    weights = ones(1, 3);
end
wE = weights(1);
wS = weights(2);
wC = weights(3);

%number of images in the stack
[r, c, col, n] = size(imageStack);

%compute the weights for each image
total  = zeros(r, c);
weight = ones(r, c, n);
for i=1:n
    if(wE > 0.0)
        weightE = MertensWellExposedness(imageStack(:,:,:,i));
        weight(:,:,i) = weight(:,:,i) .* weightE.^wE;
    end
    if(wC > 0.0)
        if(size(imageStack(:,:,:,i), 3) > 1)
            L = mean(imageStack(:,:,:,i), 3);
        else
            L = imageStack(:,:,:,i);
        end
        weightC = MertensContrast(L);
        weight(:,:,i) = weight(:,:,i) .* (weightC.^wC);
    end
    if(wS > 0.0)
        weightS = MertensSaturation(imageStack(:,:,:,i));
        weight(:,:,i) = weight(:,:,i) .* (weightS.^wS);
    end
    weight(:,:,i) = weight(:,:,i) + 1e-12;
    total = total + weight(:,:,i);
end
for i=1:n %weights normalization
    weight(:,:,i) = weight(:,:,i) ./ total;
end

%empty pyramid
pyrAcc = [];
for i=1:n
    %Laplacian pyramid: image
    pyrImg = pyrImg3(imageStack(:,:,:,i), @pyrLapGen);
    %Gaussian pyramid: weight   
    pyrW   = pyrGaussGen(weight(:,:,i));
    %image times weight
    pyrImgW = pyrLstS2OP(pyrImg, pyrW, @pyrMul);
    if(i == 1)
        pyrAcc = pyrImgW;
    else %accumulate
        pyrAcc = pyrLst2OP(pyrAcc, pyrImgW, @pyrAdd);
    end
end

%reconstruction
imgOut = zeros(r, c, col);
for i=1:col
    imgOut(:,:,i) = pyrVal(pyrAcc(i));
end

%clamp to values in [0,1]
min_i = min(imgOut(:));
max_i = max(imgOut(:));
imgOut = ClampImg((imgOut - min_i) / (max_i - min_i),  0.0, 1.0);
end

饱和度权重：

function Ws = MertensSaturation(img)
[r,c,col] = size(img);
if(col==1)
    Ws = ones(r, c);
else
    mu = zeros(r,c);
    for i=1:col
        mu = mu + img(:,:,i);
    end
    mu = mu/col;
    sumC = zeros(r,c);
    for i=1:col
        sumC = sumC + (img(:,:,i)-mu).^2;
    end
    Ws = sqrt( sumC/col );
end
end

对比度权重：

function Wc = MertensContrast(L)
H = [0 1 0; 1 -4 1; 0 1 0];
Wc = abs(imfilter(L, H, 'replicate'));
End

优曝光权重：

function We = MertensWellExposedness(img, we_mean, we_sigma)
%mean for the Well-exposedness weights.
if(~exist('we_mean', 'var'))
    we_mean = 0.5; %as in the original paper
end
%sigma for the Well-exposedness weights.
if(~exist('we_sigma', 'var'))
    we_sigma = 0.2; %as in the original paper
end
sigma2 = 2.0 * we_sigma^2;
[r, c, col] = size(img);
We = ones(r, c);
for i=1:col
    We = We .* exp(-(img(:,:,i) - we_mean).^2 / sigma2);
end
end

金字塔对应层相乘：

function pOut = pyrMul(pA, pB)
%checking lenght of the pyramids
nA = length(pA.list);
nB = length(pB.list);
%multiplying base levels
pOut.base = pA.base .* pB.base;
pOut.list = pA.list;
%multiplying the detail of each level
for i=1:nA
    pOut.list(i).detail = pA.list(i).detail .* pB.list(i).detail;
end
end

金字塔图像重构融合：

function img = pyrVal(pyramid)
list = pyramid.list;
base = pyramid.base;
 
n = length(list);
img=[];
for i=1:n
    ind = n - i + 1;
    [r, c] = size(list(ind).detail);
    if(i == 1)        
        base = imresize(base, [r, c], 'bilinear');
        img  = base + list(ind).detail;
    else
        img = imresize(img, [r, c], 'bilinear');
        img = img + list(ind).detail;        
    end
end
end