matlab/inverse_9_tomography_sobsparse/index.html


<!DOCTYPE html
  PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN">
<html xmlns:mwsh="http://www.mathworks.com/namespace/mcode/v1/syntaxhighlight.dtd">
   <head>
      <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
   
      <!--
This HTML is auto-generated from an M-file.
To make changes, update the M-file and republish this document.
      --><script type="text/javascript" src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS_HTML"></script><p style="font-size:0px">
         \[
         \newcommand{\NN}{\mathbb{N}}
         \newcommand{\CC}{\mathbb{C}}
         \newcommand{\GG}{\mathbb{G}}
         \newcommand{\LL}{\mathbb{L}}
         \newcommand{\PP}{\mathbb{P}}
         \newcommand{\QQ}{\mathbb{Q}}
         \newcommand{\RR}{\mathbb{R}}
         \newcommand{\VV}{\mathbb{V}}
         \newcommand{\ZZ}{\mathbb{Z}}
         \newcommand{\FF}{\mathbb{F}}
         \newcommand{\KK}{\mathbb{K}}
         \newcommand{\UU}{\mathbb{U}}
         \newcommand{\EE}{\mathbb{E}}
         
         \newcommand{\Aa}{\mathcal{A}}
         \newcommand{\Bb}{\mathcal{B}}
         \newcommand{\Cc}{\mathcal{C}}
         \newcommand{\Dd}{\mathcal{D}}
         \newcommand{\Ee}{\mathcal{E}}
         \newcommand{\Ff}{\mathcal{F}}
         \newcommand{\Gg}{\mathcal{G}}
         \newcommand{\Hh}{\mathcal{H}}
         \newcommand{\Ii}{\mathcal{I}}
         \newcommand{\Jj}{\mathcal{J}}
         \newcommand{\Kk}{\mathcal{K}}
         \newcommand{\Ll}{\mathcal{L}}
         \newcommand{\Mm}{\mathcal{M}}
         \newcommand{\Nn}{\mathcal{N}}
         \newcommand{\Oo}{\mathcal{O}}
         \newcommand{\Pp}{\mathcal{P}}
         \newcommand{\Qq}{\mathcal{Q}}
         \newcommand{\Rr}{\mathcal{R}}
         \newcommand{\Ss}{\mathcal{S}}
         \newcommand{\Tt}{\mathcal{T}}
         \newcommand{\Uu}{\mathcal{U}}
         \newcommand{\Vv}{\mathcal{V}}
         \newcommand{\Ww}{\mathcal{W}}
         \newcommand{\Xx}{\mathcal{X}}
         \newcommand{\Yy}{\mathcal{Y}}
         \newcommand{\Zz}{\mathcal{Z}}
         
         \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}
         \renewcommand{\phi}{\varphi}
         \renewcommand{\th}{\theta}
         \newcommand{\om}{\omega}
         \newcommand{\Om}{\Omega}
         \renewcommand{\epsilon}{\varepsilon}
         
         \newcommand{\Calpha}{\mathrm{C}^\al}
         \newcommand{\Cbeta}{\mathrm{C}^\be}
         \newcommand{\Cal}{\text{C}^\al}
         \newcommand{\Cdeux}{\text{C}^{2}}
         \newcommand{\Cun}{\text{C}^{1}}
         \newcommand{\Calt}[1]{\text{C}^{#1}}
         
         \newcommand{\lun}{\ell^1}
         \newcommand{\ldeux}{\ell^2}
         \newcommand{\linf}{\ell^\infty}
         \newcommand{\ldeuxj}{{\ldeux_j}}
         \newcommand{\Lun}{\text{\upshape L}^1}
         \newcommand{\Ldeux}{\text{\upshape L}^2}
         \newcommand{\Lp}{\text{\upshape L}^p}
         \newcommand{\Lq}{\text{\upshape L}^q}
         \newcommand{\Linf}{\text{\upshape L}^\infty}
         \newcommand{\lzero}{\ell^0}
         \newcommand{\lp}{\ell^p}
         
         
         \renewcommand{\d}{\ins{d}}
         
         \newcommand{\Grad}{\text{Grad}}
         \newcommand{\grad}{\text{grad}}
         \renewcommand{\div}{\text{div}}
         \newcommand{\diag}{\text{diag}}
         
         \newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }
         \newcommand{\pdd}[2]{ \frac{ \partial^2 #1}{\partial #2^2} }
         
         \newcommand{\dotp}[2]{\langle #1,\,#2\rangle}
         \newcommand{\norm}[1]{|\!| #1 |\!|}
         \newcommand{\normi}[1]{\norm{#1}_{\infty}}
         \newcommand{\normu}[1]{\norm{#1}_{1}}
         \newcommand{\normz}[1]{\norm{#1}_{0}}
         \newcommand{\abs}[1]{\vert #1 \vert}
         
         
         \newcommand{\argmin}{\text{argmin}}
         \newcommand{\argmax}{\text{argmax}}
         \newcommand{\uargmin}[1]{\underset{#1}{\argmin}\;}
         \newcommand{\uargmax}[1]{\underset{#1}{\argmax}\;}
         \newcommand{\umin}[1]{\underset{#1}{\min}\;}
         \newcommand{\umax}[1]{\underset{#1}{\max}\;}
         
         \newcommand{\pa}[1]{\left( #1 \right)}
         \newcommand{\choice}[1]{ \left\{  \begin{array}{l} #1 \end{array} \right. }
         
         \newcommand{\enscond}[2]{ \left\{ #1 \;:\; #2 \right\} }
         
         \newcommand{\qandq}{ \quad \text{and} \quad }
         \newcommand{\qqandqq}{ \qquad \text{and} \qquad }
         \newcommand{\qifq}{ \quad \text{if} \quad }
         \newcommand{\qqifqq}{ \qquad \text{if} \qquad }
         \newcommand{\qwhereq}{ \quad \text{where} \quad }
         \newcommand{\qqwhereqq}{ \qquad \text{where} \qquad }
         \newcommand{\qwithq}{ \quad \text{with} \quad }
         \newcommand{\qqwithqq}{ \qquad \text{with} \qquad }
         \newcommand{\qforq}{ \quad \text{for} \quad }
         \newcommand{\qqforqq}{ \qquad \text{for} \qquad }
         \newcommand{\qqsinceqq}{ \qquad \text{since} \qquad }
         \newcommand{\qsinceq}{ \quad \text{since} \quad }
         \newcommand{\qarrq}{\quad\Longrightarrow\quad}
         \newcommand{\qqarrqq}{\quad\Longrightarrow\quad}
         \newcommand{\qiffq}{\quad\Longleftrightarrow\quad}
         \newcommand{\qqiffqq}{\qquad\Longleftrightarrow\qquad}
         \newcommand{\qsubjq}{ \quad \text{subject to} \quad }
         \newcommand{\qqsubjqq}{ \qquad \text{subject to} \qquad }
         \]
         
      </p>
      <title>Tomography Inversion using Tikhonov and Sparse Regularization</title>
      <NOSCRIPT>
         <DIV STYLE="color:#CC0000; text-align:center"><B>Warning: <A HREF="http://www.math.union.edu/locate/jsMath">jsMath</A> 
               	requires JavaScript to process the mathematics on this page.<BR> 
               	If your browser supports JavaScript, be sure it is enabled.</B></DIV>
         <HR>
      </NOSCRIPT>
      <meta name="generator" content="MATLAB 8.2">
      <meta name="date" content="2014-10-20">
      <meta name="m-file" content="index">
      <LINK REL="stylesheet" HREF="../style.css" TYPE="text/css">
   </head>
   <body>
      <div class="content">
         <h1>Tomography Inversion using Tikhonov and Sparse Regularization</h1>
         <introduction>
            <p>This numerical tour explores the reconstruction from tomographic measurement with Sobolev and sparse regularization.</p>
         </introduction>
         <h2>Contents</h2>
         <div>
            <ul>
               <li><a href="#1">Installing toolboxes and setting up the path.</a></li>
               <li><a href="#8">Tomography and Regularization</a></li>
               <li><a href="#13">Tomography Acquisition over the Fourier Domain</a></li>
               <li><a href="#31">Sobolev Regularization</a></li>
               <li><a href="#41">Soft Thresholding in a Basis</a></li>
               <li><a href="#52">Sparse Regularization</a></li>
               <li><a href="#66">Translation Invariant Regularization</a></li>
               <li><a href="#71">Medical Resonnance Imaging</a></li>
            </ul>
         </div>
         <h2>Installing toolboxes and setting up the path.<a name="1"></a></h2>
         <p>You need to download the following files: <a href="../toolbox_signal.zip">signal toolbox</a> and <a href="../toolbox_general.zip">general toolbox</a>.
         </p>
         <p>You need to unzip these toolboxes in your working directory, so that you have <tt>toolbox_signal</tt> and <tt>toolbox_general</tt> in your directory.
         </p>
         <p><b>For Scilab user:</b> you must replace the Matlab comment '%' by its Scilab counterpart '//'.
         </p>
         <p><b>Recommandation:</b> You should create a text file named for instance <tt>numericaltour.sce</tt> (in Scilab) or <tt>numericaltour.m</tt> (in Matlab) to write all the Scilab/Matlab command you want to execute. Then, simply run <tt>exec('numericaltour.sce');</tt> (in Scilab) or <tt>numericaltour;</tt> (in Matlab) to run the commands.
         </p>
         <p>Execute this line only if you are using Matlab.</p><pre class="codeinput">getd = @(p)path(p,path); <span class="comment">% scilab users must *not* execute this</span>
</pre><p>Then you can add the toolboxes to the path.</p><pre class="codeinput">getd(<span class="string">'toolbox_signal/'</span>);
getd(<span class="string">'toolbox_general/'</span>);
</pre><h2>Tomography and Regularization<a name="8"></a></h2>
         <p>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.
         </p>
         <p>The tomography operator \(\Phi\) is ill-posed, so it cannot be inverted.</p>
         <p>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)\]
         </p>
         <p>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).
         </p>
         <p>Note that the parameter \(\lambda\) should be carefully chosen to fit the noise level.</p>
         <h2>Tomography Acquisition over the Fourier Domain<a name="13"></a></h2>
         <p>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 \}. \]
         </p>
         <p>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) \).
         </p>
         <p>Set the number \(q\) of rays used for the experiments.</p><pre class="codeinput">q = 32;
</pre><p>Size of the image, the number of pixels is \(N=n^2\).</p><pre class="codeinput">n = 256;
</pre><p>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) ). \]
         </p>
         <p>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\).
         </p>
         <p>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.
         </p><pre class="codeinput">Theta = linspace(0,pi,q+1); Theta(end) = [];
xi = zeros(n,n);
<span class="keyword">for</span> theta = Theta
    t = linspace(-1,1,3*n)*n;
    x = round(t.*cos(theta)) + n/2+1; y = round(t.*sin(theta)) + n/2+1;
    I = find(x&gt;0 &amp; x&lt;=n &amp; y&gt;0 &amp; y&lt;=n); x = x(I); y = y(I);
    xi(x+(y-1)*n) = 1;
<span class="keyword">end</span>
</pre><p>In Matlab, the 0 frequency is in the upper-left corner.</p><pre class="codeinput">xi = fftshift(xi);
</pre><p>Display the mask.</p><pre class="codeinput">clf;
imageplot(fftshift(xi));
</pre><img vspace="5" hspace="5" src="index_01.png"> <p>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\).
         </p>
         <p>Define a shortcut for the operator. Note the normalization of the FFT by \(1/n\) to make it orthonormal.</p><pre class="codeinput">Phi  = @(f)fft2(f).*xi / n;
</pre><p>Load the clean high resolution image.</p><pre class="codeinput">name = <span class="string">'mri'</span>;
f0 = load_image(name,n);
f0 = rescale(f0);
</pre><p>We consider here noisy tomographic measurements. \[ y = \Phi f_0 + w \] where \(w\) is a Gaussian white noise of variance
            \(\sigma^2\).
         </p>
         <p>Variance of the noise.</p><pre class="codeinput">sigma = .2;
</pre><p>Measurements.</p><pre class="codeinput">y = Phi(f0) + sigma*randn(n,n);
</pre><p>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.
         </p>
         <p>Define a shortcut for the transpose operator.</p><pre class="codeinput">PhiS = @(y)real(ifft2(y.*xi))*n;
</pre><p><i>Exercice 1:</i> (<a href="../missing-exo/">check the solution</a>) Compute and display the pseudo inverse reconstruction \( \Phi^+ y \). What do you observe ?
         </p><pre class="codeinput">exo1;
</pre><img vspace="5" hspace="5" src="index_02.png"> <h2>Sobolev Regularization<a name="31"></a></h2>
         <p>To remove some noise while inverting the operator, we can penalize high frequencies using Sobolev regularization.</p>
         <p>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\).
         </p>
         <p>Compute the Sobolev prior penalty \(S\)&nbsp;(rescale to \([0,1]\)).</p><pre class="codeinput">x = [0:n/2-1, -n/2:-1];
[Y,X] = meshgrid(x,x);
S = (X.^2 + Y.^2)*(2/n)^2;
</pre><p>We need to compute the solution of:  \[f^\star \in \text{argmin}_f \frac{1}{2}\|y-\Phi f\|^2 + \lambda J(f)\]</p>
         <p>Regularization parameter \(\lambda\):</p><pre class="codeinput">lambda = 2;
</pre><p>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) }\]
         </p>
         <p>Perform the inversion.</p><pre class="codeinput">fSob = real( ifft2( y .* xi ./ ( xi + lambda*S) ) )*n;
</pre><p>Display the result.</p><pre class="codeinput">clf;
imageplot(clamp(fSob), [<span class="string">'Sobolev inversion, SNR='</span> num2str(snr(f0,fSob),3) <span class="string">'dB'</span>] );
</pre><img vspace="5" hspace="5" src="index_03.png"> <p><i>Exercice 2:</i> (<a href="../missing-exo/">check the solution</a>) Find the optimal solution <tt>fSob</tt> by testing several value of \(\lambda\).
         </p><pre class="codeinput">exo2;
</pre><img vspace="5" hspace="5" src="index_04.png"> <p>Display optimal result.</p><pre class="codeinput">clf;
imageplot(clamp(fSob), [<span class="string">'Sobolev inversion, SNR='</span> num2str(snr(f0,fSob),3) <span class="string">'dB'</span>]);
</pre><img vspace="5" hspace="5" src="index_05.png"> <h2>Soft Thresholding in a Basis<a name="41"></a></h2>
         <p>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.
         </p>
         <p>The soft thresholding is a 1-D functional that shrinks the value of coefficients. \[ s_T(u)=\max(0,1-T/|u|)u \]</p>
         <p>Define a shortcut for this soft thresholding 1-D functional.</p><pre class="codeinput">SoftThresh = @(x,T)x.*max( 0, 1-T./max(abs(x),1e-10) );
</pre><p>Display a curve of the 1D soft thresholding.</p><pre class="codeinput">clf;
T = linspace(-1,1,1000);
plot( T, SoftThresh(T,.5) );
axis(<span class="string">'equal'</span>);
</pre><img vspace="5" hspace="5" src="index_06.png"> <p>Note that the function <tt>SoftThresh</tt> 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. \]
         </p>
         <p>In the next section, we use an orthogonal wavelet basis \(\Psi\).</p>
         <p>We set the parameters of the wavelet transform.</p><pre class="codeinput">Jmax = log2(n)-1;
Jmin = Jmax-3;
</pre><p>Shortcut for \(\Psi\) and \(\Psi^*\) in the orthogonal case.</p><pre class="codeinput">options.ti = 0; <span class="comment">% use orthogonality.</span>
Psi = @(a)perform_wavelet_transf(a, Jmin, -1,options);
PsiS = @(f)perform_wavelet_transf(f, Jmin, +1,options);
</pre><p>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 \]
         </p>
         <p>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^*\]
         </p><pre class="codeinput">SoftThreshPsi = @(f,T)Psi(SoftThresh(PsiS(f),T));
</pre><p>This soft thresholding corresponds to a denoising operator.</p><pre class="codeinput">clf;
imageplot( clamp(SoftThreshPsi(f0,.1)) );
</pre><img vspace="5" hspace="5" src="index_07.png"> <h2>Sparse Regularization<a name="52"></a></h2>
         <p>To better reconstruct edges, we consider a sparsity regularization in an orthogonal basis \(\{ \psi_m \}_m\): \[ J(f)=\sum_m
            \vert \langle f,\psi_m \rangle\vert = \| \Psi^* f \|_1. \]
         </p>
         <p>The inverse problem resolution thus require to solve: \[f^{\star} \in \text{argmin}_f \: E(f) = \frac{1}{2}\|y-\Phi f\|^2
            + \lambda \sum_m \vert \langle f,\psi_m \rangle\vert. \]
         </p>
         <p>Set up the value of the regularization parameter \(\lambda\).</p><pre class="codeinput">lambda = 0.1;
</pre><p>To solve this non-smooth optimization problem, one can use forward-backward splitting, also known as iterative soft thresholding.</p>
         <p>It computes a series of images \(f^{(\ell)}\) defined as \[ f^{(\ell+1)} = S_{\tau\lambda}^{\Psi}( f^{(\ell)} - \tau \Phi^{*}
            (\Phi f^{(\ell)} - y)  ) \]
         </p>
         <p>For \(f^{(\ell)}\) to converge to a solution of the problem, the gradient step size should be chosen as \[\tau &lt; \frac{2}{\|\Phi^*
            \Phi\|}\]
         </p>
         <p>Since \(\Phi\) is an operator of norm 1, this must be smaller than 2.</p><pre class="codeinput">tau = 1.5;
</pre><p>Initialize the solution.</p><pre class="codeinput">fSpars = PhiS(y);
</pre><p>First step: perform one step of gradient descent of the energy \( \|y-\Phi f\|^2 \).</p><pre class="codeinput">fSpars = fSpars + tau * PhiS( y-Phi(fSpars) );
</pre><p>Second step: denoise the solution by thresholding.</p><pre class="codeinput">fSpars = SoftThreshPsi( fSpars, lambda*tau );
</pre><p><i>Exercice 3:</i> (<a href="../missing-exo/">check the solution</a>) Perform the iterative soft thresholding. Monitor the decay of the energy \(E\) you are minimizing.
         </p><pre class="codeinput">exo3;
</pre><img vspace="5" hspace="5" src="index_08.png"> <p>Display the result.</p><pre class="codeinput">clf;
imageplot(clamp(fSpars), [<span class="string">'Sparsity inversion, SNR='</span> num2str(snr(f0,fSpars),3) <span class="string">'dB'</span>]);
</pre><img vspace="5" hspace="5" src="index_09.png"> <p><i>Exercice 4:</i> (<a href="../missing-exo/">check the solution</a>) Try to find the best threshold \(\lambda\). To this end, perform a <b>lot</b> of iterations, and progressively decay the threshold \(\lambda\) during the iterations. Record the best result in <tt>fBestOrtho</tt>.
         </p><pre class="codeinput">exo4;
</pre><img vspace="5" hspace="5" src="index_10.png"> <p>Display the result.</p><pre class="codeinput">clf;
imageplot(clamp(fBestOrtho), [<span class="string">'Sparsity in Orthogonal Wavelets, SNR='</span> num2str(snr(f0,fBestOrtho),3) <span class="string">'dB'</span>]);
</pre><img vspace="5" hspace="5" src="index_11.png"> <h2>Translation Invariant Regularization<a name="66"></a></h2>
         <p>Denoising artefact are reduced by using a translation invariant wavelet frame.</p>
         <p>Define new shortcut for the translation invariant wavelet transforms.</p><pre class="codeinput">options.ti = 1;
Psi = @(a)perform_wavelet_transf(a, Jmin, -1,options);
PsiS = @(f)perform_wavelet_transf(f, Jmin, +1,options);
</pre><p>Re-define the thresholding.</p><pre class="codeinput">SoftThreshPsi = @(f,T)Psi(SoftThresh(PsiS(f),T));
</pre><p><i>Exercice 5:</i> (<a href="../missing-exo/">check the solution</a>) Use the iterative thresholding but this time with the translation invariant wavelet transform. Find the best value of \(\lambda\)
            and record the best result in <tt>fBestTI</tt>.
         </p><pre class="codeinput">exo5;
</pre><img vspace="5" hspace="5" src="index_12.png"> <p>Display the result.</p><pre class="codeinput">clf;
imageplot(clamp(fBestTI), [<span class="string">'Sparsity in TI Wavelets, SNR='</span> num2str(snr(f0,fBestTI),3) <span class="string">'dB'</span>]);
</pre><img vspace="5" hspace="5" src="index_13.png"> <h2>Medical Resonnance Imaging<a name="71"></a></h2>
         <p>MRI acquisition directly gather Fourier measurement along a curve in the Fourier plane.</p>
         <p>Create a spiral for the sampling. The parameter \(\alpha\) controls the density near the origin.</p><pre class="codeinput">alpha = 3;
q = 30;
t = linspace(0,1.5,n*q*10);
x = round(.5*n*t.^alpha.*cos(2*pi*q*t)) + n/2+1;
y = round(.5*n*t.^alpha.*sin(2*pi*q*t)) + n/2+1;
I = find(x&gt;0 &amp; x&lt;=n &amp; y&gt;0 &amp; y&lt;=n); x = x(I); y = y(I);
xi = zeros(n,n); xi(x+(y-1)*n) = 1;
</pre><p>Display the sampling pattern.</p><pre class="codeinput">clf;
imageplot(xi);
</pre><img vspace="5" hspace="5" src="index_14.png"> <p><i>Exercice 6:</i> (<a href="../missing-exo/">check the solution</a>) Compare Sobolev and Sparse reconstruction for MRI imaging. For a given number of Fourier sample, compare the quality of
            the reconstruction for different \(\alpha\).
         </p><pre class="codeinput">exo6;
</pre><p class="footer"><br>
            Copyright  (c) 2010 Gabriel Peyre<br></p>
      </div>
      <!--
##### SOURCE BEGIN #####
%% Tomography Inversion using Tikhonov and Sparse Regularization
% This numerical tour explores the reconstruction from tomographic
% measurement with Sobolev and sparse regularization.

%% Installing toolboxes and setting up the path.

%%
% You need to download the following files: 
% <../toolbox_signal.zip signal toolbox> and 
% <../toolbox_general.zip general toolbox>.

%%
% You need to unzip these toolboxes in your working directory, so
% that you have 
% |toolbox_signal| and 
% |toolbox_general|
% in your directory.

%%
% *For Scilab user:* you must replace the Matlab comment '%' by its Scilab
% counterpart '//'.

%%
% *Recommandation:* You should create a text file named for instance |numericaltour.sce| (in Scilab) or |numericaltour.m| (in Matlab) to write all the
% Scilab/Matlab command you want to execute. Then, simply run |exec('numericaltour.sce');| (in Scilab) or |numericaltour;| (in Matlab) to run the commands. 

%%
% Execute this line only if you are using Matlab.

getd = @(p)path(p,path); % scilab users must *not* execute this

%%
% Then you can add the toolboxes to the path.

getd('toolbox_signal/');
getd('toolbox_general/');


%% 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.

q = 32;

%%
% Size of the image, the number of pixels is \(N=n^2\).

n = 256;

%%
% 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.

Theta = linspace(0,pi,q+1); Theta(end) = [];
xi = zeros(n,n);
for theta = Theta
    t = linspace(-1,1,3*n)*n;
    x = round(t.*cos(theta)) + n/2+1; y = round(t.*sin(theta)) + n/2+1;
    I = find(x>0 & x<=n & y>0 & y<=n); x = x(I); y = y(I);
    xi(x+(y-1)*n) = 1;
end

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

xi = fftshift(xi);

%%
% Display the mask.

clf;
imageplot(fftshift(xi));

%%
% 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.

Phi  = @(f)fft2(f).*xi / n;

%%
% Load the clean high resolution image.

name = 'mri';
f0 = load_image(name,n);
f0 = rescale(f0);


%%
% 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.

sigma = .2;

%%
% Measurements.

y = Phi(f0) + sigma*randn(n,n);

%%
% 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.

PhiS = @(y)real(ifft2(y.*xi))*n;

%%
% _Exercice 1:_ (<../missing-exo/ check the solution>)
% Compute and display the pseudo inverse reconstruction
% \( \Phi^+ y \). What do you observe ?

exo1;


%% 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]\)).

x = [0:n/2-1, -n/2:-1];
[Y,X] = meshgrid(x,x);
S = (X.^2 + Y.^2)*(2/n)^2;

%%
% 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\):

lambda = 2;

%%
% 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.

fSob = real( ifft2( y .* xi ./ ( xi + lambda*S) ) )*n;

%%
% Display the result.

clf;
imageplot(clamp(fSob), ['Sobolev inversion, SNR=' num2str(snr(f0,fSob),3) 'dB'] );

%%
% _Exercice 2:_ (<../missing-exo/ check the solution>)
% Find the optimal solution |fSob| by testing several value of
% \(\lambda\).

exo2;

%%
% Display optimal result.

clf;
imageplot(clamp(fSob), ['Sobolev inversion, SNR=' num2str(snr(f0,fSob),3) 'dB']);


%% 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.

SoftThresh = @(x,T)x.*max( 0, 1-T./max(abs(x),1e-10) );

%%
% Display a curve of the 1D soft thresholding.

clf;
T = linspace(-1,1,1000);
plot( T, SoftThresh(T,.5) );
axis('equal');

%%
% 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. 

Jmax = log2(n)-1;
Jmin = Jmax-3;

%%
% Shortcut for \(\Psi\) and \(\Psi^*\) in the orthogonal case.

options.ti = 0; % use orthogonality.
Psi = @(a)perform_wavelet_transf(a, Jmin, -1,options);
PsiS = @(f)perform_wavelet_transf(f, Jmin, +1,options);

%% 
% 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^*\]

SoftThreshPsi = @(f,T)Psi(SoftThresh(PsiS(f),T));

%%
% This soft thresholding corresponds to a denoising operator.

clf;
imageplot( clamp(SoftThreshPsi(f0,.1)) );


%% Sparse Regularization
% To better reconstruct edges, we consider a sparsity regularization in an
% orthogonal basis \(\{ \psi_m \}_m\):
% \[ J(f)=\sum_m \vert \langle f,\psi_m \rangle\vert = \| \Psi^* f \|_1. \]

%%
% The inverse problem resolution thus require to solve:
% \[f^{\star} \in \text{argmin}_f \: E(f) = \frac{1}{2}\|y-\Phi f\|^2 + \lambda \sum_m \vert \langle f,\psi_m \rangle\vert. \]

%% 
% Set up the value of the regularization parameter \(\lambda\).

lambda = 0.1;

%% 
% To solve this non-smooth optimization problem, one can use
% forward-backward splitting, also known as iterative soft thresholding.

%%
% It computes a series of images \(f^{(\ell)}\) defined as
% \[ f^{(\ell+1)} = S_{\tau\lambda}^{\Psi}( f^{(\ell)} - \tau \Phi^{*}
% (\Phi f^{(\ell)} - y)  ) \]

%%
% For \(f^{(\ell)}\) to converge to a solution of the problem, the gradient
% step size should be chosen as
% \[\tau < \frac{2}{\|\Phi^* \Phi\|}\]

%%
% Since \(\Phi\) is an operator of norm 1, this must be smaller than
% 2.

tau = 1.5;

%%
% Initialize the solution.

fSpars = PhiS(y);

%%
% First step: perform one step of gradient descent of the energy
% \( \|y-\Phi f\|^2 \).

fSpars = fSpars + tau * PhiS( y-Phi(fSpars) );

%%
% Second step: denoise the solution by thresholding.

fSpars = SoftThreshPsi( fSpars, lambda*tau );

%%
% _Exercice 3:_ (<../missing-exo/ check the solution>)
% Perform the iterative soft thresholding.
% Monitor the decay of the energy \(E\) you are minimizing.

exo3;

%% 
% Display the result.

clf;
imageplot(clamp(fSpars), ['Sparsity inversion, SNR=' num2str(snr(f0,fSpars),3) 'dB']);

%%
% _Exercice 4:_ (<../missing-exo/ check the solution>)
% Try to find the best threshold \(\lambda\). To this end, perform a *lot*
% of iterations, and progressively decay the threshold \(\lambda\) during the
% iterations. Record the best result in |fBestOrtho|.

exo4;

%% 
% Display the result.

clf;
imageplot(clamp(fBestOrtho), ['Sparsity in Orthogonal Wavelets, SNR=' num2str(snr(f0,fBestOrtho),3) 'dB']);


%% Translation Invariant Regularization
% Denoising artefact are reduced by using a translation invariant wavelet
% frame.

%%
% Define new shortcut for the translation invariant wavelet transforms.

options.ti = 1;
Psi = @(a)perform_wavelet_transf(a, Jmin, -1,options);
PsiS = @(f)perform_wavelet_transf(f, Jmin, +1,options);

%%
% Re-define the thresholding.

SoftThreshPsi = @(f,T)Psi(SoftThresh(PsiS(f),T));

%%
% _Exercice 5:_ (<../missing-exo/ check the solution>)
% Use the iterative thresholding but this time
% with the translation invariant wavelet transform.
% Find the best value of \(\lambda\) and 
% record the best result in |fBestTI|.

exo5;


%% 
% Display the result.

clf;
imageplot(clamp(fBestTI), ['Sparsity in TI Wavelets, SNR=' num2str(snr(f0,fBestTI),3) 'dB']);


%% Medical Resonnance Imaging 
% MRI acquisition directly gather Fourier measurement along a curve in the
% Fourier plane.

%%
% Create a spiral for the sampling.
% The parameter \(\alpha\) controls the density near the origin.

alpha = 3;
q = 30;
t = linspace(0,1.5,n*q*10);
x = round(.5*n*t.^alpha.*cos(2*pi*q*t)) + n/2+1; 
y = round(.5*n*t.^alpha.*sin(2*pi*q*t)) + n/2+1;
I = find(x>0 & x<=n & y>0 & y<=n); x = x(I); y = y(I);
xi = zeros(n,n); xi(x+(y-1)*n) = 1;

%%
% Display the sampling pattern.

clf;
imageplot(xi);

%%
% _Exercice 6:_ (<../missing-exo/ check the solution>)
% Compare Sobolev and Sparse reconstruction for MRI imaging.
% For a given number of Fourier sample, compare the quality of the 
% reconstruction for different \(\alpha\).

exo6;


##### SOURCE END #####
-->
   </body>
</html>