function ResIm = Image_Denoising_Global_Denoising(Im , trueIm)
% Denoise image by denoising its patches, using a pre-determined
% dictionary.
%
% Inputs : 
% Im      : image to denoise, double, [0 255]
% trueIm  : This is used (unfairly) to optimize the parameters for the algorithm
%
% Outputs :
% ResIm   : Result Image for the best parameter

% Parameters
imLen = numel(Im);

% Create the dictionary - Haar
[HaarDict , atomNorms] = Generate_Haar_Matrix(size(Im) , 2);
atomNorms = atomNorms(:);
invAtomNorms = 1 ./ atomNorms;
% Show_Haar_Dict(HaarDict , size(Im));
nAtoms = size(HaarDict , 2); % = 1 + 3 * nLevels

% Prepare the projections of the image onto the dictionary
projs = HaarDict' * Im(:);
absProjs = abs(projs);

% Test for different values of T
Tvalues = [0 : 200];
resPSNRs = zeros(1 , length(Tvalues));
h=waitbar(0,'Testing threshold values');
for Tind = 1 : length(Tvalues)

    waitbar(Tind/length(Tvalues));
	
	% Current threshold
	T = Tvalues(Tind);
	
	% Compute the hard-thresholding result
	resProjs = projs .* (absProjs > (T .* invAtomNorms));
	
	% Compute the result signal by multiplying by the dictionary
	currResIm = HaarDict * resProjs;
	currResIm = reshape(currResIm , size(Im));
	
	% Compute error result
	resPSNRs(Tind) = Compute_Error_Stats(currResIm , trueIm);
	
end
close(h);

% Find the best value and plot the results
maxPSNRind = find(max(resPSNRs)==resPSNRs,1); 
bestT = Tvalues(maxPSNRind);
if 1
	figure;
	plot(Tvalues , resPSNRs);
	title('PSNR as a function of T');
end

% Reconstruct the image for the best value of T
resProjs = projs .* (absProjs > (bestT .* invAtomNorms));
ResIm = HaarDict * resProjs;
ResIm = reshape(ResIm , size(Im));

return;

% [n,m]=size(Haar);
% W=[0.25*ones(n*4,1); 0.5*ones(n*3,1)];




function [HaarDict , atomNorms] = Generate_Haar_Matrix(imSize , nLevels)

atomNorms = [];

% Prepare the first level bands
basicH = [1 -1];
basicL = [1 1];
HH1 = kron(basicH' , basicH); 
HaarDict = Construct_One_Haar_Band(imSize , HH1);
atomNorms = [atomNorms ones(1 , prod(imSize)) * sqrt(sum(HH1(:).^2))];

LH1 = kron(basicL' , basicH); 
HaarDict = [HaarDict Construct_One_Haar_Band(imSize , LH1)];
atomNorms = [atomNorms ones(1 , prod(imSize)) * sqrt(sum(LH1(:).^2))];

HL1 = kron(basicH' , basicL); 
HaarDict = [HaarDict Construct_One_Haar_Band(imSize , HL1)];
atomNorms = [atomNorms ones(1 , prod(imSize)) * sqrt(sum(HL1(:).^2))];

LL1 = kron(basicL' , basicL);
if nLevels == 1, 
	HaarDict = [HaarDict Construct_One_Haar_Band(imSize , LL1)]; 
	atomNorms = [atomNorms ones(1 , prod(imSize)) * sqrt(sum(LL1(:).^2))];
end;

% Create the rest of the bands from the LL of the previous level and the basic ones
currLL = LL1;
for cLevel = 2 : nLevels
	currHH = kron(HH1 , currLL); 
	HaarDict = [HaarDict Construct_One_Haar_Band(imSize , currHH)];
	atomNorms = [atomNorms ones(1 , prod(imSize)) * sqrt(sum(currHH(:).^2))];
	
	currLH = kron(LH1 , currLL); 
	HaarDict = [HaarDict Construct_One_Haar_Band(imSize , currLH)];
	atomNorms = [atomNorms ones(1 , prod(imSize)) * sqrt(sum(currLH(:).^2))];
	
	currHL = kron(HL1 , currLL); 
	HaarDict = [HaarDict Construct_One_Haar_Band(imSize , currHL)];
	atomNorms = [atomNorms ones(1 , prod(imSize)) * sqrt(sum(currHL(:).^2))];
	
	currLL = kron(LL1 , currLL);
	if cLevel == nLevels, 
		HaarDict = [HaarDict Construct_One_Haar_Band(imSize , currLL)]; 
		atomNorms = [atomNorms ones(1 , prod(imSize)) * sqrt(sum(currLL(:).^2))];
	end;
	
end


function haarBand = Construct_One_Haar_Band(imSize , haarMatrix)

% Number of pixels
nPix = prod(imSize);
[c,r] = meshgrid(1:imSize(2) , 1:imSize(1));
c = c(:)'; r = r(:)';
nOnes = numel(haarMatrix);

% Prepare all the indices we will need to insert into the sparse matrix
sparseMatRowInd = zeros(numel(haarMatrix) , nPix);
spraseMatColInd = repmat([1:nPix] , numel(haarMatrix) , 1 );
sparseMatVal = zeros(numel(haarMatrix) , nPix);


% Run on the haarMatrix
insertInd = 0;
for dc = 0 : size(haarMatrix , 2)-1
	for dr = 0 : size(haarMatrix , 1)-1
		insertInd = insertInd + 1;
		sparseMatRowInd(insertInd , :) = sub2ind(imSize , mod((r + dr) - 1 , imSize(1)) + 1  , mod((c + dc) - 1 , imSize(2)) + 1);
		sparseMatVal(insertInd , :) = haarMatrix(dr + 1 , dc + 1);
	end
end
sparseMatVal = sparseMatVal / sum(abs(haarMatrix(:)),1);

% Create the sparse matrix itself
haarBand = sparse(sparseMatRowInd(:) , spraseMatColInd(:) , sparseMatVal(:) , nPix , nPix);

function Show_Haar_Dict(D , imSize)
figure;
multSize = max(min(floor(600 ./ imSize)) , 1);
for c1 = 1 : size(D , 2)
	currAtom = reshape(full(D(:,c1)) , imSize);
	currAtom2 = imresize(currAtom , multSize , 'nearest');
	cla;
	imshow(currAtom2 , [-1 1]);
	title(sprintf('Atom %d' , c1));
	drawnow; pause(0.05);
	
	
end

