HDR_rec.m

function out = HDR_rec(IM1,IM2,IM3,alpha,snr)
% Authors : Abhiram Gnanasambandam (agnanasa@purdue.edu),
%Stanley H. Chan (stanchan@purdue.edu)
%Intelligent Imaging Lab, Dept. of ECE, Purdue University
% Reference:
%[1]HDR Imaging with Quanta Image Sensors: Theoretical Limits and Optimal
%Reconstruction, Abhiram Gnanasambandam and Stanley H. Chan, 
%IEEE Transactions on Computational Imaging, 2020
%[2]High Dynamic Range Imaging using Quanta Image Sensors, Abhiram
%Gnansambandam, Jiaju Ma, Stanley H. Chan, IISW 2019.
% This function reconstructs HDR QIS images
%IM1,IM2,IM3 - 3 input images 
%alpha - Integration times corresponding to these 3 images
%snr - SNRH lookuptable
IM1(IM1<0) = 0;
IM2(IM2<0) = 0;
IM3(IM3<0) = 0;  
im1 = IM1/alpha(1);
im2 = IM2/alpha(2);
im3 = IM3/alpha(3);


wt1 = 1/3 * ones(size(im1));
wt2 = 1/3 * ones(size(im2));
wt3 = 1/3 * ones(size(im3));


for itr = 1:10
    rec = (im1.*wt1 + im2.*wt2 + im3.*wt3 ); 

    wt1 = snr(alpha(1)*rec);
    wt1(wt1<0) = 0;
    wt1 = wt1.^2;

    wt2 = snr(alpha(2)*rec);
    wt2(wt2<0) = 0;
    wt2 = wt2.^2;

    wt3 = snr(alpha(3)*rec);
    wt3(wt3<0) = 0;
    wt3 = wt3.^2;


    swt = wt1 + wt2 + wt3;


    wt1 = wt1./(swt);
    wt2 = wt2./(swt);
    wt3 = wt3./(swt);

    wt1(swt == 0) = 1/3;
    wt2(swt == 0) = 1/3;
    wt3(swt == 0) = 1/3;

end
out = reshape(rec,size(IM1));


end