function [imden,cput]=den2d_blockTI(x,qmf,thd,L,sigma)
% Block image denoising in the TI-wavelet domain.
%
% 	Input
%		x	Input image.
%		qmf	QMF filter.
%		thd	Threshold (default, 4.5).
%		L	Block size (default, floor(sqrt(2*J))).
%		sigma	Noise std (default, estimated by MAD).
%
%	Output
%		imden	Denoised image.
%		cput	A vector of CPU time of each step: 
%			+ cpu(1) UDWT transform.
%			+ cpu(2) denoise blocks.
%			+ cpu(3) inverse UDWT transform.
%

if (exist('x')~=1), 	error('Provide input image'); 		end;
if (exist('qmf')~=1), 	qmf = MakeONFilter('Symmlet',6); 	end;
if (exist('thd')~=1), 	thd = 4.5;				end;

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

n 	= nx;
J	= nextpow2(n);
if (exist('L')~=1), 	L   = floor(sqrt(2*J)); 		end; % Size of block.
scale 	= floor(log2(L));   % Coarsest scale L#2^Jc matches the size of the block.

t = cputime;
[ll,wc] = mrdwt(x,qmf,J-scale);
cput(1) = cputime - t;

% Assuming the noise is iid Gaussian process and estimate the standard deviation as: sd # MAD(x)/0.6745
% To estimate the WGN standard deviation using the MAD.
if (exist('sigma')~=1),
 HH 	 = wc(:,2*n+1:3*n);
 sigma 	 = MAD(HH(:));
 clear HH
end

t = cputime;
% Index blocks.
[dX,dY,X,Y] = ndgrid(0:L-1,0:L-1,1:L:n-L+1,1:L:n-L+1);
I = X+dX + (Y+dY-1)*n;

for j = 0:J-1-scale
 HH = wc(:,j*3*n+2*n+1:j*3*n+3*n);
 % reshape as block
 H = HH(I);
 % block soft threshold
 Hmask = max(1-thd*log(n)*sigma^2./sum(sum(H.^2)),0);
 H =  repmat(Hmask,[L L 1 1]) .* H;
 HH(I) = H;
 wc(:,j*3*n+2*n+1:j*3*n+3*n) = HH;
 
 HL = wc(:,j*3*n+n+1:j*3*n+2*n);
 % reshape as block
 H = HL(I);
 % block soft threshold
 Hmask = max(1-thd*log(n)*sigma^2./sum(sum(H.^2)),0);
 H =  repmat(Hmask,[L L 1 1]) .* H;
 HL(I) = H;
 wc(:,j*3*n+n+1:j*3*n+2*n) = HL;
 
 LH = wc(:,j*3*n+1:j*3*n+n);
 % reshape as block
 H = LH(I);
 % block soft threshold
 Hmask = max(1-thd*log(n)*sigma^2./sum(sum(H.^2)),0);
 H =  repmat(Hmask,[L L 1 1]) .* H;
 LH(I) = H;
 wc(:,j*3*n+1:j*3*n+n) = LH;
end
cput(2) = cputime - t;

t = cputime;
imden = mirdwt(ll,wc,qmf,J-scale);
cput(3) = cputime - t;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = block_partition(x,L)

n 	= length(x);
nblocks = floor(n/L);

out = [];
for i=0:nblocks-1
  out = [out x(i*L+1:(i+1)*L,:)];
end

out = reshape(out(:),L^2,nblocks^2);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = invblock_partition(x,n,L)

nblocks = floor(n/L);

buf = reshape(x,L,L*nblocks^2);

out = [];
for i=0:nblocks-1
  out = [out;buf(:,i*L*nblocks+1:(i+1)*L*nblocks)];
end


% Old block implementation.
%for j = 0:J-1-scale
% HH = wc(:,j*3*n+2*n+1:j*3*n+3*n);
% HH = block_partition(HH,L);
% %HHmask = sqrt(sum(HH.^2)) >= thd*L*sigma;
% HHmask = max(1-thd*log(n)*sigma^2./sum(HH.^2),0);
% HH = HH .* repmat(HHmask,L*L,1);
% HH = invblock_partition(HH,n,L);
% wc(:,j*3*n+2*n+1:j*3*n+3*n) = HH;
% 
% HL = wc(:,j*3*n+n+1:j*3*n+2*n);
% HL = block_partition(HL,L);
% %HLmask = sqrt(sum(HL.^2)) >= thd*L*sigma;
% HLmask = max(1-thd*log(n)*sigma^2./sum(HL.^2),0);
% HL = HL .* repmat(HLmask,L*L,1);
% HL = invblock_partition(HL,n,L);
% wc(:,j*3*n+n+1:j*3*n+2*n) = HL;
% 
% LH = wc(:,j*3*n+1:j*3*n+n);
% LH = block_partition(LH,L);
% %LHmask = sqrt(sum(LH.^2)) >= thd*L*sigma;
% LHmask = max(1-thd*log(n)*sigma^2./sum(LH.^2),0);
% LH = LH .* repmat(LHmask,L*L,1);
% LH = invblock_partition(LH,n,L);
% wc(:,j*3*n+1:j*3*n+n) = LH;
%end

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