options(warn=-1) # turns off warnings, to turn on: "options(warn=0)"
library(pracma)
library(imager)

# Importing the libraries
for (f in list.files(path="nt_toolbox/toolbox_general/", pattern="*.R")) {
    source(paste("nt_toolbox/toolbox_general/", f, sep=""))
}
for (f in list.files(path="nt_toolbox/toolbox_signal/", pattern="*.R")) {
    source(paste("nt_toolbox/toolbox_signal/", f, sep=""))
}

n = 256
name = 'nt_toolbox/data/hibiscus.png'
f0 = load_image(name, n)

options(repr.plot.width=4, repr.plot.height=4)
imageplot(f0, 'Image f0')

s = 3

x = c(c(0:((n/2) - 1)), c(-(n/2):-1))
mesh = meshgrid(x)
X = mesh$X
Y = mesh$Y
h = exp( -(X**2 + Y**2)/(2.0 * s**2) )
h = h/sum(h)

hF = Re(fft(h))

imageplot(fftshift(h), 'Filter', c(1,2,1))
imageplot(fftshift(hF), 'Fourier transform', c(1,2,2))

Phi = function(x,h){Re(fft(fft(x) * fft(h), inverse=TRUE) / length(h))}

y0 = Phi(f0[,], h)

imageplot(y0, 'Observation without noise')

sigma = .02

set.seed(123)
y = y0 + randn(n,n) * sigma

imageplot(y, 'Observation with noise')

# Normalize to 0-1
imageplot(y, paste('Observation, SNR = ',round(snr(f0, y),2),'dB'))

yF = fft(y)

Lambda = 0.02

fL2 = Re(fft(yF * hF / (abs(hF)**2 + Lambda), inverse=TRUE)) / length(yF * hF / (abs(hF)**2 + Lambda))

imageplot(fL2, paste('L2 deconvolution, SNR = ', round(snr(f0, fL2),2), 'dB'))

source("nt_solutions/inverse_2_deconvolution_variation/exo1.R")

# Insert your code here.

# Insert your code here.

S = (X**2 + Y**2) * (2/n)**2

Lambda = 0.2

fSob = Re(fft(yF * hF / (abs(hF)**2 + Lambda*S), inverse=TRUE)) / length(yF * hF / (abs(hF)**2 + Lambda*S))

imageplot(fSob, paste('Sobolev deconvolution, SNR = ', round(snr(f0, fSob),2), 'dB'))

source("nt_solutions/inverse_2_deconvolution_variation/exo2.R")

# Insert your code here.

# Insert your code here.

epsilon = 0.4*1e-2

Lambda = 0.06

tau = 1.9 / (1 + Lambda * 8 / epsilon)

fTV = y

niter = 300

#importing grad and div
grad = grad_3
div = div_2

a = randn(n,n)
b = array(0, dim=c(n,n,2))
b[,,1] = randn(n, n)
b[,,2] = randn(n, n)
dotp = function(x,y){sum(x * y)}
print(paste("Should be 0: ", dotp(grad(a),b) + dotp(a,div(b))))

Gr = grad(fTV)
d = sqrt(epsilon**2 + (Gr**2)[,,1] + (Gr**2)[,,2])
G = -div(Gr / array(rep(d, 2), dim=c(dim(d),2)) )

tv = sum(d)
print(paste('Smoothed TV norm of the image: ', round(tv,2)))

e = Phi(fTV, h) - y
fTV = fTV - tau * (Phi(e, h) + Lambda * G)

source("nt_solutions/inverse_2_deconvolution_variation/exo3.R")

# Insert your code here.

# Insert your code here.

source("nt_solutions/inverse_2_deconvolution_variation/exo4.R")

# Insert your code here.

# Insert your code here.