single_bit_recon.m

function out = single_bit_recon(frames, K ,T, q, alpha, denoiser)
% Authors : Abhiram Gnanasambandam (agnanasa@purdue.edu),
%Stanley H. Chan (stanchan@purdue.edu)
%Intelligent Imaging Lab, Dept. of ECE, Purdue University
% Reference:
%[1]Stanley H. Chan, Omar Elgendy and Xiran Wang, ‘‘Images from bits: 
%Non-iterative image reconstruction for quanta image sensors’’, 
%MDPI Sensors Special Issue on Photon-Counting Image Sensors, 
%vol. 16, no. 11, paper 1961, pp.1-21, Nov. 2016.

% frames [T,H,W]- Noisy spatio-temporal data - Should have T frames
% K - Spatial oversapling coefficient
% T - Number of frames
% q - ADC threshold for the binary data
% alpha - The brightness adjuster
% denoiser - Denoiser to be used for the reconstruction 
%           Any standard Gaussian denoiser should work. Should take image
%           in range [0,1] as input and noise level sigma in [0,1].
addpath('utils/utilities');

L = K^2*T;%Total Oversampling 
over_sampled_image = squeeze(blockfun(frames,[T,K,K],@sum));

var_normalized = sqrt(L+0.5)*asin(sqrt( (over_sampled_image + 3/8)/(L+3/4) ));% Anscombe for binomial data
sigma = 1/2;

maxv =  max(var_normalized(:));
minv =  min(var_normalized(:));

brightness_normalized = (var_normalized - minv)/(maxv - minv); %Make sure data in [0,1]
sigma_norm = 1/2/(maxv - minv);

denoised1 = denoiser(brightness_normalized,sigma_norm);

denoised2 = denoised1 * (maxv - minv) +  minv;

denoised3 = 1/(1+1/2/L) * ( (L+3/4) * sin(denoised2/(sqrt(L+1/2))).^2-1/8);

out = K/alpha * gammaincinv(1 - denoised3/L,q,'upper');

end