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., 4)
%           'ur':   undecimated orthogonal wavelet with redundant orientations, Daubechies,
%                    pyramidal version 
%                   (uses buildWUpyr_red, reconWUpyr_red).
%                   repres2 (optional): 'daub<N>', where <N> is a positive integer (e.g., 4)
%           '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

% JPM, Univ. de Granada, 5/02
% revision 31/03/2003
% Last revision 7/01/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 ((repres1=='ur') & (Nor~=5)),
    warning('For this representation Nor must be 5. Nor = 5 assumed.');
    Nor = 5;
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

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

% 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);      % We center, to allow zero padding and cropping
foo = zeros(size(im));
if (size(PS,1)<=size(im,1))&(size(PS,2)<=size(im,2)),
    foo(Ny/2+1-size(PS,1)/2:Ny/2+size(PS,1)/2,Nx/2+1-size(PS,2)/2:Nx/2+size(PS,2)/2) = PS;
elseif (size(PS,1)>size(im,1))&(size(PS,2)>size(im,2)),   
    foo = PS(Ndy/2+1-size(im,1)/2:Ndy/2+size(im,1)/2,Ndx/2+1-size(im,2)/2:Ndx/2+size(im,2)/2);
elseif (size(PS,1)<=size(im,1))&(size(PS,2)>size(im,2)),   
    foo(Ny/2+1-size(PS,1)/2:Ny/2+1+size(PS,1)/2-1,:) = PS(:,Ndx/2+1-size(im,2)/2:Ndx/2+size(im,2)/2);
elseif (size(PS,1)>size(im,1))&(size(PS,2)<=size(im,2)),   
    foo(:,Nx/2+1-size(PS,2)/2:Nx/2+1+size(PS,2)/2-1)=PS(Ndy/2+1-size(im,1)/2:Ndy/2+1+size(im,1)/2-1,:);
end 
PS = foo; clear foo
PS = fftshift(PS);

if repres1 == 'w',
    % 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)))))));
else
    % if not it is preferable to use a deterministic signal with the same sample autocorrelation as the noise
    delta = real(ifft2(sqrt(PS)));
    delta = fftshift(delta); % we center it, to avoid boundary problems
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

%Boundary handling: it extends im and delta
if boundary,
    im = bound_extension(im,By,Bx,'mirror');
    delta = bound_extension(delta,By,Bx,'mirror');    
    delta = delta*sqrt((Ny + 2*By)*(Nx + 2*Bx)/(Ny*Nx));	% to compensate for the change
												            % in average energy due to the image extension.
else
	By=0;Bx=0;
end

% 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);
elseif repres1 == 'cw', % FDCT wrapping   
    im_d = decomp_reconst_CURVWRAP(im, Nsc, [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);