function imden=den2d(x,qmf,scale,type,sigma,thd)
% Term-by-term image denoising in the orthogonal wavelet domain.
%
% 	Input
%		Image	Input image.
%		qmf	QMF filter.
%		scale	Coarsest decomposition scale.
%		type	Hard or soft thresholding ('H' (default) 'S').
%		Sigma	Noise std (default, estimated by MAD).
%		thd	Threshold (default, universal threshold).
%
%	Output
%		imden	Denoised image.
% 
if (exist('x','var')~=1), error('Provide input image'); end;
if (exist('scale','var')~=1), scale=0; end;
if (exist('qmf','var')~=1), qmf=MakeONFilter('Haar'); end;
if (exist('type','var')~=1), type='H'; end;
if (~strcmp(type,'S') & ~strcmp(type,'H')), error('Only hard or soft thresholding handled.'); end;

[nx ny] = size(x);
if nx~=ny
	disp('Matrix must be rectangular');
	return;
end

[n,J] = dyadlength(x);

wc=FWT2_PO(x,scale,qmf);

if ~exist('sigma','var')
 % Assuming the noise is iid Gaussian process and estimate the standard deviation as: sigma # MAD(x)/0.6745
 ntx=nx/2;nty=ny/2;
 hg=reshape(wc(ntx+1:nx,nty+1:ny),ntx*nty,1);
 sigma=MAD(hg);
 clear hg;
end

if ~exist('thd','var')
 thd=sigma*sqrt(2*log(nx*ny-2^scale*2^scale));
end

Coarse  = wc(1:2^scale,1:2^scale);
if(strcmp(type,'S')) Thr_wc_Noisy = SoftThresh(wc,thd);
else Thr_wc_Noisy = HardThresh(wc,thd); end;
Thr_wc_Noisy(1:2^scale,1:2^scale) = Coarse;
imden = IWT2_PO(Thr_wc_Noisy,scale,qmf);





%
% Part of DBlockToolbox Version:100
% Written by: Jalal Fadili, GREYC CNRS ENSICAEN
%             Christophe Chesneau, LMNO University of Caen
%             Jean-Luc Starck, CEA-CNRS
% Created June 2008
% This material is distributed under CeCiLL licence.
% E-mail: Jalal.Fadili@greyc.ensicaen.fr
%