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

library(imager)
library(png)

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=""))
}

source("nt_toolbox/toolbox_wavelet_meshes/meshgrid.R")
options(repr.plot.width=3.5, repr.plot.height=3.5)

n <- 256
name <- 'nt_toolbox/data/flowers.png'
f0 <- load_image(name, n)
imageplot(f0, 'Original image')

sigma <- 0.08

f <- f0 + sigma*as.cimg(rnorm(n**2))

imageplot(clamp(f), paste('Noisy, SNR=', snr(f0,f)) )

options(repr.plot.width=4, repr.plot.height=4)

thresh_hard <- function(u,t) { return(u*(abs(u)>t)) }
alpha <- seq(-3,3, length=1000)
plot(alpha, sapply(alpha, thresh_hard, t=1), 'l')

h <- c(0, .482962913145, .836516303738, .224143868042, -.129409522551)
h <- h/norm(h)
Jmin <- 2

a <- perform_wavortho_transf(f,Jmin,+1,h)

plot_wavelet(a,Jmin)

T <- 3*sigma

aT <- thresh_hard(a,T)

plot_wavelet(aT,Jmin)

fHard <- perform_wavortho_transf(aT,Jmin,-1,h)

imageplot(clamp(fHard), paste('Hard, SNR=', snr(f0,fHard)) )

thresh_soft <- function(u,t){ return(max(1-t/abs(u), 0)*u) }

alpha <- seq(-3,3, length=1000)
plot(alpha, sapply(alpha, thresh_soft, t=1), 'l')

T <- 3/2*sigma

aT <- matrix(sapply(a, thresh_soft, t=T), c(n,n))

aT[seq(1,n,2^Jmin),seq(1,n,2^Jmin)] <- a[seq(1,n,2^Jmin),seq(1,n,2^Jmin)]

fSoft <- perform_wavortho_transf(aT,Jmin,-1,h)

imageplot(clamp(fSoft), paste('Soft, SNR=', snr(f0,fSoft)) )

options(repr.plot.width=6, repr.plot.height=4)

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

m <- 4

dY <- meshgrid_2d(1:m, 1:m)$X ; dX <- meshgrid_2d(1:m, 1:m)$Y
delta <- cbind( array(dY, c(m*m,1)), array(dX, c(m*m,1)) )

fTI <- array(0, c(n,n))
T <- 3*sigma

for (i in 1:(m*m)){
    fS <- circshift(f, delta[i,])
    a <- perform_wavortho_transf(fS,Jmin,1,h)
    aT <- array(sapply(a, thresh_hard, t=T), c(n,n))
    fS <- perform_wavortho_transf(aT,Jmin,-1,h)
    fS <- circshift(fS,-delta[i,])
    fTI <- (i-1)/i*fTI + 1/i*fS
}

imageplot(clamp(fTI), paste('TI, SNR=', snr(f0,fTI)) )