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

eta = 10

f = x -> ( x[1,1]^2 + eta*x[2, 1]^2 ) / 2

include("ndgrid.jl")
t = linspace(-.7,.7,101)
(u, v) = meshgrid(t,t)
F = ( u .^ 2 + eta .* v .^ 2 ) ./ 2

contourf(t, t, F, 35) # To check

GradF = x -> [[x[1, 1]],[eta.*x[2, 1]]]

tau = 1.8/eta

## Insert your code here.

#contourf(t,t,Jmesh,35)
#plot(X[0,:], X[1,:], 'k.-')

## Insert your code here.

n = 256
name = "NtToolBox/src/data/lena.png"
x0 = load_image(name, n)

imageplot(x0)

grad = x -> cat(3, x - [x[end, :]'; x[1:end-1, :]], x - [x[:, end] x[:,1:end-1]])


v = grad(x0)

imageplot(v[:,:,1], L"\frac{d}{dx}", (1,2,1))
imageplot(v[:,:,2], L"\frac{d}{dy}", (1,2,2))

imageplot(sqrt(sum(v .* v, 3))[:, :])

diiv = v -> [v[2:end, :, 1]; v[1, :, 1]'] - v[:, :, 1] + [v[:, 2:end, 2] v[:, 1, 2]] - v[:, :, 2] # Pour ne pas avoir conflit avec la fonction div de Julia.

delta = x -> diiv(grad(x))

imageplot(delta(x0))

print("Should be 0:\n", sum(grad(x0) .* grad(x0)) + sum(delta(x0) .* x0) )

sigma = 0.1
y = x0 + sigma .* randn(n,n)

imageplot(clamP(y))

la = 0.3 / 5

epsilon = 1e-3

norm_eps = u -> sqrt(eps^2 + sum(u .* u, 3))
J = x -> sum(sum(norm_eps(grad(x))))

f = x -> 1/2 ^ norm(x - y)^2 + la .* J(x)

normalize_vec = (u, epsilon) -> u ./ repeat(reshape(NormEps(u,epsilon), (1, 1, 2)), inner = [size(NormEps(u,epsilon))[1], size(NormEps(u,epsilon))[2], 1])
gradTV = x -> - divergence(normalize_vec(grad(x)))

Gradf = x -> x - y + la .* gradTV(x)

tau = 1.8/( 1 + la*8/epsilon )
tau = tau*4

## Insert your code here.

#imageplot(clamp(x))

n = 64
radius = 0.6
t = linspace(-1,1,n)
(Y,X) = meshgrid(t,t)
x0 = (max( abs(X), abs(Y) ) .< radius) .* 1.0

a = 4
Lambda = ones(n, n)
Lambda[Int(n/2) - a : Int(n/2) + a, :] = 0

Phi  = x -> x .* Lambda
PhiS = Phi

y = Phi(x0)

imageplot(x0, "Original", (1,2,1))
imageplot(y, "Damaged", (1,2,2))

ProjH = x -> x + PhiS( y - Phi(x) )

## Insert your code here.

## Insert your code here.