Poisson Non-Local PCA

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=1).

Authors: J. Salmon, C.-A. Deledalle, R. Willett, Z. Harmany, ICASSP 2012, PDF.
Matlab toolbox for 2D images NLPCA_code.ZIP (fast v2: NLPCA_v2_code.ZIP thanks to Anthony Wang and Albert Oh). Matlab toolbox for hyperspectral images NLPCA_hyperspectral_code.ZIP (fast v2: NLPCA_hyperspectral_v2_code.ZIP thanks to Anthony Wang and Albert Oh).
Copyright (C) 2012 NL-PCA project. See The GNU Public License (GPL)

Initialization

clear all
close all

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

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=1;     %0/1 to activate a double iteration of the algorithm
anscombe=0;             % Direct Poisson (0) and  Anscombe + Gaussian (1)

Loading 'Saturn' image

ima_ori=double(imread('./data/saturn.tif'));
ima_ori= ima_ori((1:256)+70,(1:256));

Noisy image generation

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);

    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);
end
toc

Elapsed time is 165.459557 seconds.

Result display

figure('Position',[100 100   1200 600])
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);

Copyright ® 2012