samplefreenoise.m

MATLAB Denoising software for grayscale images

function delta = sampleFreeNoise(im, sig, opts, Bx, By, Npy,Npx);

% delta = sampleFreeNoise(im,sig,opts,Bx,By,Npy,Npx)
%
% function to compute an image which same covariance as a a noisy image
% with variance 'sig'. This code is taken from the BLS_GSM toolbox. We thank 
% Javier Portilla for the permission to use his code obtained from 
% http://decsai.ugr.es/~javier/denoise/index.html. It
% 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.

blSize = opts.block;
parent = opts.parent;
Nsc = opts.nLevels;
Nor = opts.nOrientations;
repres1 = opts.pyrtype;
boundary = opts.bndry;

blSzY = blSize(1);
blSzX = blSize(2);

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

if ((strcmp(repres1,'w') | strcmp(repres1,'uw')) & (Nor~=3)),
    Nor = 3;
end    

if (strcmp(repres1,'ur') & (Nor~=5)),
    warning('For this representation Nor must be 5. Nor = 5 assumed.');
    Nor = 5;
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).

PS = ones(size(im));

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 assumed 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    

if isfield(opts,'bndry') & opts.bndry
    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