function [y_est, y_hat] = denoiseVSTpoisson (z)
    % Iterative Poisson denoising algorithm based on variance stabilization and AWGN denoising (v1.00)
%
% INPUT:   z   Poisson noisy image   z ~ P(y)
% OUTPUT: yhat  denoised image,  estimate of y
%
%
% L. Azzari and A. Foi, "Variance Stabilization for Noisy+Estimate
% Combination in Iterative Poisson Denoising", submitted, March 2016
%
% http://www.cs.tut.fi/~foi/invansc/
%
%  L. Azzari and Alessandro Foi - Tampere University of Technology - 2016 - All rights reserved.
% -----------------------------------------------------------------------------------------------
%

%% Prepare functions
load Anscombe_lambda lambdaGridTimesE yGrid lambdaGrid   % load expectation precomputed on a grid for interpolation of the exact unbiased inverse

% asymInv=@(x) (x.^2/4-1/8);                      % asymptotically unbiased inverse of Anscombe (coincides with exact unbiased for large y)
f=@(x,lambda) 2*sqrt(max(0,x/lambda^2+3/8));      % forward VST (Anscombe of scaled input)

%% Apply forward VST
fz = f(z,1);

%% AWGN DENOISING

sigma_den = 1;  % Standard-deviation value assumed after variance stabiliation

% Scale the image (BM3D processes inputs in [0,1] range)
scale_range = .9;
scale_shift = (1-scale_range)/2;
maxzans = max(fz(:));
minzans = min(fz(:));
fz   =   (fz-minzans)/(maxzans-minzans);
sigma_den = sigma_den/(maxzans-minzans);
fz = scale_shift + fz*scale_range;
sigma_den = sigma_den*scale_range;

[dummy, y_est, y_hat] = BM3D(0,fz,sigma_den*255,'np',0); % denoise assuming AWGN using BM3D http://www.cs.tut.fi/~foi/GCF-BM3D/

% Scale back to the initial VST range
y_est = (y_est-scale_shift)/scale_range;
y_est = y_est*(maxzans-minzans)+minzans;

y_hat = (y_hat-scale_shift)/scale_range;
y_hat = y_hat*(maxzans-minzans)+minzans;

%% Apply the inverse VST for convex combination z_i of Poisson z and estimate yhat
y_est = Anscombe_inverse_exact_unbiased (y_est);   % exact unbiased inverse of pure Poisson  http://www.cs.tut.fi/~foi/invansc/
y_hat = Anscombe_inverse_exact_unbiased (y_hat);   % exact unbiased inverse of pure Poisson  http://www.cs.tut.fi/~foi/invansc/
end