using PyPlot
using NtToolBox
using Distributions

n = 256
N = n^2;

name = "NtToolBox/src/data/flowers.png"
x0 = load_image(name, n);

imageplot(x0)

sigma = .08

y = x0 .+ sigma.*rand(Normal(), 256, 256);

imageplot(clamP(y))

cconv = (a, b) -> real(plan_ifft((plan_fft(a)*a).*(plan_fft(b)*b))*((plan_fft(a)*a).*(plan_fft(b)*b)))

normalize = h -> h/sum(h)

X = [0:n/2; -n/2:-2]'
Y = [0:n/2; -n/2:-2]

h = mu -> normalize(exp(-(X.^2 .+ Y.^2)/(2*(mu)^2)))

mu = 10
subplot(1,2, 1)
imageplot(fftshift(h(mu)))
title("h")
subplot(1,2, 2)
imageplot(fftshift(real(plan_fft(h(mu))*h(mu))))
title(L"$\hat h$")

imageplot(h(mu))

denoise = (x, mu) -> cconv(h(mu), x)

imageplot(denoise(y, mu))

include("NtSolutions/denoisingsimp_2b_linear_image/exo1.jl")

include("NtSolutions/denoisingsimp_2b_linear_image/exo2.jl")
# mulist = linspace(.1, 3.5, 31)
# err = zeros(length(mulist))
# for i in 1 : length(mulist)
# 	mu = mulist[i]
# 	err[i] = norm(x0 - denoise(y, mu))
# end

# clf
# h1, = plot(mulist,err); axis("tight");
# # retrieve the best denoising result
# i = mapslices(indmin, err, 1)[1]
# mu = mulist[i]

imageplot(denoise(y, mu))

P = 1/N .* ( abs(plan_fft(x0)*x0).^2 );

h_w = real(plan_ifft(P ./ (P .+ sigma^2))*(P ./ (P .+ sigma^2)));

u = fftshift(h_w)
imageplot( u[Int(n/2 - 10) : Int(n/2 + 10), Int(n/2 - 10) : Int(n/2 + 10)] )

imageplot( cconv(y, h_w) )
<script>
  $(document).ready(function(){
      $('div.prompt').hide();
  });
</script>