QIS_color_non_iterative.m

%Authors : Abhiram Gnanasambandam (agnanasa@purdue.edu), 
%Omar Elgendy (oelgendy@gigajot.tech),Stanley H. Chan (stanchan@purdue.edu)
%Intelligent Imaging Lab, Dept. of ECE, Purdue University
%References :
%[1] Megapixel photon-counting color imaging using quanta image sensor,
% Abhiram Gnanasambandam, Omar A. Elgendy, Jiaju Ma, and Stanley H. Chan,
% OSA Optics Express
%This code demonstrates the color reconstruction method described in the
%above mentioned paper.

clc,clear,close all
rng('default')
addpath('utils/denoisers')
addpath('ADMM_files')
addpath('Simulation/Matlab')
%%

signalLevel =  2;
sigmaRN     =  0.25;
%Number of frames
T       =  1;
%Number of bits
Nbits = 5;

%%
xGT = im2double(imread('data/lena.jpg'));
%xGT = imresize(xGT,[128,128]);

[H,W,~] = size(xGT);
if(H>W)
    xGT = permute(xGT,[2,1,3]);
    transposed = 1;
end
[H,W,~] = size(xGT);
H1=2^ceil(log2(H));
W1=2^ceil(log2(W));
xGT = imresize(xGT,[H1,W1]);


hr      =  [0 0;0 1];
hg      =  [0 1;1 0];
hb      =  [1 0;0 0];
maskR   =  repmat(hr,H1/2,W1/2);
maskG   =  repmat(hg,H1/2,W1/2);
maskB   =  repmat(hb,H1/2,W1/2);
mask    =  cat(3,maskR,maskG,maskB);
[xL,xC1,xC2] = imsplit(RGB2LC1C2(xGT));
xGT = LC1C22RGB(cat(3,xL,xC1,xC2));
hr  = repmat([0 0;0 1],H1/2,W1/2);
hg  = repmat([0 1;1 0],H1/2,W1/2);
hb  = repmat([1 0;0 0],H1/2,W1/2);
xCFA= sum(xGT.*cat(3,hr,hg,hb),3);
xGT255 = xGT*255.0;




%%

peak   =  signalLevel/mean(xCFA(:));
x      =  squeeze(sum(simulate_noisy_images_lite(xCFA*peak,sigmaRN,T,Nbits),1)); %poissrnd(xCFA*peak) +  sigmaRN*randn(H1,W1);
%%
% Anscombe Transform for stabilizing the noise variance
% Would be generally appropriate to replace the Poisson VST with the
% corresponding VST for respective bet depth, which is outside the scope of
% this paper.
Z        =  2*sqrt(max(3/8 + x + T*sigmaRN^2,0));

sigmaVST =  1;  % Standard deviation after Anscombe Transform


Bn = Z(1:2:end,1:2:end);
Gn = Z(2:2:end,1:2:end) + Z(1:2:end,2:2:end);
Rn = Z(2:2:end,2:2:end);

maxr = max(Rn(:));
minr = min(Rn(:));
Rn = (Rn - minr)/(maxr - minr);
sigmar = sigmaVST/(maxr-minr);

Rd = wrapper_GBM3D(Rn,sigmar) * (maxr-minr) + minr;

maxg = max(Gn(:));
ming = min(Gn(:));
Gn = (Gn - ming)/(maxg - ming);
sigmag = sigmaVST*sqrt(2)/(maxg-ming);

Gd = wrapper_GBM3D(Gn,sigmag) * (maxg-ming) + ming;

maxb = max(Bn(:));
minb = min(Bn(:));
Bn = (Bn - minb)/(maxb - minb);
sigmab = sigmaVST/(maxb-minb);

Bd = wrapper_GBM3D(Bn,sigmar) * (maxb-minb) + minb;
Z = zeros(size(Bd,1),size(Bd,2),3);

Z(:,:,1) = Rd;
Z(:,:,2) = Gd/2;
Z(:,:,3) = Bd;

Z = (Z/2).^2 - 1/8 - T*sigmaRN^2;


%%
figure;
subplot(1,2,1);imshow(x/T/peak);title('Noisy input');
subplot(1,2,2);imshow(Z/peak/T); title('Color Reconstruction');