# Volumetric wavelet Data Processing¶


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

In [2]:
using PyPlot
using NtToolBox
#arequire("NtToolBox")


## 3D Volumetric Datasets¶

In [2]:
M = NtToolBox.read_bin("NtToolBox/src/data/vessels.bin", 3);

In [3]:
M = NtToolBox.rescale(M);


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

In [4]:
n = size(M)[2];


We can display some horizontal slices.

In [5]:
slices = Array{Int64,1}(round(linspace(10, n-10, 4)))
figure(figsize = (10,10))

for i in 1:length(slices)
s = slices[i]
NtToolBox.imageplot(M[:, :, s], @sprintf("Z = %i", s), [2, 2, i])
end


We can display an isosurface of the dataset (here we sub-sample to speed up the computation). You need to have PyCall package installed in Julia in order to use some functions in python, besides you also have to install the package skimage in python.

In [6]:
include("NtToolBox/src/isosurface.jl")
isosurface(M, .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 [7]:
MW = copy(M);


Average/difference along X

In [8]:
MW = cat(1, (MW[1:2:n, :, :] + MW[2:2:n, :, :])./sqrt(2), (MW[1:2:n, :, :] - MW[2:2:n, :, :])./sqrt(2) );


Average/difference along Y

In [9]:
MW = cat(2, (MW[:, 1:2:n, :] + MW[:, 2:2:n, :])./sqrt(2), (MW[:, 1:2:n, :] - MW[:, 2:2:n, :])./sqrt(2) );


Average/difference along Z

In [10]:
MW = cat(3, (MW[:, :, 1:2:n] + MW[:, :, 2:2:n])./sqrt(2), (MW[:, :, 1:2:n] - MW[:, :, 2:2:n])./sqrt(2) );


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

In [11]:
figure(figsize = (10, 5))
imageplot(MW[:, :, 30], "Horizontal slice", [1,2,1])
imageplot((MW[:, 30, :]), "Vertical slice", [1,2,2])

Out[11]:
PyObject <matplotlib.text.Text object at 0x000000002F2CD208>

Exercise 1

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

In [12]:
include("NtSolutions/multidim_2_volumetric/exo1.jl")

In [13]:
## 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 [13]:
m = Int(round(.01*n^3))
MWT = NtToolBox.perform_thresholding(MW, m, "largest");


Exercise 2

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

In [14]:
include("NtSolutions/multidim_2_volumetric/exo2.jl")

In [16]:
## Insert your code here.


Display the approximation as slices.

In [15]:
s = 30

figure(figsize = (10, 5))
imageplot(M[:, :, s], "Original", [1, 2, 1])
imageplot(clamP(M1[:, :, s]), "Approximation", [1,2,2])

Out[15]:
PyObject <matplotlib.text.Text object at 0x000000002F36B160>

Display the approximated isosurface.

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


## Linear Volumetric Denoising ¶

Linear denoising is obtained by low pass filtering.

We add a Gaussian noise to the image.

In [18]:
sigma = .06
Mnoisy = M + sigma.*randn(n, n, n);


Display slices of the noisy data.

In [19]:
figure(figsize = (10, 5))
imageplot(Mnoisy[:, :, Base.div(n, 2)], "X slice", [1,2,1])
imageplot(Mnoisy[:, Base.div(n, 2), :], "Y slice", [1,2,2])