1.1 Load and plot a signal from the MakeSignal function

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)

1.2 Plot a noisy version of the original function

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))

1.3 Forward Wavelet Transform (periodized, orthogonal)

qmf <- MakeONFilter('Daubechies', 10)
j0 <- 3
wc <- FWT_PO(y, j0, qmf)
wch <- wcs <- wcbs <- wcbh <- wc
wcf <- FWT_PO(f, j0, qmf)

1.4 Spike-plot display of wavelet coefficients

PlotWaveCoeff(wc, j0, 0.5)
title("Noisy Wavelet Coefs")

PlotWaveCoeff(wcf, j0, 0.5)
title("Original Wavelet Coefs")

1.5 Universal hard thresholding

# 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)

1.6 Universal hard thresholding (Block version)

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)

1.7 Universal soft thresholding

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)

1.8 Universal soft thresholding (Block version)

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)

2.1 Load and plot an image

library(imager)
name <- '../inst/extdata/lena.png'
f <- load.image(name)
plot(f, axes=F, interpolate=F, xlab="", ylab="")

2.2 Plot a noisy version of the original image

ssig <- sd(f)
sdnoise <- ssig/SNR
y <- f + rnorm(ncol(f)*nrow(f), mean=0, sd=sdnoise)
snrin <- SNR(f,y)
plot(y, axes=F, interpolate=F, xlab="", ylab="",
     main=format(round(snrin,2), nsmall = 2))

2.3 Image denoising using hard thresholding

wc <- FWT2_PO(as.array(squeeze(y)), j0, qmf)
wcb <- wc
thr <- 3*sdnoise
aT <- wc*(abs(wc)>thr)
fest <- IWT2_PO(aT, j0, qmf)
snrout <- SNR(f, fest)
plot(as.cimg(fest), axes=FALSE, xlab="", ylab="",
     main=format(round(snrout,2), nsmall=2))

3.1 Load and plot a 3d image

library(misc3d)
library(rgl)
data("SLphantom")
n <- dim(SLphantom)[1]
contour3d(SLphantom,0)
rglwidget()

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)

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