Reinhard的HDR色调映射算法参数自动选择版本MatLab程序

Reinhard的HDR色调映射算法参数自动选择版本的具体原理及公式请参见：Erik Reinhard (2002) Parameter Estimation for Photographic Tone Reproduction, Journal of Graphics Tools, 7:1, 45-51, DOI: 10.1080/10867651.2002.10487554。这里我就不再赘述了。下面直接给出代码及其处理效果（算法会根据Zone数，即动态范围，自动选择是用全局映射算法还是局部映射算法）：

%% HDR算法复现
% based on the tone mapping technology as described in Erik Reinhard, Michael Stark, 
% Peter Shirley, James Ferwerda, “Photographic Tone Reproduction for Digital Images”
% ACM, 2002
% and Reinhard's “Parameter Estimation for Photographic Tone Reproduction”
% Journal of Graphics Tools, 7:1, 45-51, DOI: 10.1080/10867651.2002.10487554
% author: Zhuangzhuang Zhang
% time: 2020/10/15
close all;clc
filename = 'data/HDR images2/nave.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/1.6;
figure,imshow(hdr.^gamma);
imwrite(hdr.^gamma,'原始HDR图像.png');
%%  预处理，log域均值取指数结果映射到中性灰
Lw = 0.27*hdr(:,:,1) + 0.67*hdr(:,:,2) + 0.06*hdr(:,:,3) + 1e-6; %论文推荐的world luminance计算公式,小偏移量为避免被0除的情况
R = hdr(:,:,1) ./ Lw;
G = hdr(:,:,2) ./ Lw;
B = hdr(:,:,3) ./ Lw;

% maxL=max(L(:));
% settedMaxL=20;
% thL=10;
% mask=L>thL;
% settedMaxL=min(maxL,settedMaxL);
% L(mask)=(settedMaxL-thL)*(L(mask)-thL)/(maxL-thL)+thL;%高曝光部分压缩动态范围

LL=log(Lw+10^-9);
meanLw=exp(mean(LL(:))); %实践中发现将对数均值再反对数变换（取指数）结果映射到key值，通常可以得到更好的效果（主要是大区域之间的对比度更加合理而不会过低）
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);
L=Lw*a/meanLw;
if zone<11
    Lwhite=1.5*2^(log2(maxLw)-log2(minLw)-5);
    Ld=(L+L.^2/Lwhite^2)./(L+1);
else
    picNum=8;
    ratio=1.3;
    phi=8.0;
    eps=0.05;
    alpha=1/sqrt(8);
    refV=1.0;
    V1=zeros([size(L),picNum+1]); % V2的当前尺度与V1的下一尺度相同
    %-----------------------空域实现-------------------------------%
    tic
    for i=1:picNum+1 % picNum+1个尺度，以ratio为倍数递增
        s=ratio^(i-1);
        sigma=alpha*s;
        kernelRadius = ceil(2*sigma);
        kernelSize = 2*kernelRadius+1;
        gaussKernelHorizontal = fspecial('gaussian', [kernelSize 1], sigma);
        temp= conv2(L, gaussKernelHorizontal, 'same');
        gaussKernelVertical = fspecial('gaussian', [1 kernelSize], sigma);
        V1(:,:,i)= conv2(temp, gaussKernelVertical, 'same');
    end
    V=zeros([size(L),picNum]);
    for i=1:picNum
        s=ratio^(i-1);
        V(:,:,i)=abs(V1(:,:,i+1)-V1(:,:,i))./(a*2^phi/s^2+V1(:,:,i));
    end
    
    Vsm=V1(:,:,picNum); %给Vsm初始化，作为8个尺度都不满足条件下的预设值
    for i=1:size(L,1)
        for j=1:size(L,2)
            for k=1:picNum
                if V(i,j,k)>eps
                    Vsm(i,j)=V1(i,j,k);
                    break;
                end
            end
        end
    end
    
    toc
    figure,imshow(Vsm);
    imwrite(Vsm,'各点参考值图像.png');
    Ld=L./(refV+Vsm);
end
Ld(Ld>1)=1;
figure,imshow(Ld);
%% 色彩还原
rgb=zeros(size(L,1),size(L,2),3);
rgb(:,:,1)=Ld.*R;
rgb(:,:,2)=Ld.*G;
rgb(:,:,3)=Ld.*B;
rgb=rgb.^gamma;
figure,imshow(rgb);

输入图像：

输出（选择了局部色调映射）：

输入：

输出（选择了全局色调映射）：

可以看出，参数自动选择算法所选的参数和映射策略表现值得肯定。