<h2><center> Tomography Reconstruction using CIL </center></h2>

### In this demo, we reconstruct simulated tomographic data using:

* the Filtered Back Projection (FBP) algorithm, 
* Total variation (TV) regularisation under a non-negativity constraint

    <a id='tv_reg'></a>
    $$\begin{equation}
    \underset{u}{\operatorname{argmin}} \frac{1}{2} \| A u - g\|^{2}_{2} + \alpha\,\mathrm{TV}(u) + \mathbb{I}_{\{u\geq0\}}(u)
    \tag{1}
    \end{equation}$$

    where,

    * g is a noisy data corrupted with Gaussian noise and
    * $A$ is the Projection operator.


```python
# Import libraries

from cil.framework import  AcquisitionGeometry

from cil.optimisation.functions import L2NormSquared, BlockFunction, MixedL21Norm, IndicatorBox
from cil.optimisation.operators import GradientOperator, BlockOperator
from cil.optimisation.algorithms import PDHG
from cil.plugins.astra.operators import ProjectionOperator
from cil.plugins.astra.processors import FBP
from cil.plugins import TomoPhantom
from cil.utilities.display import show2D, show_geometry
from cil.utilities import noise

import matplotlib.pyplot as plt

import numpy as np
```

We first configure our acquisition geometry, e.g., **2D parallel geometry**. Then, the image geometry is extracted and used to configure our phantom. To create our simulated phantoms, we use the [Tomophantom](https://github.com/dkazanc/TomoPhantom) library.


```python
# Detectors
N = 256
detectors =  N

# Angles
angles = np.linspace(0,180,180, dtype='float32')

# Setup acquisition geometry
ag = AcquisitionGeometry.create_Parallel2D()\
                        .set_angles(angles)\
                        .set_panel(detectors, pixel_size=0.1)
# Get image geometry
ig = ag.get_ImageGeometry()

# Get phantom
phantom = TomoPhantom.get_ImageData(12, ig)
```

Next, we create our simulated tomographic data by projecting our noiseless `phantom` to the acquisition space. Using the image geometry `ig` and acquisition geometry `ag`, we define the `ProjectionOperator` with `device=cpu` or `device=gpu`. Finally,  Gaussian noise is added to the noiseless sinogram.


```python
# Create projection operator using Astra-Toolbox. Available CPU/CPU
A = ProjectionOperator(ig, ag, device = 'cpu')

# Create an acqusition data (numerically)
sino = A.direct(phantom)

# Simulate Gaussian noise for the sinogram
gaussian_var = 0.5
gaussian_mean = 0

n1 = np.random.normal(gaussian_mean, gaussian_var, size = ag.shape)
                      
noisy_sino = ag.allocate()
noisy_sino.fill(n1 + sino.array)
noisy_sino.array[noisy_sino.array<0]=0

# Show numerical and noisy sinograms
show2D([phantom, sino, noisy_sino], title = ['Ground Truth','Sinogram','Noisy Sinogram'], num_cols=3, cmap = 'inferno')

```

For our first reconstruction, we use the Filtered BackProjection algorithm, i.e., `FBP` applied to our noisy sinogram `noisy_sino`.


```python
# Setup and run the FBP algorithm
fbp_recon = FBP(ig, ag,  device = 'cpu')(noisy_sino)

# Show reconstructions
show2D([phantom, fbp_recon], 
       title = ['Ground Truth','FBP reconstruction'], 
       cmap = 'inferno', fix_range=(0,1.), size=(10,10))
```

In the above reconstruction noise is not penalised. In order to remove noise artifacts, we will use the TotalVariation regularisation as shown in the minimisation problem above.

We solve [(1)](#tv_reg), using the **Primal-Dual Hybrid Gradient (PDHG)** algorithm introduced in [ChambollePock](https://link.springer.com/article/10.1007/s10851-010-0251-1). We need to write [(1)](#tgv_reg) in the following general form $$\underset{x\in \mathbb{X} }{\operatorname{argmin}} f(Kx) + g(x).$$

We let $x=u$ and

* define an operator $K:\mathbb{X}\rightarrow\mathbb{Y}$ as$\\[10pt]$

    $$\begin{equation}
    K = 
    \begin{bmatrix}
    A \\
    D 
    \end{bmatrix} \quad\Rightarrow\quad
    Kx = Ku = 
    \begin{bmatrix}
    \mathcal{A}u\\
    Du\\
    \end{bmatrix} =
    \begin{bmatrix}
    y_{1}\\
    y_{2}
    \end{bmatrix} = y\in \mathbb{Y},
    \label{def_K}
    \end{equation}$$ 
    
* define a function $f:\mathbb{Y}\rightarrow\mathbb{R}$ as$\\[10pt]$

    $$\begin{equation}
    \begin{aligned}
    & f(y)  := f(y_{1}, y_{2}) = f_{1}(y_{1}) +  f_{2}(y_{2}) , \mbox{ where},\\[10pt]
    & f_{1}(y_{1}) :=  \frac{1}{2}\| y_{1} - g\|_{2}^{2},\, f_{2}(y_{2}) :=  \alpha \|y_{2}\|_{2,1}
    \end{aligned}
    \label{def_f}
    \end{equation}\\[10pt]$$
    
* and define a function $g(x) = g(u) = \mathbb{I}_{\{u\geq0\}}(u)$ for the non-negativity constraint.




```python
# Define BlockOperator K
Grad = GradientOperator(ig)
K = BlockOperator(A, Grad)

# Define BlockFunction f
alpha = 0.1
f1 = 0.5 * L2NormSquared(b=noisy_sino)
f2 = alpha * MixedL21Norm()
f = BlockFunction(f1, f2)

# Define Function g
g = IndicatorBox(lower=0)

# Primal/Dual stepsizes
normK = K.norm()
sigma = 1./normK
tau = 1./normK

# Setup and run PDHG
pdhg = PDHG(f = f, g = g, operator = K, sigma = sigma, tau = tau, 
            max_iteration = 200,
            update_objective_interval = 50)
pdhg.run(verbose=2)
```

Finally, we compare the PDHG and FBP reconstructions and plot the middle line profiles.


```python
show2D([pdhg.solution,fbp_recon, phantom], title = ['TV regularisation','FBP','Ground Truth'], cmap = 'inferno', num_cols=3, fix_range=(0,1.))

# Plot middle line profile
plt.figure(figsize=(30,15))
plt.plot(phantom.subset(horizontal_y = int(N/2)).as_array(), label = 'Ground Truth', linewidth=5)
plt.plot(fbp_recon.subset(horizontal_y = int(N/2)).as_array(), label = 'FBP', linewidth=5, linestyle='dashed')
plt.plot(pdhg.solution.subset(horizontal_y = int(N/2)).as_array(), label = 'TV', linewidth=5)
plt.legend()
plt.title('Middle Line Profiles')
plt.show()
```


```python

```

350825df28b1acf1766b50f679fd842cdb474886
ipynb
Jupyter Notebook
binder/TomographyReconstruction.ipynb
Sophilyplum/CIL-Demos
f1a2d4b7675540f5e56accb255a2ea66a1931079
Apache-2.0
2019-11-13T05:20:46.000Z
2021-03-18T14:28:59.000Z
binder/TomographyReconstruction.ipynb
Sophilyplum/CIL-Demos
f1a2d4b7675540f5e56accb255a2ea66a1931079
Apache-2.0
2019-06-25T20:44:51.000Z
2021-05-07T12:25:32.000Z
binder/TomographyReconstruction.ipynb
Sophilyplum/CIL-Demos
f1a2d4b7675540f5e56accb255a2ea66a1931079
Apache-2.0
2021-06-18T10:40:57.000Z
2021-11-09T10:32:22.000Z
Qwen/Qwen-72B
1. YES
2. YES