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

% Decompose image into subbands (undecimated wavelet), denoise, and recompose again.
%		fh = decomp_reconst_CURVWRAP(im,Nsc,block,noise,parent,covariance,optim,sig);
%       im :         image
%       Nsc:         Number of scales
%       block:       size of neighborhood within each undecimated subband.
%       noise:       image having the same autocorrelation as the noise (e.g., a delta, for white noise)
%       parent:      are we including the coefficient at the central location at the next coarser scale?
%       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)

% MJF, GREYC CNRS UMR 6072, 3/06

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=4;
bands={'NORTH','EAST','SOUTH','WEST'};
[cimg,cpyr]=CurvWrapPyramid(im,Nsc+1);
[cimgn,cpyrn]=CurvWrapPyramid(noise,Nsc+1);

% Apply denoiser
cpyrh = cpyr;
for s = 1:Nsc
 for nor=1:Nor
    com=sprintf('aux=cpyr{%d}.%s;',s,bands{nor});
    eval(com);
    com=sprintf('auxn=cpyrn{%d}.%s;',s,bands{nor});
    eval(com);
    prnt = parent & (s > 2);   % 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;
    if prnt,
	[m,n] = size(aux);
	com=sprintf('aux=expandft(cpyr{%d}.%s,max(size(cpyr{%d}.%s)./size(cpyr{%d}.%s)));',...
	    	     s-1,bands{nor},s,bands{nor},s-1,bands{nor});
	eval(com);
	com=sprintf('auxn=expandft(cpyrn{%d}.%s,max(size(cpyrn{%d}.%s)./size(cpyrn{%d}.%s)));',...
	    	     s-1,bands{nor},s,bands{nor},s-1,bands{nor});
	eval(com);   
	BL(:,:,2) = aux;
	BLn(:,:,2) = auxn;
    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_mbkf_band(BL,block,BLn,prnt,covariance,optim,sig);
    com=sprintf('cpyrh{%d}.%s=BL;',s,bands{nor});
    eval(com);
 end 
end

% Take inverse curvelet transform 
fh = InvCurvWrapPyramid(cpyrh);