function imden=den2d_curvwrap(Image,scale,type,sigma,thresh)

if (exist('Image','var')~=1), error('Provide input image'); end;

[M,N]=size(Image);

if (exist('scale','var')~=1), scale=floor(log2(min(M,N)))-3; end;
if (exist('type','var')~=1), type='H'; end;
if (~strcmp(type,'S') & ~strcmp(type,'H')), error('Only hard or soft thresholding handled.'); end;
J=floor(log2(min(M,N)));
if (exist('sigma','var')~=1),
	% To estimate the WGN standard deviation using the MAD.
	% The qmf is quite arbitrary, sufficiently regular to approach a good band-pass.
	qmf=MakeONFilter('Daubechies',4);
	wc = FWT2_PO(Image,J-1,qmf);
	hh = wc(N/2+1:N,N/2+1:N);
	sigma = MAD(hh(:));
end

% Take curvelet transform
C = fdct_wrapping(Image,1,scale);

if ~exist('thresh','var')
 thresh=3*ones(length(C),1)*sigma;
 thresh(end)=4*sigma;
elseif length(thresh)==1
  thresh=thresh*ones(length(C),1)*sigma;
end

% Compute L2 norms of curvelets (for noise std computations)
E = computeL2norm(C);

% Apply thresholding
Ct = C;
for s = 2:length(C)
  for w = 1:length(C{s})
    if type=='H'
  	  Ct{s}{w} = C{s}{w}.* (abs(C{s}{w}) > thresh(s)*E{s}{w});
    else
  	  Ct{s}{w} = sign(C{s}{w}) .* max(abs(C{s}{w}) - thresh(s)*E{s}{w},0);
    end
  end
end

% Take inverse curvelet transform 
imden = real(ifdct_wrapping(Ct,1,scale));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function E=computeL2norm(coeffs)
% Compute norm of curvelets (exact)

F = ones(coeffs{1}{2});
X = fftshift(ifft2(F)) * sqrt(prod(size(F))); 	% Appropriately normalized Dirac
C = fdct_wrapping(X,1,length(coeffs)); 		% Get the curvelets
nb=prod(size(C{1}{1}));

E = cell(size(C));
for j=1:length(C)
  E{j} = cell(size(C{j}));
  for w=1:length(C{j})
    A = C{j}{w};
    E{j}{w} = sqrt(sum(sum(A.*conj(A))) / prod(size(A)));
    nb=nb+prod(size(C{j}{w}));
  end
end
%nb/prod(size(X))