# Volumetric wavelet Data Processing¶


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¶

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