Tomography Inversion using Tikhonov and Sparse Regularization¶

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 tomographic measurement with Sobolev and sparse regularization.

Tomography and Regularization¶

We consider here tomography measurements $ y = \Phi f_0 + w $, where $f_0$ is the (unknown) image to recover, and $w$ is an additive noise.

The tomography operator $\Phi$ is ill-posed, so it cannot be inverted.

We consider variational reconstruction methods, that finds a solution through a convex optimization: $$f^\star \in \text{argmin}_f \frac{1}{2}\|y-\Phi f\|^2 + \lambda J(f)$$

Where $J(f)$ is a prior energy. In this tour we consider a a Sobolev prior (the image is uniformly smooth) and a sparse wavelet prior (the image is expected to be compressible in a wavelet basis).

Note that the parameter $\lambda$ should be carefully chosen to fit the noise level.

Tomography Acquisition over the Fourier Domain¶

The tomography operator computes the projection of the signal along rays. It corresponds to the Radon transform $$ (T_\theta f)( t ) = \int_{\Delta_{\theta,t}} f(x) d x, $$ where $\Delta_{\theta,t}$ is the line $$ \Delta_{\theta,t} = { x = (x_1,x_2) ::: x_1 \cos(\theta)

• x_2 \sin(\theta) = t }. $$

A tomography measurment computes $$ \{ T_{\theta_i} f \}_{i=0}^{q-1}$$ for a small sub-set of $q$ orientation $ \theta_i \in [0,\pi) $.

Set the number $q$ of rays used for the experiments.

Size of the image, the number of pixels is $N=n^2$.

The Fourier Slice Theorem shows that obtaining the measurements $ T_{\theta} f $ is equivalent to computing slice of the 2D Fourier transform of $f$ along a ray of direction $\theta$ $$ \hat T_{\theta} f(\omega) = \hat f( r \cos(\theta), r \cos(\theta) ). $$

In this tour, we thus consider a discrete Fourier version of the Tomography operator. It corresponds to sampling the Fourier transform of the image along discretized rays $\Omega$.

Compute the set of rays $ \Omega $ by computing the mask $ \xi $ so that $ \xi(\omega)=1 $ if $ \omega \in \Omega $ and $ \xi(\omega) = 0 $ otherwise.

In Matlab, the 0 frequency is in the upper-left corner.

Display the mask.

Sampling the Fourier transform for points is equivalent to masking the Fourier transform $\hat f$ of $f$. We thus define our Fourier tomography measurements operator as $$ \Phi f(\omega) = \hat f(\omega) \xi(\omega) $$ where $ \hat f $ is the 2-D discrete Fourier transform of $f$.

Define a shortcut for the operator. Note the normalization of the FFT by $1/n$ to make it orthonormal.

Load the clean high resolution image.

We consider here noisy tomographic measurements. $$ y = \Phi f_0 + w $$ where $w$ is a Gaussian white noise of variance $\sigma^2$.

Variance of the noise.

Measurements.

The pseudo inverse $\Phi^+$ is equal to the transposed operator $\Phi^*$. It corresponds to inverting the Fourier transform after multiplication by the mask. $$ \Phi^* y = \Phi^+ y = \mathcal{F}^{-1}( y \xi ), $$ where $ \mathcal{F}^{-1}(\hat f) = f $ is the inverse Fourier transform.

Define a shortcut for the transpose operator.

Exercise 1

Compute and display the pseudo inverse reconstruction $ \Phi^+ y $. What do you observe ?

Sobolev Regularization¶

To remove some noise while inverting the operator, we can penalize high frequencies using Sobolev regularization.

The Sobolev prior reads: $$J(f) = \frac{1}{2} \sum_x \|\nabla f(x)\|^2 = \frac{1}{2}\sum_{\omega} S(\omega) \|\hat f(\omega)\|^2 $$ where $S(\omega)=\|\omega\|^2$.

Compute the Sobolev prior penalty $S$ (rescale to $[0,1]$).

We need to compute the solution of: $$f^\star \in \text{argmin}_f \frac{1}{2}\|y-\Phi f\|^2 + \lambda J(f)$$

Regularization parameter $\lambda$:

Since both the prior $J$ and the operator $\Phi$ can be written over the Fourier domain, one can compute the solution to the inversion with Sobolev prior simply with the Fourier coefficients: $$\hat f^\star(\omega) = \frac{\hat y(\omega) \xi(\omega)}{ \xi(\omega) + \lambda S(\omega) }$$

Perform the inversion.

Display the result.

Exercise 2

Find the optimal solution |fSob| by testing several value of $\lambda$.

Display optimal result.

Soft Thresholding in a Basis¶

The soft thresholding operator is at the heart of $\ell^1$ minimization schemes. It can be applied to coefficients $a$, or to an image $f$ in an ortho-basis.

The soft thresholding is a 1-D functional that shrinks the value of coefficients. $$ s_T(u)=\max(0,1-T/|u|)u $$

Define a shortcut for this soft thresholding 1-D functional.

Display a curve of the 1D soft thresholding.

Note that the function |SoftThresh| can also be applied to vector (because of Matlab/Scilab vectorialized computation), which defines an operator on coefficients: $$ S_T(a) = ( s_T(a_m) )_m. $$

In the next section, we use an orthogonal wavelet basis $\Psi$.

We set the parameters of the wavelet transform.

Shortcut for $\Psi$ and $\Psi^*$ in the orthogonal case.

The soft thresholding opterator in the basis $\Psi$ is defined as $$S_T^\Psi(f) = \sum_m s_T( \langle f,\psi_m \rangle ) \psi_m $$

It thus corresponds to applying the transform $\Psi^*$, thresholding the coefficients using $S_T$ and then undoing the transform using $\Psi$. $$ S_T^\Psi(f) = \Psi \circ S_T \circ \Psi^*$$

This soft thresholding corresponds to a denoising operator.