var_sigma.m

function x = var_sigma(L, theta,p)
% 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.
x0 = L^2;
tmp = 0;
for q = 0:L-1
    tmp = tmp + (2*q+1)*gammainc(theta,q+1,'upper');
end

x0 = x0 - tmp;

del_x = 0;

for m = -1:1
    tmp = 0;
    for q = 0:L-1
        tmp = tmp + q^2*(exp(-theta)*theta^(max(q-m,0)) /factorial(max(q-m,0)) - exp(-theta) * theta^q/factorial(q)) ;
    end
    tmp = tmp + L^2*(gammainc(theta,L,'upper') - gammainc(theta,max(L-m,0),'upper'));
    del_x = del_x + tmp*p(m+2);
end

x = x0 + del_x;

x = x - mean_mu(L, theta,p).^2;
%x(x<0) = 0;
end