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.
using PyPlot
using NtToolBox
using Autoreload
arequire("NtToolBox")
areload()
3D Volumetric Datasets¶
We load a volumetric data.
M = NtToolBox.read_bin("NtToolBox/src/data/vessels.bin", 3);
M = NtToolBox.rescale(M);
Size of the image (here it is a cube).
n = size(M)[2];
We can display some horizontal slices.
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.
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
MW = copy(M);
Average/difference along X
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
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
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.
figure(figsize = (10, 5))
imageplot(MW[:, :, 30], "Horizontal slice", [1,2,1])
imageplot((MW[:, 30, :]), "Vertical slice", [1,2,2])
Exercise 1
Implement the forward wavelet transform by iteratively applying these transform steps to the low pass residual.
include("NtSolutions/multidim_2_volumetric/exo1.jl")
## 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
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.
include("NtSolutions/multidim_2_volumetric/exo2.jl")
## Insert your code here.
Display the approximation as slices.
s = 30
figure(figsize = (10, 5))
imageplot(M[:, :, s], "Original", [1, 2, 1])
imageplot(clamP(M1[:, :, s]), "Approximation", [1,2,2])
Display the approximated isosurface.
isosurface(M1, .5, 2, "")
Linear Volumetric Denoising ¶
Linear denoising is obtained by low pass filtering.
We add a Gaussian noise to the image.
sigma = .06
Mnoisy = M + sigma.*randn(n, n, n);
Display slices of the noisy data.
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])