using NtToolBox
using PyPlot
using Autoreload
arequire("NtToolBox")

p = 4;

h = [0, .482962913145, .836516303738, .224143868042, -.129409522551];
h = h/norm(h);
g = cat(1, 0, h[length(h):-1:2]) .* ( (-1).^(1:length(h)) );

println("h filter = ", h);
println("g filter = ", g);

N = 256;
f = rand(N, 1);

a = subsampling( cconvol(f,h) );
d = subsampling( cconvol(f,g) );

e0 = norm(f[:])^2;
e1 = norm(a[:])^2 + norm(d[:])^2;
println("Energy before : $e0");
println("Energy before : $e1");

f1 =  cconvol(upsampling(a), reverse(h)) + cconvol(upsampling(d), reverse(g));

println("Error |f-f1|/|f| = ", norm(f[:]-f1[:])/norm(f[:]) );

n = 256;
name = "NtToolBox/src/data/hibiscus.png";
f = load_image(name, N);
f = rescale(sum(f,3));
f = f[:, :, 1];
imageplot(f);

j = log2(n) - 1;

fW = copy(f);

A = fW[1 : Int(2^(j + 1)), 1 : Int(2^(j + 1))];

Coarse = subsampling(cconvol(A, h, 1), 1);
Detail = subsampling(cconvol(A, g, 1), 1);

norm(A[:])^2 - norm(Coarse[:])^2 - norm(Detail[:])^2

A = cat(1, Coarse, Detail );

clf;
imageplot(f,"Original image", [1, 2, 1]);
imageplot(A,"Vertical transform", [1, 2, 2]);

Coarse = subsampling(cconvol(A, h, 2), 2);
Detail = subsampling(cconvol(A, g, 2), 2);

A = cat(2, Coarse, Detail );

fW[1 : Int(2^(j + 1)), 1 : Int(2^(j+1))] = A;

clf;
imageplot(f,"Original image", [1, 2, 1]);
subplot(1,2,2);
plot_wavelet(fW, Int(log2(n) - 1)); 
title("Transformed");

print( norm(f[:]) - norm(fW[:]) )

include("NtSolutions/wavelet_4_daubechies2d/exo1.jl");

e0 = norm(f[:])^2; 
e1 = norm(fW[:])^2; 
println("Energy of the signal       = $e0");
println("Energy of the coefficients = $e1");

clf;
imageplot(f, "Original", [1, 2, 1]);
subplot(1, 2, 2);
plot_wavelet(fW, Jmin); 
title("Transformed");

f1 = copy(fW);
j = 0;

A = f1[1 : Int(2^(j + 1)), 1 : Int(2^(j + 1))];

Coarse = A[1 : Int(2^j), :];
Detail = A[Int(2^j + 1) : Int(2^(j + 1)), :];

Coarse = cconvol(upsampling(Coarse, 1), reverse(h), 1);
Detail = cconvol(upsampling(Detail, 1), reverse(g), 1);

A = Coarse + Detail;

Coarse = A[:, 1 : Int(2^j)];
Detail = A[:, Int(2^j + 1) : Int(2^(j + 1))];

Coarse = cconvol(upsampling(Coarse, 2), reverse(h), 2);
Detail = cconvol(upsampling(Detail, 2), reverse(g), 2);

A = Coarse + Detail;

f1[1 : Int(2^(j + 1)), 1 : Int(2^(j + 1))] = A;

include("NtSolutions/wavelet_4_daubechies2d/exo2.jl");

e = norm(f[:] - f1[:])/norm(f[:]);
print("Error |f-f1|/|f| = $e");

eta = 4;
fWLin = zeros(n, n);
fWLin[1 : Int(n/eta), 1 : Int(n/eta)] = fW[1 : Int(n/eta), 1 : Int(n/eta)];

include("NtSolutions/wavelet_4_daubechies2d/exo3.jl");

T = .2;

fWT = fW .* float(abs(fW) .> T);

clf;
subplot(1,2,1);
plot_wavelet(fW,Jmin); 
title("Original coefficients");
subplot(1,2,2);
plot_wavelet(fWT,Jmin); 
title("Thresholded coefficients");

include("NtSolutions/wavelet_4_daubechies2d/exo4.jl");