Volumetric wavelet Data Processing

Important: Please read the installation page for details about how to install the toolboxes. $\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\si}{\sigma}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$

This numerical tour explores volumetric (3D) data processing.

In [3]:
options(warn=-1) # turns off warnings, to turn on: "options(warn=0)"


library(imager)
library(png)
library(misc3d)
library(SynchWave)

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)

3D Volumetric Datasets

We load a volumetric data.

In [4]:
M <- read_bin("nt_toolbox/data/vessels.bin", ndims=3)
In [5]:
M <- rescale(M)

Size of the image (here it is a cube).

In [6]:
n <- dim(M)[2]

We can display some horizontal slices.

In [7]:
options(repr.plot.width=6, repr.plot.height=6)

slices <- seq(10,n-10, length = 4)
for (i in 1:length(slices)){
    s <- slices[i]
    s <- as.integer(s)
    imageplot(M[,,s], paste("Z = ",s), c(2,2,i)) }

We can display an isosurface of the dataset (here we sub-sample to speed up the computation).

In [8]:
options(repr.plot.width=8, repr.plot.height=6)

isosurface(M,0.5,3)

3D Haar Transform

An isotropic 3D Haar transform recursively extracts details wavelet coefficients by performing local averages/differences along the X/Y/Z axis.

We apply a step of Haar transform in the X/Y/Z direction

Initialize the transform

In [9]:
MW <- M

Average/difference along X

In [10]:
a <- (MW[seq(1,n,2),,] + MW[seq(2,n,2),,])/sqrt(2)
b <- (MW[seq(1,n,2),,] - MW[seq(2,n,2),,])/sqrt(2)
c <- array(0, c(n, n, n))
c[1:(n/2),,] <- a
c[(n/2+1):n,,] <- b
MW <- c

Average/difference along Y

In [11]:
a <- (MW[,seq(1,n,2),] + MW[,seq(2,n,2),])/sqrt(2)
b <- (MW[,seq(1,n,2),] - MW[,seq(2,n,2),])/sqrt(2)
c <- array(0, c(n, n, n))
c[,1:(n/2),] <- a
c[,(n/2+1):n,] <- b
MW <- c

Average/difference along Z

In [12]:
a <- (MW[,,seq(1,n,2)] + MW[,,seq(2,n,2)])/sqrt(2)
b <- (MW[,,seq(1,n,2)] - MW[,,seq(2,n,2)])/sqrt(2)
c <- array(0, c(n, n, n))
c[,,1:(n/2)] <- a
c[,,(n/2+1):n] <- b
MW <- c

Display a horizontal and vertical slice to see the structure of the coefficients.

In [13]:
options(repr.plot.width=8, repr.plot.height=6)

imageplot(MW[,,30], "Horizontal slice", c(1,2,1))
imageplot((MW[,30,]), "Vertical slice", c(1,2,2))

Exercise 1

Implement the forward wavelet transform by iteratively applying these transform steps to the low pass residual.

In [14]:
source("nt_solutions/multidim_2_volumetric/exo1.R")
In [15]:
## Insert your code here.

Volumetric Data Haar Approximation

An approximation is obtained by keeping only the largest coefficients.

We threshold the coefficients to perform $m$-term approximation.

number of kept coefficients

In [16]:
m <- round(0.01*n**3)
MWT <- perform_thresholding(MW, m, type="largest")

Exercise 2

Implement the backward transform to compute an approximation $M_1$ from the coefficients MWT.

In [17]:
source("nt_solutions/multidim_2_volumetric/exo2.R")

Display the approximation as slices.

In [18]:
s <- 30

imageplot(M[, ,s], 'Original', c(1,2,1))
imageplot(clamp(M1[, ,s]), 'Approximation', c(1,2,2))

Display the approximated isosurface.

In [19]:
isosurface(M1,.5,2)

Linear Volumetric Denoising

Linear denoising is obtained by low pass filtering.

We add a Gaussian noise to the image.

In [20]:
sigma <- 0.06
Mnoisy <- M + array(rnorm(n**3, sd = sigma), c(n,n,n))

Display slices of the noisy data.

In [21]:
imageplot(Mnoisy[,,n/2],"X slice",c(1,2,1))
imageplot(Mnoisy[,n/2,],"Y slice",c(1,2,2))