close all;
clear all;
imagedata = double(imread('barbara.pgm'));
% imagedata = rotplane(32,32, 45);
[Faf, Fsf] = Lift53;
Faf{1}(:,1) = Faf{1}(:,1)/sqrt(2); 
Faf{1}(:,2) = Faf{1}(:,2)*sqrt(2);
Faf{2}(:,1) = Faf{2}(:,1)/sqrt(2);
Faf{2}(:,2) = Faf{2}(:,2)*sqrt(2);


Fsf{1}(:,1) = Fsf{1}(:,1)*sqrt(2); 
Fsf{1}(:,2) = Fsf{1}(:,2)/sqrt(2);
Fsf{2}(:,1) = Fsf{2}(:,1)*sqrt(2);
Fsf{2}(:,2) = Fsf{2}(:,2)/sqrt(2);

af = Faf;
sf = Fsf;

[Fafk, Fsfk] = FSfarras;
Fafk{1}(:,1) = Fafk{1}(:,1)/sqrt(2); 
Fafk{1}(:,2) = Fafk{1}(:,2)*sqrt(2);
Fafk{2}(:,1) = Fafk{2}(:,1)/sqrt(2);
Fafk{2}(:,2) = Fafk{2}(:,2)*sqrt(2);


Fsfk{1}(:,1) = Fsfk{1}(:,1)*sqrt(2); 
Fsfk{1}(:,2) = Fsfk{1}(:,2)/sqrt(2);
Fsfk{2}(:,1) = Fsfk{2}(:,1)*sqrt(2);
Fsfk{2}(:,2) = Fsfk{2}(:,2)/sqrt(2);

[afk sfk]= dualfilt1;
afk{1}(:,1) = afk{1}(:,1)/sqrt(2); 
afk{1}(:,2) = afk{1}(:,2)*sqrt(2);
afk{2}(:,1) = afk{2}(:,1)/sqrt(2);
afk{2}(:,2) = afk{2}(:,2)*sqrt(2);

sfk{1}(:,1) = sfk{1}(:,1)*sqrt(2); 
sfk{1}(:,2) = sfk{1}(:,2)/sqrt(2);
sfk{2}(:,1) = sfk{2}(:,1)*sqrt(2);
sfk{2}(:,2) = sfk{2}(:,2)/sqrt(2);

J = 5;
T1 = 128;
T2 = 1;
Step = -1;
T = T1:Step:T2;
k=1.6;
TT = [50 40 30 20 10 5];

wdt = dualtree2D(imagedata, J, Faf, af);
wdtns = dtns2D(imagedata, wdt, J, Faf, af, Fsf, sf, T, k);

wdtk = dualtree2D(imagedata, J, Fafk, afk);
wdtnsk = dtns2D(imagedata, wdtk, J, Fafk, afk, Fsfk, sfk, T, k);

% Configuration 1
% [coeff, realpart, imagepart] = dtcell2arr(wdtns, J)
% Configuration 2
% step1: my = dtmerge(wdtns, J);
% step2: allceoff = dwtcell2arr(my, J);
my = dtmerge(wdtns, J);
allceoff = dwtcell2arr(my, J);

myk = dtmerge(wdtnsk, J);
allceoffk = dwtcell2arr(myk, J);

% fid = fopen('e:\\coeff','w');
% fwrite(fid, allceoff', 'int32');
% fclose(fid);
% 
% fid = fopen('e:\\coeff_r','r');
% coeff_r = fread(fid,[1024,512],'int32');
% coeff_r=coeff_r';
% fclose(fid);

wdt_r = idtmerge(allceoff,J);
wdtk_r = idtmerge(allceoffk,J);

% [coeff_r, realpart, imagpart] = dtcell2arr(wdt_r, J)
% coeff = [realpart, imagepart];
% y = dtarr2cell(coeff_r, J);

for i = 1:length(TT)

w1 = DTThrshld(wdtns,J,TT(i));
w2 = DTThrshld(wdtnsk,J,TT(i));

Ry1{i} = idualtree2D(w1, J, Fsf, sf);
Ry2{i} = idualtree2D(w2, J, Fsfk, sfk);

PSNR1(i) = calcpsnr(Ry1{i}, imagedata);
PSNR2(i) = calcpsnr(Ry2{i}, imagedata);

end

plot(TT,PSNR1,'-bd',TT,PSNR2,'-rs');
legend('DDWT 5/3','DDWT kingsbury');
xlabel('Threshold');
ylabel('PSNR');

% recimage = idualtree2D(wdt_r, J, Fsf, sf);
% 
% psnr = calcpsnr(recimage, imagedata)
% max = max(max((abs(recimage-imagedata))))
% min = min(min((abs(recimage-imagedata))))
% figure;
% imshow(recimage, [])