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 <- 512
f <- load_image("nt_toolbox/data/wood.jpg", n, grayscale=0)

imageplot(f)

z <- array(0, c(1,1,1,dim(f)[4]))

fe_high <- array(0, c(1,n+2,1,dim(f)[4]))
fe_high[1,2:(n+1),1,] <- f[1,,,]

fe_middle <- array(0, c(n,n+2,1,dim(f)[4]))
fe_middle[,1,,] <- f[,1,,]
fe_middle[,2:(n+1),,] <- f
fe_middle[,n+2,,] <- f[,n,,]

fe_low <- array(0, c(1,n+2,1,dim(f)[4]))
fe_low[1,2:(n+1),1,] <- f[n,,,]

fe <- array(0, c(n+2,n+2,1,dim(f)[4]))
fe[1,,,] <- fe_high
fe[2:(n+1),,,] <- fe_middle
fe[n+2,,,] <- fe_low

circshift <- function(x,v){
    x <- roll(x, v[1], axis=1)
    x <- roll(x, v[2], axis=2)
    return(x)
}

laplacian <- function(x){
    return( 4*x - (circshift(x, c(0,1)) + circshift(x, c(1,0)) + circshift(x, c(-1,0)) + circshift(x, c(0,-1))) )
}

d <- laplacian(fe)
d <- d[2:(n+1),2:(n+1),,,drop=F]

grid <- meshgrid_2d(0:(n-1),0:(n-1))
X <- grid$X ; Y <- grid$Y
U <- 4 - 2*cos(2*X*pi/n) - 2*cos(2*Y*pi/n)

U_extended <- array(0, c(n,n,1,dim(f)[4]))
for (i in 1:dim(f)[4]){
    U_extended[,,,i] <- U
}

FFT <- apply(d, c(3,4), fft)
P <- array(0, dim(f))
for (i in 1:n){ for (j in 1:n){ P[i,j,,] <- FFT[(i-1)*n + j,,] } }
P <- P / U_extended

P[1,1,1,] <- apply(f, 4, sum)

IFFT <- apply(P, c(3,4), fft, inverse=TRUE) / (n*n)
p <- array(0, dim(f))
for (i in 1:n){ for (j in 1:n){ p[i,j,,] <- Re(IFFT[(i-1)*n + j,,]) } }

mydisp <- function(x){
    x_4 <- array(0, c(2*dim(x)[1], 2*dim(x)[2], dim(x)[3:length(dim(x))]))
    x_4[1:dim(x)[1],1:dim(x)[2],,] <- x
    x_4[(dim(x)[1]+1):dim(x_4)[1],1:dim(x)[2],,] <- x
    x_4[1:dim(x)[1],(dim(x)[2]+1):dim(x_4)[2],,] <- x
    x_4[(dim(x)[1]+1):dim(x_4)[1],(dim(x)[2]+1):dim(x_4)[2],,] <- x
    return(x_4)
}

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

imageplot(as.cimg(mydisp(f)), 'Original, periodized', c(1,2,1))
imageplot(as.cimg(mydisp(p)), 'Periodic layer, periodized', c(1,2,2))

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

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

## Insert your code here.

f0 <- apply(p, c(2,1), mean)

u <- f0
u[1,1] <- 0
u[2,1] <- 1

options(repr.plot.width=3.5, repr.plot.height=3.5)
imageplot(clamp(u))

w <- array(rnorm(n*n), c(n,n)) / n
w <- w - mean(w) + 1/n**2

f <- Re( fft(fft(f0)*fft(w), inverse=TRUE)/length(f0) )

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

u <- f0
u[1,1] <- 0
u[2,1] <- 1
imageplot(clamp(u), 'Input', c(1,2,1))

u <- f
u[1,1] <- 0
u[2,1] <- 1
imageplot(clamp(u), 'Synthesized', c(1,2,2))

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

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

## Insert your code here.

f0 <- p

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

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

## Insert your code here.

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

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

## Insert your code here.

n0 <- n*2

a <- 10/n

phi <- function(t){ sin(t/a*pi/2)**2 }
psi <- function(t){ phi(t)*(t < a) + (t >= a)*(t <= 1-a) + phi(1-t)*(t > 1-a) }

t <- seq(0, 1, length=n)
g <- psi(t)*t(psi(t))

g_extended <- array(0, c(n,n,1,dim(f0)[4]))
for (k in 1:dim(f0)[4]){
   g_extended[,,,k] <- g
  }
g <- g_extended     

mu0 <- array(rep(apply(f0, 4, mean), each=n*n), c(n,n,1,3))
f1 <- array(rep(apply(f0, 4, mean), each=n0*n0), c(n0,n0,1,3))
f1[(round(n0/2) - round(n/2) + 1) : (round(n0/2) + round(n/2))
   ,(round(n0/2) - round(n/2) + 1) : (round(n0/2) + round(n/2)),,] <- mu0 + n0/n*g*(f0-mu0)

u <- f1
u[1,1,,] <- c(0,0,0)
u[2,1,,] <- c(1,1,1)

options(repr.plot.width=4.5, repr.plot.height=4.5)
imageplot(as.cimg(clamp(u)))

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

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

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