基于视觉显著性度量的HDR图像色调映射MatLab程序（论文复现系列）

最新推荐文章于 2023-03-27 15:39:42 发布

光电学子

最新推荐文章于 2023-03-27 15:39:42 发布

阅读量2.3k

点赞数 5

分类专栏：信号与图像处理

本文链接：https://blog.csdn.net/u014230360/article/details/109157891

版权

信号与图像处理专栏收录该内容

45 篇文章 62 订阅

订阅专栏

该文主要是提供论文"Visual-Salience-Based Tone Mapping for High Dynamic Range Images"中基于视觉显著度的HDR色调映射算法的MatLab实现代码。算法的基本思想基本思想并不复杂：首先将HDR图像分解为背景层和细节层，对背景层进行亮度压缩，对细节层进行对比度增强，然后重新组合形成映射后的亮度图，接着进行色彩重现完成整个处理。所不同的是使用了保边效果更好的加权引导滤波对HDR图像进行分解以抑制光晕效应，并将显著度嵌入到背景层的亮度压缩中。论文还提出将显著度也嵌入到HDR图像的分解中（加权引导滤波的权重），提出的权重公式为：

其中分子分母分别为edge-aware的权重与背景层saliency-aware的权重。这里背景层saliency-aware的权重为什么会出现在分母中论文中没讲，我也没想到有什么合理的解释（或许这一步本身就缺乏合理性）。由于计算 $\Gamma_b(p')$ 的过程论文中语焉不详（使用什么算子计算梯度及其方向，梯度的两个方向分量均为0时方向如何认定等都没有说明），这部分也难以准确复现。我个人使用简单一阶差分计算GX、GY及atan(GY./GX)进行方向角度计算的复现结果并不理想，所以权重部分我使用的就只是 $W(p')=\Gamma_e(p')$ 。另外还尝试了使用 $W(p')=\Gamma_e(p')*\Gamma_b(p')$ ，即认为显著区域及边界区域都应该给更大的权重。

论文中关于量化的部分也是语焉不详，量化等级K的值也没给出。我的实现中是将对数域亮度进行线性映射，最小值和最大值分别映射到1和K，而不是设定具体的量化区间，且只针对亮度通道进行显著性度量的计算（分通道计算的结果不如只在亮度通道计算的效果理想）。

下面给出代码及复现结果（代码在最下方）。

原图（伽马校正后）：

权重为 $W(p')=\Gamma_e(p')$ 时分解得的背景层图像及细节层图像：

权重为 $W(p')=\Gamma_e(p')*\Gamma_b(p')$ 时分解得的背景层图像及细节层图像：

两种权重的最终色调映射结果分别为：

可以看出，将显著性度量作为权重项中的乘数因子使得映射结果中显著区域更加突出。

作为参考，Reinhard色调映射算法(参数自动选择版本)的处理结果为：

同样地，对于下面的输入图像：

两种权重的映射结果分别为：

作为参考，Reinhard色调映射算法(参数自动选择版本)的处理结果为：

输入图像：

两种权重的处理结果：

作为参考，Reinhard色调映射算法(参数自动选择版本)的处理结果为：

从前面的输入图像的处理结果可以看出，一般来说，论文中权重为 $W(p')=\Gamma_e(p')$ 时的映射结果与Reinhard色调映射算法处理结果的结果比较接近，不太好分出优劣，而 $W(p')=\Gamma_e(p')*\Gamma_b(p')$ 时的映射结果通常亮度偏暗或者偏暗，因此，建议一般情况下根据期望的输出图像亮度水平来选择Reinhard色调映射算法或者 $W(p')=\Gamma_e(p')$ 情况下的的论文算法来实施色调映射。

接下来我们再考察一下加权引导滤波和不加权的引导滤波之间的差异：

对于前面的第一张输入图像，加权与不加权引导滤波分解出的背景层图像分别为：

动态范围压缩后分别为：

可以看到，加权引导滤波具有更好的保边性能，背景层动态范围压缩后对比度更高，边缘更加突出。

加权及不加权引导滤波分解的色调映射处理结果分别为：

可以看出，由于加权引导滤波的背景层动态范围压缩得更好，重新组合后的映射结果同样具有更加优秀的动态范围压缩效果，因为整体处理质量看起来更好一些。

以下是论文算法的MatLab实现代码（Reinhard色调映射算法的参数自动选择版本的程序见下一篇博客）：

% "Visual-Salience-Based Tone Mapping for High Dynamic Range Images"论文算法复现
% author: Zhuangzhuang Zhang
% 2020/10/19
% clear all;
close all;clc
filename = 'data/HDR images2/rosette.hdr';
try
    hdr = double(hdrread(filename));
catch
    disp('Can not open file!');
    return;
end
if max(size(hdr))>1500
    hdr=imresize(hdr,0.5);
    hdr(hdr<0)=0;
end

gamma=1/2.0;
figure,imshow(hdr.^gamma);
imwrite(hdr.^gamma,'原始HDR图像.png');
%%  预处理，log域均值取指数结果映射到中性灰
Lw = 0.299*hdr(:,:,1) + 0.587*hdr(:,:,2) + 0.144*hdr(:,:,3) + 1e-9; %论文推荐的world luminance计算公式,小偏移量为避免被0除的情况
R = hdr(:,:,1) ./ Lw;
G = hdr(:,:,2) ./ Lw;
B = hdr(:,:,3) ./ Lw;

%% 具有边缘意识的色调映射
Lh=log(Lw);
meanLw=exp(mean(Lh(:))); 
maxLw=quantile(Lw(:),0.99);
minLw=quantile(Lw(:),0.01);
zone=log2(maxLw)-log2(minLw);
a=0.18*4^((2*log2(meanLw)-log2(minLw)-log2(maxLw))/zone);
r=15;
eps=1;
Ge=varBasedWeight(Lh); %\Gamma_e
Lb=WGIF(Lh,Lh,r,eps,Ge); %加权引导滤波获得基底层
Ld=Lh-Lb; %细节层
figure,imshow(exp(Lb));
figure,imshow(exp(Ld));
meanLh=mean(Lh(:));
compLb=Lb+log(a)-meanLh-log(1+a*exp(Lb-meanLh));
figure,imshow(exp(compLb));
theta=1.5; %细节缩放因子
Lt=exp(compLb+theta*Ld);
figure,imshow(Lt);
% 色彩还原
rgb=zeros(size(Lw,1),size(Lw,2),3);
rgb(:,:,1)=Lt.*R;
rgb(:,:,2)=Lt.*G;
rgb(:,:,3)=Lt.*B;
rgb=rgb.^gamma;
figure,imshow(rgb);

%% 基于显著度的色调映射与色彩还原
Gb=SaliencyBasedICH(Lh,20,4); %\Gamma_b
temp=1./(Gb+10^-9);
Gb=(Gb+10^-9)*mean(temp(:));
lambda=0.75;
Gb=Gb.^lambda;
Gb(Gb<1)=1;
figure,imshow(Gb/max(Gb(:)));
W2=Ge.*Gb;
Lb2=WGIF(Lh,Lh,r,eps,W2); %加权引导滤波获得基底层
Ld2=Lh-Lb2; %细节层
figure,imshow(exp(Lb2));
figure,imshow(exp(Ld2));

wLh=Gb.*Lh;
meanLh2=sum(wLh(:))/sum(Gb(:));
compLb2=Lb2+log(a)-meanLh2-log(1+a*exp(Lb2-meanLh2));
theta=1.5; %细节缩放因子
Lt2=exp(compLb2+theta*Ld2);
rgb2=zeros(size(Lw,1),size(Lw,2),3);
rgb2(:,:,1)=Lt2.*R;
rgb2(:,:,2)=Lt2.*G;
rgb2(:,:,3)=Lt2.*B;
rgb2=rgb2.^gamma;
figure,imshow(rgb2);

function wM = varBasedWeight(img,phi)
if ~exist('phi','var')
    phi=0.75;
end
img=img-min(img(:));
L=max(img(:))-min(img(:));
r=15;
paddedImg=padarray(img,[r r],'replicate','both');
sq_paddedImg=paddedImg.^2;
intMap=intergalMap(paddedImg);
sq_intMap=intergalMap(sq_paddedImg);
sumImg=zeros(size(img));
meanImg=zeros(size(img));
sq_sumImg=zeros(size(img));
varImg=zeros(size(img));
wM=zeros(size(img)); %edge-aware weight image
num=(2*r+1)^2;
for i=1:size(img,1)
    for j=1:size(img,2)
        sumImg(i,j)=intMap(i+2*r+1,j+2*r+1)+intMap(i,j)-intMap(i,j+2*r+1)-intMap(i+2*r+1,j);
        meanImg(i,j)=sumImg(i,j)/num;
        sq_sumImg(i,j)=sq_intMap(i+2*r+1,j+2*r+1)+sq_intMap(i,j)-sq_intMap(i,j+2*r+1)-sq_intMap(i+2*r+1,j);
        varImg(i,j)= sq_sumImg(i,j)-(sumImg(i,j)^2)/num;
    end
end
v1=(0.001*L)^2;
v2=10^-9;
wM=((varImg+v1)./(meanImg.^2+v2)).^phi;
temp=1./wM;
wM=wM*mean(temp(:));


function q = WGIF(I, p, r, eps, W)
%   WEIGHTEDGUIDEDFILTER   O(1) time implementation of weighted guided filter.
%
%   - guidance image: I (should be a gray-scale/single channel image)
%   - filtering input image: p (should be a gray-scale/single channel image)
%   - local window radius: r
%   - regularization parameter: eps
%   - weight image: W 

[hei, wid] = size(I);
N = boxfilter(ones(hei, wid), r); % the size of each local patch; N=(2r+1)^2 except for boundary pixels.

mean_I = boxfilter(I, r) ./ N;
mean_p = boxfilter(p, r) ./ N;
mean_Ip = boxfilter(I.*p, r) ./ N;
cov_Ip = mean_Ip - mean_I .* mean_p; % this is the covariance of (I, p) in each local patch.

mean_II = boxfilter(I.*I, r) ./ N;
var_I = mean_II - mean_I .* mean_I;

a = W.*cov_Ip ./ (W.*var_I + eps); % Eqn. (5) in the paper;
b = mean_p - a .* mean_I; % Eqn. (6) in the paper;

mean_a = boxfilter(a, r) ./ N;
mean_b = boxfilter(b, r) ./ N;

q = mean_a .* I + mean_b; % Eqn. (8) in the paper;
end

function S=SaliencyBasedICH(img,K,r)
%基于灰度共生矩阵的显著度矩阵计算程序
%K为量化等级
%r为衡量共生性质的邻域半径，灰度m,n在(2r+1,2r+1)的窗口内同时出现过即视为共生
M=ICH(img,K,r);
nM=M/sum(M(:));
nzMask=nM~=0;
nzMat=double(nzMask);
nzNum=sum(nzMat(:));
invPMF=nM;
mask=nzMask & nM<1/nzNum;
invPMF(mask)=1/nzNum-nM(mask);
invPMF(~mask)=0;

%求窗口w*w内的平均显著度
interval=(max(img(:))-min(img(:)))/(K-0.5);
Q=floor((img-min(img(:)))/interval)+1;
pQ=padarray(Q,[r r],'replicate','both');
num=(2*r+1)^2;
S=zeros(size(img,1),size(img,2));
for i=r+1:size(pQ,1)-r
    for j=r+1:size(pQ,2)-r
        sumS=0;
        for p=-r:r
            for q=-r:r
                sumS=sumS+invPMF(pQ(i,j),pQ(i+p,j+q));
            end
        end
        S(i-r,j-r)=sumS/num;
    end
end
S=S/max(S(:));

function M=ICH(img,K,r)
%求灰度共生矩阵程序
%K为量化等级
%r为衡量共生性质的邻域半径，灰度m,n在(2r+1,2r+1)的窗口内同时出现过即视为共生
if ~exist('K','var')
    K=20;
end
if ~exist('r','var')
    r=4;
end
interval=(max(img(:))-min(img(:)))/(K-0.5);
Q=floor((img-min(img(:)))/interval)+1;
rowNum=size(Q,1);
M=zeros(K,K);
idxShift=zeros(r+1,r+1); %参考像素相对于当前像素的索引偏移量
for i=0:r
    for j=0:r
        idxShift(i+1,j+1)=i+j*rowNum;
    end
end
idxShift=idxShift(:);
for p=2:length(idxShift) %idxShift(1)为0，表示当前进行统计的两像素实际为同一像素，跳过
    for i=1:size(Q,1)-r
        for j=1:size(Q,2)-r
            idx=(j-1)*rowNum+i;
            s=idxShift(p);
            M(Q(idx),Q(idx+s))=M(Q(idx),Q(idx+s))+1;
            M(Q(idx+s),Q(idx))=M(Q(idx+s),Q(idx))+1;
        end
    end
end

，

光电学子

关注

5
点赞
踩
39

收藏

觉得还不错? 一键收藏
9
评论
基于视觉显著性度量的HDR图像色调映射MatLab程序（论文复现系列）

该文主要是提供论文"Visual-Salience-Based Tone Mapping for High Dynamic Range Images"的MatLab
复制链接

扫一扫