# Volumetric wavelet Data Processing¶


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

In [1]:
from __future__ import division

import numpy as np
import scipy as scp
import pylab as pyl
import matplotlib.pyplot as plt

from nt_toolbox.general import *
from nt_toolbox.signal import *

import warnings
warnings.filterwarnings('ignore')

%matplotlib inline


## 3D Volumetric Datasets¶

In [2]:
from nt_toolbox.read_bin import *

In [3]:
M = rescale(M)


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

In [4]:
n = np.shape(M)[1]


We can display some horizontal slices.

In [5]:
slices = np.round(np.linspace(10, n-10, 4))
plt.figure(figsize = (10,10))

for i in range(len(slices)):
s = slices[i]
imageplot(M[:,:,s], "Z = %i" %s, [2,2,i+1])


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

In [6]:
from nt_toolbox.isosurface import *
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 = np.copy(M)


Average/difference along X

In [8]:
MW = np.concatenate(((MW[0:n:2,:,:] + MW[1:n:2,:,:])/np.sqrt(2), (MW[0:n:2,:,:] - MW[1:n:2,:,:])/np.sqrt(2)),0)


Average/difference along Y

In [9]:
MW = np.concatenate(((MW[:,0:n:2,:] + MW[:,1:n:2,:])/np.sqrt(2), (MW[:,0:n:2,:] - MW[:,1:n:2,:])/np.sqrt(2)),1)


Average/difference along Z

In [10]:
MW = np.concatenate(((MW[:,:,0:n:2] + MW[:,:,1:n:2])/np.sqrt(2), (MW[:,:,0:n:2] - MW[:,:,1:n:2])/np.sqrt(2)),2)


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

In [11]:
plt.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.

In [12]:
run -i nt_solutions/multidim_2_volumetric/exo1

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 [14]:
from nt_toolbox.perform_thresholding import *

m = round(.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 [15]:
run -i nt_solutions/multidim_2_volumetric/exo2

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


Display the approximation as slices.

In [17]:
s = 30

plt.figure(figsize=(10,5))
imageplot(M[:, :,s], 'Original', [1,2,1])
imageplot(clamp(M1[:,:,s]), 'Approximation', [1,2,2])


Display the approximated isosurface.

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


## Linear Volumetric Denoising ¶

Linear denoising is obtained by low pass filtering.

We add a Gaussian noise to the image.

In [19]:
from numpy import random

sigma = .06
Mnoisy = M + sigma*random.randn(n, n, n)


Display slices of the noisy data.

In [20]:
plt.figure(figsize=(10,5))
imageplot(Mnoisy[:,:,n//2],"X slice",[1,2,1])
imageplot(Mnoisy[:,n//2,:],"Y slice",[1,2,2])


A simple denoising method performs a linear filtering of the data.

We build a Gaussian filter of width $\sigma$.

Construct a 3D grid

In [21]:
x = np.arange(-n//2,n//2)
[X, Y, Z] = np.meshgrid(x, x, x)


Gaussian filter

In [22]:
s = 2 #width
h = np.exp(-(X**2 + Y**2 + Z**2)/(2*s**2))
h = h/np.sum(h)


The filtering is computed over the Fourier domain.

In [23]:
Mh = np.real(pyl.ifft2(pyl.fft2(Mnoisy,axes=(0,1,2)) * pyl.fft2(pyl.fftshift(h,axes= (0,1,2)),axes=(0,1,2)),axes= (0,1,2)))


Display denoised slices.

In [24]:
i = 40
plt.figure(figsize=(10,5))
imageplot(Mnoisy[:,:,i], "Noisy", [1,2,1])
imageplot(Mh[:,:,i], "Denoised", [1,2,2])


Display denoised iso-surface.

In [25]:
isosurface(M,.5,3)