hadamard_noise.Rmd

---
title: "Denoising Images"
author: "Derek Wayne"
date: "`r format(Sys.time(), '%d %B %Y')`"
output:
  html_notebook: default
  html_document:
    df_print: paged
  pdf_document: default
subtitle: Applications of the Walsh-Hadamard Transform
abstract: |
  The Hadamard transform is an example of a generalized Fourier transform. It has
  numerous applications in signal processing, data compression, and image/video coding to
  name a few. The Fast Walsh-Hadamard transform (FWHT) is an efficient algorithm for
  denoising images by using various methods of pixel thresholding. In this note, we
  explore how this might be implented on a noisy image.
---

# Hadamard Matrices

For *Hadamard* matrices $H_m$ and $H_n$ we have the following useful recursive properties (note that since m and n are powers of 2, the matrices are symmetric),

\begin{equation}
  H_mH_m^T = 2^kI  
\end{equation}

\begin{equation}
  H_m^{-1} = \frac{1}{m}H_m
\end{equation}

### Proposition 1.
If we let $\hat{Z} = H_mZH_n$, then $Z =  H_m \hat{Z} H_n/(mn)$

*Proof*
$$
\begin{aligned}
\hat{Z} &= H_mZH_n \\
\iff (1/m)H_m \hat{Z} (1/n)H_n &= H_m^{-1}H_m Z H_nH_n^{-1} \\
\iff H_m \hat{Z} H_n/(mn) &= Z
\end{aligned}
$$

Loading in the image as a matrix:

```{r, fig.height=5, fig.width=5}
boats <- matrix(scan("boats.txt"), ncol=256, byrow=T)
image(boats, axes=F, col=grey(seq(0,1,length=256)), main = "Original")
```


Below is an implementation of the FWHT.

```{r}
fwht2d <- function(x) {
  h <- 1
  length <- ncol(x)
  while (h < length) {
    for (i in seq(1, length, by=h*2)) {
      for (j in seq(i,i+h-1)) {
        a <- x[,j]
        b <- x[,j+h]
        x[,j] <- a + b
        x[,j+h] <- a - b
      }
    }
    h <- 2*h
  }
  h <- 1
  length <- nrow(x)
  while (h < length) {
    for (i in seq(1, length, by=h*2)) {
      for (j in seq(i,i+h-1)) {
        a <- x[j, ]
        b <- x[j+h, ]
        x[j, ] <- a + b
        x[j+h, ] <- a - b
      }
    }
    h <- 2*h
  }
  x
}
```


Below are the functions to perform soft and hard thresholding.

```{r}
hard.thresh <- function(x, lambda) {ifelse(abs(x) > lambda, x, 0)}
soft.thresh <- function(x, lambda) {
  sign(x) * (abs(x) >= lambda) * (abs(x) - lambda)
}
```

The function below will divide the original image up into subimages depending on the function parameter "size"; referring to $n$ where $n=2^k, \quad 2\leq k \leq 8$ so that the subimages are n by n pixel sub-images.


```{r}
sub.im.transform <- function(image, size, Th.type="soft", thresh) {
  # image size must be square
  N <- ncol(image)
  
  if ((N %% size) != 0) {
    stop("image size is not divisible by specified sub-image size")
  }
  # create a list of submatrices to apply the transform to
  flat <- as.vector(image)
  subimages <- list()
  A <- N^2/size^2 # index span for each subimage 
  for (j in 0:(A-1)) {
    subimages[[j+1]] <- matrix(flat[(j*(size^2)+1):((j+1)*(size^2))],ncol=size)
  }
  # apply transform to each sub-image and reconstruct the original matrix/image
  subimages.hat <- list()
  i <- 1
  for (im in subimages) {
    Zhat <- fwht2d(im)
    if (Th.type == "soft") {
      Zhat.star <- soft.thresh(Zhat, thresh)
    } else {
      Zhat.star <- hard.thresh(Zhat, thresh)
    }
    Zstar <- fwht2d(Zhat.star) / size
    subimages.hat[[i]] <- Zstar
    i <- i + 1
  }
  
  # reconstruct original image
  original <- c()
  for (i in 1:A) {
    original <- c(original, as.vector(subimages.hat[[i]]))
  }
  original <- matrix(original, ncol = 256)
  original
}
```

# Results

```{r, fig.align="center", echo=F, fig.height=5}
par(mfrow=c(1,2))
soft <- sub.im.transform(boats, 64, Th.type = "soft", thresh = 7)
image(boats, axes=F, col=grey(seq(0,1,length=256)), main="Original")
image(soft, axes=F, col=grey(seq(0,1,length=256)), main="Soft Threshold: lam=7")
```

```{r, fig.align="center", echo=F, fig.height=5}
par(mfrow=c(1,2))
hard <- sub.im.transform(boats, 32, Th.type = "hard", thresh = 6)
image(boats, axes=F, col=grey(seq(0,1,length=256)), main="Original")
image(hard, axes=F, col=grey(seq(0,1,length=256)), main="Hard Threshold: lam=6")
```

```{r, fig.align="center", echo=F, fig.height=5}
par(mfrow=c(1,2))
soft2 <- sub.im.transform(boats, 8, Th.type = "soft", thresh = 0.8)
image(boats, axes=F, col=grey(seq(0,1,length=256)), main="Original")
image(soft2, axes=F, col=grey(seq(0,1,length=256)), main="Soft Threshold: lam=0.8")
```