Volumetric Radon Inversion

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 the reconstruction from 3D tomographic measurement with TV regularization.

In [2]:
addpath('toolbox_signal')
addpath('toolbox_general')
addpath('solutions/multidim_6_tomography')

3D Volumetric Datasets

We load a volumetric data.

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

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

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

Reduce dimensionality

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

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 Tomograpic Measurements

Tomographic measurements corresponds to the computation of orthogonal projections (integration along liners) of the 3D datasets on set of 2D planes. Thanks to the Fourier-slice theorem, this is equivalent to performing a sub-sampling of the 3D Fourier transform along planes (orthogonal to the projection directions).

Number of projections

In [ ]:
P = 12;

Either uniform sampling for P==12, or randomized projection directions, on the sphere.

In [ ]:
if P==12
    tau = 0.8506508084;
    one = 0.5257311121;
    S = [  tau,  one,    0;
        -tau,  one,    0
        -tau, -one,    0;
        tau, -one,    0;
        one,   0 ,  tau;
        one,   0 , -tau;
        -one,   0 , -tau;
        -one,   0 ,  tau;
        0 ,  tau,  one;
        0 , -tau,  one;
        0 , -tau, -one;
        0 ,  tau, -one ]';
else
    S = randn(3,P);
    S = S ./ repmat( sqrt( sum(S.^2,1) ), [3 1]);
end

Display the directions on the sphere.

In [ ]:
clf;
plot3(S(1,:), S(2,:), S(3,:), '.');
axis equal;

The Fourier mask.

In [ ]:
x = [0:n/2-1, -n/2:-1];
[X,Y,Z] = ndgrid(x,x,x);
mask = zeros(n,n,n);
epsilon = .5;
for i=1:P
    q = S(:,i);
    d = q(1)*X + q(2)*Y + q(3)*Z;
    mask( abs(d)<=epsilon ) = 1;
end

Tomographic measurement can thus be intepreted as a selection of a few Fourier frequencies.

In [ ]:
F = fftn(M);
y = F(mask==1);

Number of measures.

In [ ]:
Q = length(y);
disp(strcat(['Number of measurements Q=' num2str(Q) '.']));
disp(strcat(['Sub-sampling Q/N=' num2str(length(y)/n^3,2) '.']));

The transposed operator corresponds to the pseudo inverse reconstruction (because the measurement operator is in fact an orthogonal projection). It is similar to the filtered back-projection (excepted that the Fourier sub-sampling is now on a discrete grid, which is not really faithful to the geometry of tomographic acquisition).

In [ ]:
F1 = zeros(n,n,n);
F1(mask==1) = y;
M1 = real( ifftn(F1) );

Display.

In [ ]:
clf;
for i=1:length(slices)
    s = slices(i);
    imageplot( clamp(M1(:,:,s)), strcat(['Z=' num2str(s)]), 2,2,i );
end

SNR of the pseudo inverse reconstruction.

In [ ]:
disp(['Pseudo-inverse, SNR=' num2str(snr(M,M1),4) 'dB.']);

3D Total Variation

Since medical images often have a cartoon morphology, it makes sense to perform the reconstruction while minimizing the TV norm of the image.

To be able to use a gradient descent for the reconstruction, we use a smoothed TV-norm, using a smoothing parameter epsilon. The smaller, the closer to TV, but the slower the method is.

In [ ]:
epsilon = 1e-2;

The gradient of the volum is a 3D vector field, and hence a |(n,n,n,3)| matrix.

In [ ]:
G = grad(M);

Display the gradient as a color image.

In [ ]:
clf;
for i=1:length(slices)
    s = slices(i);
    imageplot( squeeze(G(:,:,s,:)), strcat(['Z=' num2str(s)]), 2,2,i );
end

Compute the smoothed norm of the gradient. This corresponds to a edge detector.

In [ ]:
d = sqrt(sum(G.^2,4)+epsilon^2);

Display.

In [ ]:
clf;
for i=1:length(slices)
    s = slices(i);
    imageplot( squeeze(d(:,:,s,:)), strcat(['Z=' num2str(s)]), 2,2,i );
end

Compute the regularized TV norm.

In [ ]:
tv = sum(d(:));
disp(['TV norm=' num2str(tv) '.']);

TV Reconstruction from Partial 3D Tomoraphic Measurements

We perform the reconstruction using a projected gradient descent.

Gradient descent step. Should be proportional to |epsilon|.

In [ ]:
tau = epsilon*.2;

Initialize the solution using the pseudo inverse.

In [ ]:
Mtv = M1;

The gradient of the smoothed total variation is minus the divergence of the normalized gradient of the image.

In [ ]:
G = grad(Mtv);
d = sqrt(sum(G.^2,4)+epsilon^2);
dG = -div( G ./ repmat(d, [1 1 1 3]) );

Display it.

In [ ]:
clf;
for i=1:length(slices)
    s = slices(i);
    imageplot( dG(:,:,s), strcat(['Z=' num2str(s)]), 2,2,i );
end

Perform the gradient step.

In [ ]:
Mtv = Mtv - tau*dG;

Perform the projection step to impose the known values.

In [ ]:
F = fftn(Mtv);
F(mask==1) = y;
Mtv = real(ifftn(F));

Exercise 1

Perform the projected TV gradient desent, and monitor the epsilon-TV norm during the iterations. isplay TV decay

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

Display.

In [ ]:
clf;
for i=1:length(slices)
    s = slices(i);
    imageplot( clamp(Mtv(:,:,s)), strcat(['Z=' num2str(s)]), 2,2,i );
end

SNR of the resulting reconstruction.

In [ ]:
disp(['Total variation, SNR=' num2str(snr(M,Mtv),4) 'dB.']);