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*2

f0 <- as.matrix(load_image("nt_toolbox/data/hibiscus.png",n))

imageplot(f0)

cconv <- function(f, h){
    c <- fft(f)*fft(h)
   return( Re(fft(c, inverse=T )/length(c)) ) }

t <- c(0:round(n/2), (round(-n/2)+1):(-1))
X2 <- meshgrid_2d(t, t)$X ; X1 <- meshgrid_2d(t, t)$Y
normalize <- function(h){ h/sum(h) }
h <- function(sigma){ normalize(exp(-(X1**2 + X2**2)/ (2*sigma**2))) }

blur <- function(f, sigma){ cconv(f, h(sigma)) }

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

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

## Insert your code here.

s <- c(n, 1:(n-1))

nabla <- function(f){
    dx <- f - f[s,]
    dy <- f - f[,s]
    grad <- array(rep(0, length=2*n*n), c(n,n,2))
    grad[,,1] <- dx ; grad[,,2] <- dy
    return(grad)
}

v <- nabla(as.matrix(f0))

options(repr.plot.width=7, repr.plot.height=3.5)

imageplot(v[,,1], "d/dx", c(1,2,1))
imageplot(v[,,2], "d/dy", c(1,2,2))

sigma <- 1
d <- sqrt(apply(nabla(blur(as.matrix(f0), sigma))**2, c(1,2), sum))

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

imageplot(d)

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

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

## Insert your code here.

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

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

## Insert your code here.

div <- function(v){
    v1 <- v[,,1]
    v2 <- v[,,2]
    t <- c(2:n,1)
    return( v1[t,] - v1 + v2[,t] - v2 )
}

delta <- function(f){ div(nabla(f)) }

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

imageplot(delta(as.matrix(f0)))

dotp <- function(a,b){ sum(a*b) }
print(paste("Should be 0:", (dotp(nabla(f0), nabla(f0)) + dotp(delta(f0), f0))))

sigma = 4

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

plot_levelset(delta(blur(f0, sigma)), f0)

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

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

## Insert your code here.

dx <- function(f){ (f[s,] - f[t,])/2 }
dy <- function(f){t(dx(t(f))) }

s <- c(2:n, 1)
t <- c(n, 1:(n-1))

d2x <- function(f){ f[s,] + f[t,] -2*f }
d2y <- function(f){ t(d2x(t(f))) }
dxy <- function(f){ dy(dx(f)) }

hessian <- function(f){
    d2x_f <- d2x(f)
    dxy_f <- dxy(f)
    d2y_f <- d2y(f)
    hess <- array(rep(0, length=3*n*n), c(n,n,3))
    hess[,,1] <- d2x_f ; hess[,,2] <- dxy_f ; hess[,,3] <- d2y_f
    return(hess) }

sigma <- 6
g <- grad(blur(f0, sigma))

H <- hessian(blur(f0, sigma))

a <- H[,,c(1,2)] * array(rep(g[,,1],2), c(n,n,2)) + H[,,c(2,3)] * array(rep(g[,,2],2), c(n,n,2))

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

plot_levelset(apply(a*g, c(1,2), sum), f0)

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

source("nt_solutions/segmentation_1_edge_detection/exo5.R")

## Insert your code here.
<script>
  $(document).ready(function(){
      $('div.prompt').hide();
  });
</script>