---
title: 'Rwavelet Vignette'
author: "Christophe Chesneau and Fabien Navarro"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_document:
    theme: cerulean
    number_sections: yes
    toc: yes
    toc_depth: 5
    toc_float: true
vignette: >
  %\VignetteIndexEntry{Rwavelet Vignette}
  %\VignetteEngine{knitr::rmarkdown}
  %\usepackage[utf8]{inputenc}
---
<style type="text/css">

body{ /* Normal  */
      font-size: 12px;
  }
td {  /* Table  */
  font-size: 8px;
}
h1.title {
  font-size: 38px;
  color: DarkRed;
}
h1 { /* Header 1 */
  font-size: 28px;
  color: DarkBlue;
}
h2 { /* Header 2 */
    font-size: 18px;
  color: DarkBlue;
}
h3 { /* Header 3 */
  font-size: 15px;
  font-family: "Times New Roman", Times, serif;
  color: DarkBlue;
}
code.r{ /* Code block */
    font-size: 12px;
}
pre { /* Code block - determines code spacing between lines */
    font-size: 14px;
}
</style>

```{r, echo=FALSE}
#options(vignetteDocumentFormat=rmarkdown::all_output_formats("rwaveletvignette.Rmd"))
```

```{r Setup, include=FALSE}
library(rwavelet)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width=5,
  fig.height=5, 
  fig.align="center"
)
set.seed(0)
par(mgp = c(2,0.5,0),
    mar = c(4,3,2,1),
    oma = c(1,1,0,0))
```

# 1d examples

## Load and plot a signal from the MakeSignal function {#example1}
```{r, echo=TRUE}
J <- 10; n <- 2^J; t <- (1:n) / n
name <- c('Bumps')
f <- MakeSignal(name, n)
plot(t, f, xlab="t", ylab="f(t)",
     type='l', lwd=1.2, main=name)
```

## Plot a noisy version of the original function {#example2}
```{r, echo=TRUE}
SNR <- 4
set.seed(1)
ssig <- sd(f)
sigma <- ssig / SNR
y <- f + rnorm(n, mean=0, sd=sigma)
plot(t, y, xlab="t", ylab="y", main=paste("Noisy", name))
```

## Forward Wavelet Transform (periodized, orthogonal) {#example3}
```{r, echo=TRUE}
qmf <- MakeONFilter('Daubechies', 10)
j0 <- 3
wc <- FWT_PO(y, j0, qmf)
wch <- wcs <- wcbs <- wcbh <- wc
wcf <- FWT_PO(f, j0, qmf)
```

## Spike-plot display of wavelet coefficients {#example4}
```{r, echo=TRUE}
PlotWaveCoeff(wc, j0, 0.5)
title("Noisy Wavelet Coefs")
PlotWaveCoeff(wcf, j0, 0.5)
title("Original Wavelet Coefs")
```

## Universal hard thresholding 
```{r, echo=TRUE}
# estimate sigma using the Median Absolute Deviation
# using only the fine scale of wc
hatsigma <- MAD(wc[(2^(J-1)+1):2^J])
thr <- sqrt(2*log(length(y)))*hatsigma
# apply hard thresholding 
wch[(2^(j0)+1):n] <- HardThresh(wc[(2^(j0)+1):n], thr)
# plot the resulting coeficients estimates
PlotWaveCoeff(wch, j0, 0.5)
title("Estimated Wavelet Coefs")
fest <- IWT_PO(wch, j0, qmf)
snrout <- SNR(f, fest)
plot(t, fest, type='l', lwd=1.4, col='red', xlab="t", ylab="hat_f(t)",
	   main=format(round(snrout,2), nsmall=2))
matlines(t, f, type='l', lty=2)
```

## Universal hard thresholding (Block version)
```{r}
L <- 2 # block size
wcbh <- BlockThresh(wc, j0, hatsigma, L, qmf, thresh = "hard")
fest_bh <- IWT_PO(wcbh,j0,qmf)
PlotWaveCoeff(wcbh, j0, 0.5)
title("Estimated Wavelet Coefs")
plot(t, fest_bh, type='l', lwd=1.4, col='red', xlab="t", ylab="hat_f(t)",
	   main=format(round(SNR(f, fest_bh),2), nsmall=2))
matlines(t, f, type='l', lty=2)
```

## Universal soft thresholding 
```{r, echo=TRUE}
wcs[(2^(j0)+1):n] <- SoftThresh(wc[(2^(j0)+1):n], thr)
PlotWaveCoeff(wcs, j0, 0.5)
title("Estimated Wavelet Coefs (soft)")
f_soft <- IWT_PO(wcs, j0, qmf)
plot(t, f_soft, type='l', lwd=1.4, col='red', xlab="t", ylab="hat_f(t)",
	   main=format(round(SNR(f, f_soft),2), nsmall=2))
matlines(t, f, type='l', lty=2)
```

## Universal soft thresholding (Block version)
```{r}
wcbs <- BlockThresh(wc, j0, hatsigma, L, qmf, thresh = "soft")
PlotWaveCoeff(wcbs, j0, 0.5)
fest_bs <- IWT_PO(wcbs,j0,qmf)
plot(t, fest_bs, type='l', lwd=1.4, col='red', xlab="t", ylab="hat_f(t)",
	   main=format(round(SNR(f, fest_bs),2), nsmall=2))
matlines(t, f, type='l', lty=2)
```

# 2d examples
## Load and plot an image

```{r, echo=TRUE, message=FALSE}
name <- '../inst/extdata/lena.png'
if (requireNamespace("imager", quietly = TRUE)) {
      f <- imager::load.image(name)
      plot(f, axes=F, interpolate=F, xlab="", ylab="")
   } else {
      ## use e.g. image from graphics package
   }
```

## Plot a noisy version of the original image
```{r, echo=TRUE, message=FALSE}
ssig <- sd(f)
sdnoise <- ssig/SNR
y <- f + rnorm(ncol(f)*nrow(f), mean=0, sd=sdnoise)
snrin <- SNR(f,y)
if (requireNamespace("imager", quietly = TRUE)) {
      plot(y, axes=F, interpolate=F, xlab="", ylab="",
     main=format(round(snrin,2), nsmall = 2))
   } else {
      ## use e.g. image from graphics package
   }
```

## Image denoising using hard thresholding
```{r, echo=TRUE, message=FALSE}
wc <- FWT2_PO(as.array(imager::squeeze(y)), j0, qmf)
wcb <- wc
thr <- 3*sdnoise
aT <- wc*(abs(wc)>thr)
fest <- IWT2_PO(aT, j0, qmf)
snrout <- SNR(f, fest)
if (requireNamespace("imager", quietly = TRUE)) {
      plot(imager::as.cimg(fest), axes=FALSE, xlab="", ylab="",
     main=format(round(snrout,2), nsmall=2))
   } else {
      ## use e.g. image from graphics package
   }
```


# 3d example
## Load and plot a 3d image

```{r, eval=FALSE}
library(misc3d)
library(rgl)
data("SLphantom")
n <- dim(SLphantom)[1]
contour3d(SLphantom,0)
rglwidget()
```

```{r, eval=FALSE}
sig <- sd(SLphantom)
sdnoise <- sig/SNR
y <- SLphantom + rnorm(n^3, mean=0, sd=sdnoise)
wcf <- FWT3_PO(SLphantom,j0,qmf)
wc <- FWT3_PO(y,j0,qmf)
wcn <- wc
thr <- 3*sdnoise
wc <- wc*(abs(wc)>thr)
fhat <- IWT3_PO(wc, j0, qmf)
```


```{r, eval=FALSE}
op <- par(mfrow=c(3,3), mgp = c(1.2, 0.5, 0), tcl = -0.2,
mar = .1 + c(0.1,0.1,0.1,0.1), oma = c(0,0,0,0))
ll=screen=list(z = 130, x = -80)
contour3d(SLphantom,0,color="gray", engine = "standard")
contour3d(y,0.05,color="gray", engine = "standard")
contour3d(fhat,0.2,color="gray", engine = "standard")
contour3d(wcf,0.1, color="gray", engine = "standard")
contour3d(wcn,0.05,color="gray", engine = "standard")
contour3d(wc,0.1,color="gray", engine = "standard")
plot(as.cimg(SLphantom[n/2,,]), axes=FALSE, xlab="", ylab="")
plot(as.cimg(y[n/2,,]), axes=FALSE, xlab="", ylab="")
plot(as.cimg(fhat[n/2,,]), axes=FALSE, xlab="", ylab="")
```