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 [2]:
addpath('toolbox_signal')
addpath('toolbox_general')
addpath('solutions/multidim_2_volumetric')

For Scilab user (otherwise there is not enough memory). WARNING: you should extend the stack size only ONCE.

In [3]:
extend_stack_size(10);

3D Volumetric Datasets

We load a volumetric data.

In [4]:
name = 'vessels';
options.nbdims = 3;
M = read_bin(name, options);
M = rescale(M);

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

In [5]:
n = size(M,1);

We can display some horizontal slices.

In [ ]:
slices = round(linspace(10,n-10,4));
clf;
for i=1:length(slices)
    s = slices(i);
    imageplot( M(:,:,s), strcat(['Z=' num2str(s)]), 2,2,i );
end

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

In [ ]:
sel = 1:2:n;
clf;
isosurface( M(sel,sel,sel), .5);
axis('off');

3D Haar Transform

An isotropic 3D Haar transform recursively extract 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 [ ]:
MW = M;

average/difference along X

In [ ]:
MW = cat3(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 [ ]:
MW = cat3(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 [ ]:
MW = cat3(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 [ ]:
clf;
imageplot(MW(:,:,30), 'Horizontal slice', 1,2,1);
imageplot(squeeze(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. nitialize the transform

In [ ]:
exo1()
In [ ]:
%% 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 [ ]:
m = round(.01*n^3);
MWT = perform_thresholding(MW, m, 'largest');

Exercise 2

Implement the backward transform to compute an approximation |M1| from the coefficients |MWT|.

In [ ]:
exo2()
In [ ]:
%% Insert your code here.

Display the approximation as slices.

In [ ]:
s = 30;
clf;
imageplot( M(:,:,s), 'Original', 1,2,1 );
imageplot( clamp(M1(:,:,s)), 'Approximation', 1,2,2 );

Display the approximated isosurface.

In [ ]:
sel = 1:2:n;
clf;
isosurface( M1(sel,sel,sel), .5);
axis('off');

Linear Volumetric Denoising

Linear denoising is obtained by low pass filtering.

We add a Gaussian noise to the image.

noise level.

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

Display slices of the noisy data.

In [ ]:
clf;
imageplot(Mnoisy(:,:,n/2),'X slice', 1,2,1);
imageplot(squeeze(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 [ ]:
x = -n/2:n/2-1;
[X,Y,Z] = ndgrid(x,x,x);

gaussian filter

In [ ]:
s = 2; % width
h = exp( -(X.^2 + Y.^2 + Z.^2)/(2*s^2) );
h = h/sum(h(:));

The filtering is computed over the Fourier domain.

In [ ]:
Mh = real( ifftn( fftn(Mnoisy) .* fftn(fftshift(h)) ) );

Display denoised slices.

In [ ]:
i = 40;
clf;
imageplot( Mnoisy(:,:,i), 'Noisy', 1,2,1 );
imageplot( Mh(:,:,i), 'Denoised', 1,2,2 );

Display denoised iso-surface.

In [ ]:
sel = 1:2:n;
clf;
isosurface( Mh(sel,sel,sel), .5);
axis('off');

Exercise 3

Select the optimal blurring width |s| to reach the smallest possible SNR. Keep the optimal denoising |Mblur|

In [ ]:
exo3()
In [ ]:
%% Insert your code here.

Display optimally denoised iso-surface.

In [ ]:
sel = 1:2:n;
clf;
isosurface( Mblur(sel,sel,sel), .5);
axis('off');
title(['Filtering, SNR=' num2str(snr(M,Mblur),3) 'dB']);

Non-Linear Wavelet Volumetric Denoising

Denoising is obtained by removing small amplitude coefficients that corresponds to noise.

Exercise 4

Perforn Wavelet denoising by thresholding the wavelet coefficients of Mnoisy. Test both hard thresholding and soft thresholding to determine the optimal threshold and the corresponding SNR. Record the optimal result |Mwav|.

In [ ]:
exo4()
In [ ]:
%% Insert your code here.

Display denoised iso-surface with optimal soft thresholding.

In [ ]:
sel = 1:2:n;
clf;
isosurface( Mwav(sel,sel,sel), .5);
title(['Soft thresholding, SNR=' num2str(snr(M,Mwav),3) 'dB']);
axis('off');

Orthogonal wavelet thresholdings suffers from blocking artifacts. This can be aleviated by performing a cycle spinning denoising, which average the denosing result of translated version of the signal.

A typical cycle spinning process is like this.

maximum translation

In [ ]:
w = 4;

list of translations

In [ ]:
[dX,dY,dZ] = ndgrid(0:w-1,0:w-1,0:w-1);

initialize spinning process

In [ ]:
Mspin = zeros(n,n,n);

spin

In [ ]:
for i=1:w^3
    % shift the image
    MnoisyC = circshift(Mnoisy, [dX(i) dY(i) dZ(i)]);
    % denoise the image to get a result M1
    M1 = MnoisyC; % replace this line by some denoising
    % shift inverse
    M1 = circshift(M1, -[dX(i) dY(i) dZ(i)]);
    % average the result
    Mspin = Mspin*(i-1)/i + M1/i;
end

Exercise 5

Implement cycle spinning hard thresholding with |T=3*sigma|.

In [ ]:
exo5()
In [ ]:
%% Insert your code here.

Display denoised iso-surface.

In [ ]:
sel = 1:2:n;
clf;
isosurface( Mspin(sel,sel,sel), .5);
title(['Cycle spinning, SNR=' num2str(snr(M,Mspin),3) 'dB']);
axis('off');

Higher Order Daubechies Wavelets

Similarely to the Haar transform, wavelets with more vanishing moments are obtained with filtering and subsampling (pyramidal algorihtm).

First we create the filters for 4 vanishing moments.

In [ ]:
[h,g] = compute_wavelet_filter('Daubechies',2*4);

Then we initialize the wavelet transform with the 3D image itself, and set the current scale.

In [ ]:
MW = M;
j = log2(n)-1;

We transform by filtering + sub-sampling the low pass residual along the three directions.

In [ ]:
A = MW(1:2^(j+1),1:2^(j+1),1:2^(j+1));
for d=1:3
    a = cat(d, subsampling(cconv(A,h,d),d), subsampling(cconv(A,g,d),d) );
end
MW(1:2^(j+1),1:2^(j+1),1:2^(j+1)) = A;

Exercise 6

Implement the full 3D forward wavelet transform by applying these steps for decaying scales |j| toward 0.

In [ ]:
exo6()
In [ ]:
%% Insert your code here.

Threshold the coefficients.

In [ ]:
MWT = perform_thresholding(MW,m,'largest');

Display the coefficients and thresholded coefficients for one slice.

In [ ]:
clf;
subplot(1,2,1);
plot_wavelet(MW(:,:,n/8));
subplot(1,2,2);
plot_wavelet(MWT(:,:,n/8));

Initialize the backward transform.

In [ ]:
M1 = MWT;
j = 0;

Undo one step of the wavelet transform. Note: |subselectdim(A,sel,1)| is equivalent to |A(sel,:,:)| while |subselectdim(A,sel,2)| is equivalent to |A(:,sel,:)|.

In [ ]:
A = M1(1:2^(j+1),1:2^(j+1),1:2^(j+1));
for d=1:3
    W = subselectdim(A,2^j+1:2^(j+1),d);
    A = subselectdim(A,1:2^j,d);
    A = cconv(upsampling(A,d),reverse(h),d) + cconv(upsampling(W,d),reverse(g),d);
end
M1(1:2^(j+1),1:2^(j+1),1:2^(j+1)) = A;

Exercise 7

Implement the full 3D backward wavelet transform by applying these steps for increasing scales |j|.

In [ ]:
exo7()
In [ ]:
%% Insert your code here.

Display approximated iso-surface.

In [ ]:
sel = 1:2:n;
clf;
isosurface( M1(sel,sel,sel), .5);
title(['Soft thresholding, SNR=' num2str(snr(M,M1),3) 'dB']);
axis('off');

Denoising Daubechies Wavelets

Better denoising results are obtined by thresholding orthogonal wavelet coefficients.

Exercise 8

Implement denoising by soft and hard thresholding Daubechies wavelet coefficients.

In [ ]:
exo8()
In [ ]:
%% Insert your code here.

Display denoised iso-surface with optimal soft thresholding.

In [ ]:
sel = 1:2:n;
clf;
isosurface( Mwav(sel,sel,sel), .5);
title(['Soft thresholding, SNR=' num2str(snr(M,Mwav),3) 'dB']);
axis('off');

Exercise 9

Implement cycle spinning hard thresholding with Daubechies wavelets with |T=3*sigma|.

In [ ]:
exo9()
In [ ]:
%% Insert your code here.

Display denoised iso-surface.

In [ ]:
sel = 1:2:n;
clf;
isosurface( Mspin(sel,sel,sel), .5);
title(['Cycle spinning, SNR=' num2str(snr(M,Mspin),3) 'dB']);
axis('off');