DEMO_NLPCA.m

%% Simple example on Saturn
%
% We present here the NL-PCA algorithm to denoise Poisson corrupted
% images.
% Two versions are proposed. The default one is the one dealing directly
% with the Poisson structure (anscombe=0). The second one is the one
% performing a variance stabilization step/ Anscombe transform 
% (anscombe=0).
% 
%   Copyright (C) 2012 NL-PCA project
%   Joseph Salmon, Charles-Alban Deledalle, Rebecca Willet, Zachary Harmany
%
%   See The GNU Public License (GPL)

%---------------------------------------------------------------------
%
%   This file is part of NL-PCA.
%
%   NL-PCA is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as
%   published by the Free Software Foundation, either version 3 of
%   the License, or (at your option) any later version.
%
%   NL-PCA is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public
%   License along with NL-PCA.  If not, see
%   <http://www.gnu.org/licenses/>.
%   Joseph Salmon, Charles-Alban Deledalle, Rebecca Willet, Zachary Harmany
%
%   See The GNU Public License (GPL)


clear all
close all

addpath('functions')
addpath('tools')
addpath('Anscombe')

tic
%% Parameters:
Patch_width=20;         % Patch width 
nb_axis=4;              % Number of axis in the "PCA" 
nb_iterations=20;       % Max number of iteration 
nb_clusters=14;         % Number of clusters in the kmeans step
eps_stop=1e-1;          % Stoping criterion
epsilon_cond=1e-3;      % Condition number for Hessian inversion
double_iteration=0;     %0/1 to activate a double iteration of the algorithm
anscombe=0;             % Direct Poisson (0) and  Anscombe + Gaussian (1)
newton=0;               % Newton's method (1) and no Newton's method (0)
parallel=1;             % Parallelization (1) and no parallelization (0)

%% Loading image part
%  'Saturn'
ima_ori=double(imread('./data/saturn.tif'));
ima_ori= ima_ori((1:256)+70,(1:256));

%% Noisy image generation for the simulation case
peak=0.1;
sd=1;
rng(sd)
Q = max(max(ima_ori)) /peak;
ima_lambda = ima_ori / Q;
ima_lambda(ima_lambda == 0) = min(min(ima_lambda(ima_lambda > 0)));
ima_nse_poiss = knuth_poissrnd(ima_lambda);
[m,n]=size(ima_nse_poiss);

func_clustering=@(X) clustering_litekmeans(X,Patch_width,nb_clusters,m,n);
func_thresholding = @(ima_ppca) no_thresholding(ima_ppca);

%% Denoising Part
tic
if anscombe==1
    
    eps_stop=1e-3;
    epsilon_cond=1e-5;
    func_denoising_patches=@(X)...
        gaussian_NL_PCA(X{1},nb_axis,nb_iterations,...
        X{2},X{3},eps_stop,epsilon_cond);    
    func_recontruction=@(X) reconstruction_gaussian(X);    
    ima_nse_poiss_anscombe = 2*sqrt(ima_nse_poiss + 3/8);
    [ima_fil,ima_int,~,~]=NL_PCA(ima_nse_poiss_anscombe,...
        Patch_width,nb_axis,nb_clusters,func_thresholding,...
        func_recontruction,func_denoising_patches,func_clustering,...
        double_iteration,newton,parallel);
    
    ima_fil = Anscombe_inverse_exact_unbiased_Foi(ima_fil);
    ima_int = Anscombe_inverse_exact_unbiased_Foi(ima_int);    
else
    
    func_recontruction=@(X) reconstruction_poisson(X);
    func_denoising_patches=@(X)...
         poisson_NL_PCA(X{1},nb_axis,nb_iterations,X{2},X{3},...
         eps_stop,epsilon_cond);
    
    [ima_fil,ima_int,~,~]=NL_PCA(ima_nse_poiss,...
        Patch_width,nb_axis,nb_clusters,func_thresholding,...
        func_recontruction,func_denoising_patches,func_clustering,...
        double_iteration,newton,parallel);          
end
toc

%% Result display
figure
ax(1) = subplot(1, 3, 1);
plotimage(Q * ima_nse_poiss);
title(sprintf('Noisy PSNR = %f',psnr(Q*ima_nse_poiss, Q*ima_lambda, 255)));
ax(2) = subplot(1, 3, 2);
plotimage(Q * ima_int);
title(sprintf('First iteration = %f',psnr(Q*ima_int, Q*ima_lambda, 255)));
ax(3) = subplot(1, 3, 3);
plotimage(Q * ima_fil);
title(sprintf('Second iteration = %f',psnr(Q*ima_fil, Q*ima_lambda, 255)));
linkaxes(ax);