function [imden,cput]=den2d_blockcurvwrap(x,scale,thd,L,sigma)
% Block image denoising in the curvelet 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) FDCT transform.
%			+ cpu(2) l_2 norms computation.
%			+ cpu(3) denoise blocks.
%			+ cpu(4) inverse FDCT transform.
%

if (exist('x')~=1), 	error('Provide input image'); 			  end;
if (exist('scale')~=1), scale = floor(2*log2(sqrt(2*nextpow2(length(x))))); end; % Coarsest scale L#2^Jc matches the size of the block.
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.

% 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),
 qmf 	= MakeONFilter('Symmlet',6);
 wc	= FWT2_PO(x,1,qmf);
 HH	= wc(n/2+1:n,n/2+1:n);
 sigma	= MAD(HH(:));
 clear HH
end

% main
% Take curvelet transform
t = cputime;
C = fdct_wrapping(x,1,J-scale);
cput(1) = cputime - t;

% Compute L2 norms of curvelets (FDCT implements a PTF).
% This is necessary for thresholding with non-normalized dictionaries.
t 	= cputime;
E 	= computeL2norm(C);
cput(2) = cputime - t;

Cden	= C;
t 	= cputime;
for sc=2:length(C)
 Nor = length(C{sc});
 for nor=1:Nor
  [nx,ny] = size(C{sc}{nor});
  nbound  = L*ceil([nx ny]/L);
  subband = zeros(nbound);
  subband(1:nx,1:ny) = C{sc}{nor};
  [dX,dY,X,Y] = ndgrid(0:L-1,0:L-1,1:L:nbound(1)-L+1,1:L:nbound(2)-L+1);
  I = X+dX + (Y+dY-1)*nbound(1);
  coeffs = subband(I);
  norml2 = sum(sum(coeffs.^2));
  cmask  = zeros(size(norml2));
  cmask(find(norml2))  = max(1-thd*log(n)*(sigma*E{sc}{nor})^2./norml2(find(norml2)),0);
  coeffs = coeffs .* repmat(cmask,[L L 1 1]);
  subband(I) = coeffs;
  Cden{sc}{nor} = subband(1:nx,1:ny);
 end
end
cput(3) = cputime - t;

t 	= cputime;
imden = real(ifdct_wrapping(Cden,1,J-scale));
cput(4) = cputime - t;

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

n	= size(x);
nbound 	= L*ceil(size(x)/L);
nblocks = floor(nbound/L);
buf	= zeros(nbound);
buf(1:n(1),1:n(2)) = x;	

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

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


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

nblocks = floor(nbound/L);

buf = reshape(x,L,L*prod(nblocks));

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

out = out(1:n(1),1:n(2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [E,Cdelta]=computeL2norm(C)
% Compute norm of curvelets (exact)

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function E=computeL2normfast(n,scale)
% Compute norm of curvelets (fast)

[c,m0] = fdct_wrapping_range(n,nextpow2(n)-scale+1);
L0 = m0/(n*n);
E  = 1/sqrt(L0);

% Old block implementation.
% for sc=2:length(C)
%  Nor = length(C{sc});
%  for nor=1:Nor
%   coeffs = C{sc}{nor};
%   coeffs = block_partition(coeffs,L);
%   norml2 = sum(coeffs.^2);
%   cmask  = zeros(size(norml2));
%   cmask(find(norml2))  = max(1-thd*log(n)*(sigma*E{sc}{nor})^2./norml2(find(norml2)),0);
%   coeffs = coeffs .* repmat(cmask,L*L,1);
%   coeffs = invblock_partition(coeffs,L*ceil(size(C{sc}{nor})/L),size(C{sc}{nor}),L);
%   Cden{sc}{nor} = coeffs;
%  end
% 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
%