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=""))
}
options(repr.plot.width=3.5, repr.plot.height=3.5)

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

u <- (-rep(1,length(h)-1))**(seq(1,length(h)-1)) # alternate +1/-1
g <- array(c(0,rev(h[2:length(h)])*u))

print(h)
print(g)

N <- 256
f <- matrix(runif(N),N)

a <- subsampling( cconv(f,h,1),1 )
d <- subsampling( cconv(f,g,1),1 )

f1 <-  cconv(upsampling(a,1),reverse(h),1) + cconv(upsampling(d,1),reverse(g),1)

print( norm(f-f1)/norm(f) )

n <- 256
name <- 'nt_toolbox/data/flowers.png'
f <- load_image(name, n)

imageplot(f, 'Original image')

j = as.integer(log(n,2)-1)

fW <- f

A <- fW[1:2**(j+1),1:2**(j+1)]

Coarse <- subsampling(cconv(A,h,1),1)
Detail <- subsampling(cconv(A,g,1),1)

A <- rbind(Coarse, Detail)

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

imageplot(f, 'Original image', c(1,2,1))
imageplot(A, 'First step', c(1,2,2))

Coarse = subsampling(cconv(A,h,2),2)
Detail = subsampling(cconv(A,g,2),2)

A <- cbind( Coarse, Detail )

fW[1:2**(j+1),1:2**(j+1)] <- A

plot_wavelet(fW,j)

options(repr.plot.width=7, repr.plot.height=7)
layout(matrix(1:12, c(3,4), byrow = FALSE))

Jmin <- 0
Jmax <- as.integer(log(n,2)-1)
fW <- f
for ( j in Jmax:Jmin ){
    A <- fW[1:2**(j+1), 1:2**(j+1)]
    for (d in c(1,2)){
        Coarse <- subsampling(cconv(A,h,d),d)
        Detail <- subsampling(cconv(A,g,d),d)
        if (d==1){ A <- rbind( Coarse, Detail ) }
        else{ A <- cbind( Coarse, Detail ) } }
    fW[1:2**(j+1), 1:2**(j+1)] <- A
    j1 <- Jmax-j
    if (j1<4){
        imageplot(A[1:2**j, (2**j + 1):2**(j+1)], paste("H,j=",j,sep="") )
        imageplot(A[(2**j + 1):2**(j+1),1:2**j], paste("V,j=",j,sep="") )
        imageplot(A[(2**j + 1):2**(j+1),(2**j + 1):2**(j+1)], paste("D,j=",j,sep=""))
    }
}

paste( 'Energy of the signal/coefficients = ', (norm(f)/norm(fW)) )

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

f1 <- fW
j <- 0

A <- f1[1:2**(j+1),1:2**(j+1)]

Coarse <- A[1:2**j, , drop=F]
Detail <- A[(2**j + 1):2**(j+1), , drop=F]

Coarse <- cconv(upsampling(Coarse,1),reverse(h),1)
Detail <- cconv(upsampling(Detail,1),reverse(g),1)

A <- Coarse + Detail

Coarse <- A[,1:2**j, drop=F]
Detail <- A[,(2**j + 1):2**(j+1), drop=F]

Coarse <- cconv(upsampling(Coarse,2),reverse(h),2)
Detail <- cconv(upsampling(Detail,2),reverse(g),2)

A <- Coarse + Detail

f1[1:2**(j+1),1:2**(j+1)] <- A

options(repr.plot.width=5, repr.plot.height=5)
layout(matrix(4:1, c(2,2), byrow = TRUE))

f1 <- fW
for (j in Jmin:Jmax){
    A <- f1[1:2**(j+1),1:2**(j+1)]
    for (d in 1:2){
        if (d==1){
            Coarse <- A[1:2**j, , drop = F]
            Detail <- A[(2**j + 1):2**(j+1), , drop = F]
        }
        else{
            Coarse <- A[,1:2**j, drop = F]
            Detail <- A[,(2**j + 1):2**(j+1), drop = F]
        }
        Coarse <- cconv(upsampling(Coarse,d),reverse(h),d)
        Detail <- cconv(upsampling(Detail,d),reverse(g),d)
        A <- Coarse + Detail
    }
    j1 <- Jmax-j
    if (j1>0 & j1<5){
        imageplot(A, paste("j=",j,sep=""))
        }
    f1[1:2**(j+1),1:2**(j+1)] <- A
}

paste( 'Error |f-f1|/|f| = ' , norm(f-f1)/norm(f) ) 

fW <- perform_wavortho_transf(f,Jmin,+1,h)
eta <- 4
fWLin <- matrix(0,n,n)
fWLin[1:n/eta,1:n/eta] <- fW[1:n/eta,1:n/eta]
fLin <- perform_wavortho_transf(fWLin,Jmin,-1,h)

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

elin <- snr(f,fLin)
imageplot(clamp(fLin), paste('Linear, SNR=',elin))

T <- 0.2

fWT <- fW * (abs(fW)>T)

options(repr.plot.width=5, repr.plot.height=3)
layout(matrix(1:2, c(1,2), byrow = TRUE))

plot_wavelet(fW,Jmin)
plot_wavelet(fWT,Jmin)

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

m <- round((n**2)/(eta**2))
v <- reverse( sort( ravel(abs(fW)) ) )
plot(log(v, 10), type="l", col="blue", style = "d")

T <- v[m]

fWT <- fW * (abs(fW)>=T)
fnl <- perform_wavortho_transf(fWT,Jmin,-1,h)
enl <- snr(f,fnl)
imageplot(clamp(fnl), paste('Non-linear, SNR=', enl))
<script>
  $(document).ready(function(){
      $('div.prompt').hide();
  });
</script>