perform_blsgsm_denoising.m

A toolbox for the non-local means algorithm

function y = perform_blsgsm_denoising(x, options)

% perform_blsgsm_denoising - denoise an image using BLS-GSM
%
% y = perform_blsgsm_denoising(x, options);
%
%   BLS-GSM stands for "Bayesian Least Squares - Gaussian Scale Mixture". 
%
%   This function is a wrapper for the code of J.Portilla.
%
%   You can change the following optional parameters
%   options.Nor: Number of orientations (for X-Y separable wavelets it can
%       only be 3)
%   options.repres1: Type of pyramid ('uw': shift-invariant version of an orthogonal wavelet, in this case)
%   options.repres2: Type of wavelet (daubechies wavelet, order 2N, for 'daubN'; in this case, 'Haar')
%
%   Other parameter that should not be changed unless really needed.
%   options.blSize	    	n x n coefficient neighborhood (n must be odd): 
%   options.parent			including or not (1/0) in the
%       neighborhood a coefficient from the same spatial location
%   options.boundary        Boundary mirror extension, to avoid boundary artifacts 
%   options.covariance     	Full covariance matrix (1) or only diagonal elements (0).
%   options.optim          	Bayes Least Squares solution (1), or MAP-Wiener solution in two steps (0)
%
%   Default parameters favors speed.
%   To get the optimal results, one should use:
%       options.Nor = 8;                8 orientations
%       options.repres1 = 'fs';         Full Steerable Pyramid, 5 scales for 512x512
%       options.repres2 = '';           Dummy parameter when using repres1 = 'fs'   
%       options.parent = 1;             Include a parent in the neighborhood
%
%   J Portilla, V Strela, M Wainwright, E P Simoncelli, 
%   "Image Denoising using Scale Mixtures of Gaussians in the Wavelet Domain,"
%   IEEE Transactions on Image Processing. vol 12, no. 11, pp. 1338-1351,
%   November 2003 
%
%   Javier Portilla, Universidad de Granada, Spain
%   Adapted by Gabriel Peyre in 2007

options.null = 0;

if isfield(options, 'sigma')
    sig = options.sigma;
else
    sig = perform_noise_estimation(x)
end


[Ny,Nx] = size(x);
% Pyramidal representation parameters
Nsc = ceil(log2(min(Ny,Nx)) - 4);  % Number of scales (adapted to the image size)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% parameters
if isfield(options, 'Nor')
    Nor = options.Nor;
else
    Nor = 3;				            % Number of orientations (for X-Y separable wavelets it can only be 3)
end
if isfield(options, 'repres1')
    repres1 = options.repres1;
else
    repres1 = 'uw';                     % Type of pyramid (shift-invariant version of an orthogonal wavelet, in this case)
end

if isfield(options, 'repres2')
    repres2 = options.repres2;
else
    repres2 = 'daub1';                  % Type of wavelet (daubechies wavelet, order 2N, for 'daubN'; in this case, 'Haar')
end

% Model parameters (optimized: do not change them unless you are an advanced user with a deep understanding of the theory)
if isfield(options, 'blSize')
    blSize = options.blSize;
else
    blSize = [3 3];	    % n x n coefficient neighborhood of spatial neighbors within the same subband
end

if isfield(options, 'parent')
    parent = options.parent;
else
    parent = 0;			% including or not (1/0) in the neighborhood a coefficient from the same spatial location
end
if isfield(options, 'boundary')
    boundary = options.boundary;
else
    boundary = 1;		% Boundary mirror extension, to avoid boundary artifacts
end
if isfield(options, 'covariance')
    covariance = options.covariance;
else
    covariance = 1;     % Full covariance matrix (1) or only diagonal elements (0).
end
if isfield(options, 'optim')
    optim = options.optim;
else
    optim = 1;          % Bayes Least Squares solution (1), or MAP-Wiener solution in two steps (0)
end

PS = ones(size(x));	% power spectral density (in this case, flat, i.e., white noise)
seed = 0;

% Uncomment the following 4 code lines for reproducing the results of our IEEE Trans. on Im. Proc., Nov. 2003 paper
% This configuration is slower than the previous one, but it gives slightly better results (SNR)
% on average for the test images "lena", "barbara", and "boats" used in the
% cited article.
% Nor = 8;                           % 8 orientations
% repres1 = 'fs';                    % Full Steerable Pyramid, 5 scales for 512x512
% repres2 = '';                      % Dummy parameter when using repres1 = 'fs'   
% parent = 1;                        % Include a parent in the neighborhood


% Call the denoising function
y = denoi_BLS_GSM(x, sig, PS, blSize, parent, boundary, Nsc, Nor, covariance, optim, repres1, repres2, seed);




function im_d = denoi_BLS_GSM(im, sig, PS, blSize, parent, boundary, Nsc, Nor, covariance, optim, repres1, repres2, seed);

% [im_d,im,SNR_N,SNR,PSNR] = denoi_BLS_GSM(im, sig, ft, PS, blSize, parent, boundary, Nsc, Nor, covariance, optim, repres1, repres2, seed);
%
%	im:	input noisy image
%	sig:	standard deviation of noise
%	PS:	Power Spectral Density of noise ( fft2(autocorrelation) )
%			NOTE: scale factors do not matter. Default is white.
%	blSize: 2x1 or 1x2 vector indicating local neighborhood
%		([sY sX], default is [3 3])
%   parent: Use parent yes/no (1/0)
%	Nsc:	Number of scales
%   Nor:  Number of orientations. For separable wavelets this MUST be 3.
%   covariance: Include / Not Include covariance in the GSM model (1/0)
%   optim: BLS / MAP-Wiener(2-step) (1/0)
%   repres1: Possible choices for representation:
%           'w':    orthogonal wavelet
%                   (uses buildWpyr, reconWpyr)
%                   repres2 (optional):
%                      haar:              - Haar wavelet.
%                      qmf8, qmf12, qmf16 - Symmetric Quadrature Mirror Filters [Johnston80]
%                      daub2,daub3,daub4  - Daubechies wavelet [Daubechies88] (#coef = 2N, para daubN).
%                      qmf5, qmf9, qmf13: - Symmetric Quadrature Mirror Filters [Simoncelli88,Simoncelli90]
%           'uw':   undecimated orthogonal wavelet, Daubechies, pyramidal version
%                   (uses buildWUpyr, reconWUpyr).
%                   repres2 (optional): 'daub<N>', where <N> is a positive integer (e.g., 2)
%           's':    steerable pyramid [Simoncelli&Freeman95].
%                   (uses buildSFpyr, reconSFpyr)
%           'fs':   full steerable pyramid [Portilla&Simoncelli02].
%                   (uses buildFullSFpyr2, reconsFullSFpyr2)
%   seed (optional):    Seed used for generating the Gaussian noise (when ft == 0)
%                       By default is 0.
%
%   im_d: denoising result

% Javier Portilla, Univ. de Granada, 5/02
% revision 31/03/2003
% revision 7/01/2004
% Last revision 15/11/2004

if ~exist('blSize'),
    blSzX = 3;	% Block size
    blSzY = 3;
else
    blSzY = blSize(1);
    blSzX = blSize(2);
end

if (blSzY/2==floor(blSzY/2))|(blSzX/2==floor(blSzX/2)),
    error('Spatial dimensions of neighborhood must be odd!');
end

if ~exist('PS'),
    no_white = 0;   % Power spectral density of noise. Default is white noise
else
    no_white = 1;
end

if ~exist('parent'),
    parent = 1;
end

if ~exist('boundary'),
    boundary = 1;
end

if ~exist('Nsc'),
    Nsc = 4;
end

if ~exist('Nor'),
    Nor = 8;
end

if ~exist('covariance'),
    covariance = 1;
end

if ~exist('optim'),
    optim = 1;
end

if ~exist('repres1'),
    repres1 = 'fs';
end

if ~exist('repres2'),
    repres2 = '';
end

if (((repres1=='w') | (repres1=='uw')) & (Nor~=3)),
    warning('For X-Y separable representations Nor must be 3. Nor = 3 assumed.');
    Nor = 3;
end

if ~exist('seed'),
    seed = 0;
end

[Ny Nx] = size(im);

% We ensure that the processed image has dimensions that are integer
% multiples of 2^(Nsc+1), so it will not crash when applying the
% pyramidal representation. The idea is padding with mirror reflected
% pixels (thanks to Jesus Malo for this idea).

Npy = ceil(Ny/2^(Nsc+1))*2^(Nsc+1);
Npx = ceil(Nx/2^(Nsc+1))*2^(Nsc+1);

if Npy~=Ny | Npx~=Nx,
    Bpy = Npy-Ny;
    Bpx = Npx-Nx;
    im = bound_extension(im,Bpy,Bpx,'mirror');
    im = im(Bpy+1:end,Bpx+1:end);	% add stripes only up and right
end


% size of the extension for boundary handling
if (repres1 == 's') | (repres1 == 'fs'),
    By = (blSzY-1)*2^(Nsc-2);
    Bx = (blSzX-1)*2^(Nsc-2);
else
    By = (blSzY-1)*2^(Nsc-1);
    Bx = (blSzX-1)*2^(Nsc-1);
end

if ~no_white,       % White noise
    PS = ones(size(im));
end

% As the dimensions of the power spectral density (PS) support and that of the
% image (im) do not need to be the same, we have to adapt the first to the
% second (zero padding and/or cropping).

PS = fftshift(PS);
isoddPS_y = (size(PS,1)~=2*(floor(size(PS,1)/2)));
isoddPS_x = (size(PS,2)~=2*(floor(size(PS,2)/2)));
PS = PS(1:end-isoddPS_y, 1:end-isoddPS_x);          % ensures even dimensions for the power spectrum
PS = fftshift(PS);

[Ndy,Ndx] = size(PS);   % dimensions are even

delta = real(ifft2(sqrt(PS)));
delta = fftshift(delta);
aux = delta;
delta = zeros(Npy,Npx);
if (Ndy<=Npy)&(Ndx<=Npx),
    delta(Npy/2+1-Ndy/2:Npy/2+Ndy/2,Npx/2+1-Ndx/2:Npx/2+Ndx/2) = aux;
elseif (Ndy>Npy)&(Ndx>Npx),
    delta = aux(Ndy/2+1-Npy/2:Ndy/2+Npy/2,Ndx/2+1-Npx/2:Ndx/2+Npx/2);
elseif (Ndy<=Npy)&(Ndx>Npx),
    delta(Npy/2+1-Ndy/2:Npy/2+1+Ndy/2-1,:) = aux(:,Ndx/2+1-Npx/2:Ndx/2+Npx/2);
elseif (Ndy>Npy)&(Ndx<=Npx),
    delta(:,Npx/2+1-Ndx/2:Npx/2+1+Ndx/2-1) = aux(Ndy/2+1-Npy/2:Ndy/2+1+Npy/2-1,:);
end

if repres1 == 'w',
    PS = abs(fft2(delta)).^2;
    PS = fftshift(PS);
    % noise, to be used only with translation variant transforms (such as orthogonal wavelet)
    delta = real(ifft2(sqrt(PS).*exp(j*angle(fft2(randn(size(PS)))))));
end

%Boundary handling: it extends im and delta
if boundary,
    im = bound_extension(im,By,Bx,'mirror');
    if repres1 == 'w',
        delta = bound_extension(delta,By,Bx,'mirror');
    else
        aux = delta;
        delta = zeros(Npy + 2*By, Npx + 2*Bx);
        delta(By+1:By+Npy,Bx+1:Bx+Npx) = aux;
    end
else
    By=0;Bx=0;
end

delta = delta/sqrt(mean2(delta.^2));    % Normalize the energy (the noise variance is given by "sig")
delta = sig*delta;                      % Impose the desired variance to the noise


% main
t1 = clock;
if repres1 == 's',  % standard steerable pyramid
    im_d = decomp_reconst(im, Nsc, Nor, [blSzX blSzY], delta, parent,covariance,optim,sig);
elseif repres1 == 'fs', % full steerable pyramid
    im_d = decomp_reconst_full(im, Nsc, Nor, [blSzX blSzY], delta, parent, covariance, optim, sig);
elseif repres1 == 'w',  % orthogonal wavelet
    if ~exist('repres2'),
        repres2 = 'daub1';
    end
    filter = repres2;
    im_d = decomp_reconst_W(im, Nsc, filter, [blSzX blSzY], delta, parent, covariance, optim, sig);
elseif repres1 == 'uw',    % undecimated daubechies wavelet
    if ~exist('repres2'),
        repres2 = 'daub1';
    end
    if repres2(1:4) == 'haar',
        daub_order = 2;
    else
        daub_order = 2*str2num(repres2(5));
    end
    im_d = decomp_reconst_WU(im, Nsc, daub_order, [blSzX blSzY], delta, parent, covariance, optim, sig);
else
    error('Invalid representation parameter. See help info.');
end
t2 = clock;
elaps = t2 - t1;
elaps(4)*3600+elaps(5)*60+elaps(6); % elapsed time, in seconds

im_d = im_d(By+1:By+Npy,Bx+1:Bx+Npx);
im_d = im_d(1:Ny,1:Nx);

function fh = decomp_reconst_full(im,Nsc,Nor,block,noise,parent,covariance,optim,sig);

% Decompose image into subbands, denoise, and recompose again.
%		fh = decomp_reconst(im,Nsc,Nor,block,noise,parent,covariance,optim,sig);
%       covariance:	 are we considering covariance or just variance?
%       optim:		 for choosing between BLS-GSM (optim = 1) and MAP-GSM (optim = 0)
%       sig:        standard deviation (scalar for uniform noise or matrix for spatially varying noise)
% Version using the Full steerable pyramid (2) (High pass residual
% splitted into orientations).

% JPM, Univ. de Granada, 5/02
% Last Revision: 11/04

if (block(1)/2==floor(block(1)/2))|(block(2)/2==floor(block(2)/2)),
    error('Spatial dimensions of neighborhood must be odd!');
end

if ~exist('parent'),
    parent = 1;
end

if ~exist('covariance'),
    covariance = 1;
end

if ~exist('optim'),
    optim = 1;
end

if ~exist('sig'),
    sig = sqrt(mean(noise.^2));
end

[pyr,pind] = buildFullSFpyr2(im,Nsc,Nor-1);
[pyrN,pind] = buildFullSFpyr2(noise,Nsc,Nor-1);
pyrh = real(pyr);
Nband = size(pind,1)-1;
for nband = 2:Nband, % everything except the low-pass residual
    fprintf('%d % ',round(100*(nband-1)/(Nband-1)))
    aux = pyrBand(pyr, pind, nband);
    auxn = pyrBand(pyrN, pind, nband);
    [Nsy,Nsx] = size(aux);
    prnt = parent & (nband < Nband-Nor);   % has the subband a parent?
    BL = zeros(size(aux,1),size(aux,2),1 + prnt);
    BLn = zeros(size(aux,1),size(aux,2),1 + prnt);
    BL(:,:,1) = aux;
    BLn(:,:,1) = auxn*sqrt(((Nsy-2)*(Nsx-2))/(Nsy*Nsx));     % because we are discarding 2 coefficients on every dimension
    if prnt,
        aux = pyrBand(pyr, pind, nband+Nor);
        auxn = pyrBand(pyrN, pind, nband+Nor);
        if nband>Nor+1,     % resample 2x2 the parent if not in the high-pass oriented subbands.
            aux = real(expand(aux,2));
            auxn = real(expand(auxn,2));
        end
        BL(:,:,2) = aux;
        BLn(:,:,2) = auxn*sqrt(((Nsy-2)*(Nsx-2))/(Nsy*Nsx)); % because we are discarding 2 coefficients on every dimension
    end

    sy2 = mean2(BL(:,:,1).^2);
    sn2 = mean2(BLn(:,:,1).^2);
    if sy2>sn2,
        SNRin = 10*log10((sy2-sn2)/sn2);
    else
        disp('Signal is not detectable in noisy subband');
    end

    % main
    BL = denoi_BLS_GSM_band(BL,block,BLn,prnt,covariance,optim,sig);
    pyrh(pyrBandIndices(pind,nband)) = BL(:)';
end
fh = reconFullSFpyr2(pyrh,pind);


function fh = decomp_reconst_W(im,Nsc,filter,block,noise,parent,covariance,optim,sig);

% Decompose image into subbands, denoise using BLS-GSM method, and recompose again.
%		fh = decomp_reconst(im,Nsc,filter,block,noise,parent,covariance,optim,sig);
%       im:         image
%       Nsc:        number of scales
%       filter:     type of filter used (see namedFilters)
%       block:      2x1 vector indicating the dimensions (rows and columns) of the spatial neighborhood 
%       noise:      signal with the same autocorrelation as the noise
%       parent:     include (1) or not (0) a coefficient from the immediately coarser scale in the neighborhood
%       covariance:	 are we considering covariance or just variance?
%       optim:		 for choosing between BLS-GSM (optim = 1) and MAP-GSM (optim = 0)
%       sig:        standard deviation (scalar for uniform noise or matrix for spatially varying noise)
% Version using a critically sampled pyramid (orthogonal wavelet), as implemented in MatlabPyrTools (Eero).

% JPM, Univ. de Granada, 3/03

if (block(1)/2==floor(block(1)/2))|(block(2)/2==floor(block(2)/2)),
   error('Spatial dimensions of neighborhood must be odd!');
end   

if ~exist('parent'),
        parent = 1;
end

if ~exist('covariance'),
        covariance = 1;
end

if ~exist('optim'),
        optim = 1;
end

if ~exist('sig'),
        sig = sqrt(mean(noise.^2));
end

Nor = 3;    % number of orientations: vertical, horizontal and mixed diagonals (for compatibility)

[pyr,pind] = buildWpyr(im,Nsc,filter,'circular');
[pyrN,pind] = buildWpyr(noise,Nsc,filter,'circular');
pyrh = pyr;