% clc;
clear; close all;
YES = 0;NO = 1;

Im0 = double(imread('Heart.png'));
Im0=imresize(Im0,[256 256]);
f=1.5; % scale for function truesize
[Ny,Nx] = size(Im0); % image size
Im0 = double(Im0); Im0 = Im0/ max(Im0(:)); % normalize image
k = 0.05; % noise level
ImOrig = Im0; Im0 = Im0 + k*randn(size(Im0)); % add noise
Min1 = min(Im0(:)); Max1 = max(Im0(:));
Im0 = ( Im0-Min1 )/ ( Max1-Min1 ); % noisy image
ImOrig = ( ImOrig-Min1 )/ ( Max1-Min1 ); % original image

% parameters for Weight function
m = 5; % patch size
w = 11; % window search size
h = (0.25)^2; % scale parameter for weight function
NbNeigh = 10; % number of neighbors
includeCloseNeigh = YES; % include 4 closest neighbors
includeCloseNeigh = NO;
VecParametersW = [ Ny; Nx; m; w; h; NbNeigh; includeCloseNeigh; ];

% parameters for Split-Bregman Nonlocal TV
mu = 50; % 
lambdaInv = 0.05; %  
lambda = 1/lambdaInv;
nbIters = 4; % nomber of outer iterations
NbInnerIter = 2; % number of inner iterations
VecParametersNLTV = [ Ny; Nx; m; w; NbNeigh; lambda; mu;...
    nbIters; NbInnerIter; includeCloseNeigh;];

% other parameters
cpt_fig = 0; % figure number

% % other parameters
cpt_fig = 0; % figure number

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute NL-Weights
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%file = 'weights.mat';
[W,Y,SY] = compute_fastNLWeights_mex(single(Im0),...
    single(VecParametersW));
%save(file,'W','Y','SY');
%load(file,'W','Y','SY');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Perform NL-TV
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialization
NbNeigh2 = VecParametersNLTV(5);
includeCloseNeigh2 = VecParametersNLTV(10);
if includeCloseNeigh2==YES; i=4; else; i=0; end
d = zeros(Nx,Ny,2*(NbNeigh2+i));
b = zeros(Nx,Ny,2*(NbNeigh2+i));
u = Im0;

[u_new,d_new,b_new,Temp] = SBNLTV_mex(single(u),single(d),single(b),...
    single(Im0),single(W),int32(Y),int32(SY),single(VecParametersNLTV));

% [u_new,d_new,b_new,Temp] = SBNLTV_mex(single(u),single(d),single(b),...
%     single(Im0),single(W),int32(Y),int32(SY),single(VecParametersNLTV));

u_new=double(u_new);
Lh=1;
wav = daubcqf(2,'min');
xwt1=mdwt(u,wav,Lh);
xwt2=mdwt(u_new,wav,Lh);

im2=midwt(xwt2,wav,Lh);

figure;subplot(2,3,1);imshow(u,[]);
subplot(2,3,2);imshow(u_new,[]);
subplot(2,3,3);imshow(im2,[]);
subplot(2,3,4);imshow(xwt1,[]);
subplot(2,3,5);imshow(xwt2,[]);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute SNR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
var_u = var(u_new(:));
var_u_uRef = var(u_new(:)-ImOrig(:));
SNR = 10*log10(var_u/var_u_uRef)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Display results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cpt_fig = cpt_fig + 1;
figure(cpt_fig); clf;
subplot(1,2,1); imagesc(Im0,[0 1]); colormap(gray); title('Im0'); hold on; %colorbar;
subplot(1,2,2); imagesc(u_new,[0 1]); colormap(gray); title(['u, SNR= ',num2str(SNR)]); hold on; %colorbar;
truesize(cpt_fig,[round(f*Ny) round(f*Nx)]);