denoise_main.m

MATLAB Denoising software for grayscale images

function im_d = denoise_main(im,noiseVar,opts);

% function im_d = denoise_main(im,noiseVar,options);
%
% The main program of the OE denoiser. The image 'im' which is contaminated
% with gaussian random pixel noise of variance 'noiseVar' is being
% denoised. The parameter settings are encapsulated in 'options'. The
% function 'newOptions' will generate this settings. The denoised image
% will be returned in 'im_d'. As paramters one can choose:
%
% 'opts.pyrtype'        The pyramid to use, possible values are
%       'w'         QMF-9 pyramid
%       'sf'        steerable pyramid defined in the fourier domain
%       'fullsf'    Full steerable pyramid with oriented highpass residuals
%
% 'opts.alphadef'       the default value of alpha
% 'opts.adaptalpha'      1 if alpha is estimated with moment matching
% 'opts.Nsize2'         the neighbourhood size for variance estimation
% 'opts.Nsize4'         the neighbourhood size for estimation of fourth moment
% 'opts.covariance'     use the complete noise covariance in the wavelet domain
%                       instead of just the 'noiseVar'
% 'opts.nOrientations'  only if oriented transforms are used, otherwise it
%                       is defined by the wavelet pyramid
% 'opts.adaptneighbours'Use adaptive neighbourhood selection otherwise
%                       set 'opts.block'
% 'opts.thresh'         the threshold for neighbourhood selection. See the
%                       function 'adaptiveNeighbourhood' for details
% 'opts.block'          use a block neighborhood surrounding the center
% 'opts.parent'         use the parent coefficient if existent
% 'opts.mapsteps'       how many steps of the IWF should be used?
% 'opts.daub_order'     if the WU pyramid is used, specify the length of
%                       the daubechies wavelet filter
%
% Example :
%   options = newOptions;
%   options.thresh = 0.3;
%   im_denoised = denoise_main(im_noisy,15^2,options);
%   PSNR(im,im_denoised)
%
% We thank Javier Portilla for the permission to use pieces of his code
% obtained from http://decsai.ugr.es/~javier/denoise/index.html. The code
% was used for the results in the following paper:
% 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.
%
% Peter-Vincent Gehler 15.March 2005
% last updated 31.May 2005

if nargin~=3,error('Wrong number of input arguments'); end

% set verbosity level
verbose = 0;

% Get the size of the image ...
[Ny Nx] = size(im);

% ... and compute the height of the wp ...
if ~isfield(opts,'nLevels')
    opts.nLevels = ceil(log2(min(Ny,Nx))-4);
end

if ~isfield(opts,'bigblock')
    opts.bigblock = [7 7];
end


Npy = ceil(Ny/2^(opts.nLevels+1))*2^(opts.nLevels+1);
Npx = ceil(Nx/2^(opts.nLevels+1))*2^(opts.nLevels+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   

if opts.adaptneighbours
    extndX = max(opts.bigblock(1),7);
    extndY = max(opts.bigblock(2),7);
else
    extndX = max(opts.block(1),7);
    extndY = max(opts.block(2),7);
end

% ... extend the image ...
if isfield(opts,'bndry') & opts.bndry
    if strcmp(opts.pyrtype,'s') | strcmp(opts.pyrtype,'sf') | ...
            strcmp(opts.pyrtype,'fullsf')
        By = max((extndY-1)*2^(opts.nLevels-2),2^(opts.nLevels-1));
        Bx = max((extndX-1)*2^(opts.nLevels-2),2^(opts.nLevels-1));
    else
        By = max((extndY-1)*2^(opts.nLevels-1),2^(opts.nLevels-1));
        Bx = max((extndX-1)*2^(opts.nLevels-1),2^(opts.nLevels-1));
    end
    im = bound_extension(im,By,Bx,'mirror');
else
    Bx = 0;
    By = 0;
end

% ... needed for noise covariance estimation ...
noise = sampleFreeNoise(im,sqrt(noiseVar),opts,Bx,By,Npy,Npx);

% call depending on the type of representation ... 
switch lower(opts.pyrtype)
    case 'sf' % steerable - fourier domain
        im_d = denoiseSFpyr(im,noiseVar,noise,opts);
    case 'fullsf' % as 'sf' with oriented high pass residuals
        im_d = denoiseFullSFpyr(im,noiseVar,noise,opts);
    case 'w' % QMF-9 
        im_d = denoiseWpyr(im,noiseVar,noise,opts);
    case 'wu' % QMF-9 
        im_d = denoiseWUpyr(im,noiseVar,noise,opts);
    otherwise
        error(sprintf('Pyramid type ''%s'', not supported',opts.pyrtype));
end

% ... and remove the boundary if existent
if opts.bndry
    im_d = im_d(By+1:end-By,Bx+1:end-Bx);
end