<jupyter_start><jupyter_text>Primal Dual Hybrid Gradient Algorithm: Tomography In this demo, we learn how to use the Primal Dual hubrid algorithm introduced by [<NAME>](https://hal.archives-ouvertes.fr/hal-00490826/document). We focus on Tomography Reconstruction under an edge-preserving prior, i.e., the __Total variation regularisation__, see [ROF](https://en.wikipedia.org/wiki/Total_variation_denoising). The minimisation problem for Total variation Tomography reconstuction is $$ \underset{u}{\operatorname{argmin}} \frac{1}{2} \| \mathcal{A} u - g\|^{2} + \alpha\,\mathrm{TV}(u) +  \mathbb{I}_{\{u>0\}}(u) $$where,1. The total variation is defined as $$\mathrm{TV}(u) = \|\nabla u \|_{2,1} = \sum \sqrt{ (\partial_{y}u)^{2} + (\partial_{x}u)^{2} }$$1. g is the Acqusisition data obtained from the detector.1. $\mathcal{A}$ is the projection operator ( _Radon transform_ ) that maps from an image-space to an acquisition space, i.e., $\mathcal{A} : X \rightarrow Y, $ where X is an __ImageGeometry__ and Y is an __AcquisitionGeometry__.1. $\alpha$: regularising parameter that measures a trade-off between the fidelity and the regulariser terms.1. $\mathbb{I}_{\{u>0\}}(u) : = \begin{cases}0, & \mbox{ if } u>0\\\infty , & \mbox{ otherwise}\quad\end{cases}$, $\quad$ a positivity constraint for the minimiser $u$. Learning objectives1. Brief intro for non-smooth minimisation problem using PDHG algorithm.1. Setup and run PDHG with (__non-smooth__) $L^{1}$ norm regulariser. __(No BlockFramework)__1. Use __BlockFunction__ and setup PDHG within the __Block Framework__.1. Run Total variation reconstruction with different regularising parameters and compared with FBP and SIRT reconstructions.At the end of this demo, we will be able to reproduce all the reconstructions presented in the figure below. One can observe that the __Tikhonov regularisation__ with $L = \nabla$ was able to remove the noise but could not preserve the edges. However, this can be achieved with the the total variation reconstruction. Prerequisites- AcquisitionData, AcquisitionGeometry, AstraProjectorSimple.- BlockOperator, Gradient.- FBP, SIRT, CGLS, Tikhonov.In order to use the PDHG algorithm for the problem above, we need to express our minimisation problem into the following form:$$\min_{u} \mathcal{F}(K u) + \mathcal{G}(u)$$where we assume that:1. $\mathcal{F}$, $\mathcal{G}$ are __convex__ functionals        - $\mathcal{F}: Y \rightarrow \mathbb{R}$         - $\mathcal{G}: X \rightarrow \mathbb{R}$        2. $K$ is a continuous linear operator acting from a space X to another space Y :$$K : X \rightarrow Y \quad $$ with operator norm  defined as $$\| K \| = \max\{ \|K x\|_{Y} : \|x\|_{X}\leq 1 \}.$$  **Note**: The Gradient operator  has  $\|\nabla\| = \sqrt{8} $ and for the projection operator we use the [Power Method](https://en.wikipedia.org/wiki/Power_iteration) to approximate the greatest eigenvalue of $K$.<jupyter_code># Imports
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

from ccpi.framework import ImageData, TestData, ImageGeometry, AcquisitionGeometry, AcquisitionData, BlockDataContainer

from ccpi.optimisation.functions import L2NormSquared, ZeroFunction, L1Norm, BlockFunction, MixedL21Norm, IndicatorBox, FunctionOperatorComposition
from ccpi.optimisation.operators import Gradient, BlockOperator
from ccpi.optimisation.algorithms import PDHG, SIRT, CGLS

from ccpi.astra.operators import AstraProjectorSimple, AstraProjector3DSimple
from ccpi.astra.processors import FBP

import tomophantom
from tomophantom import TomoP2D
import os, sys

import matplotlib.pyplot as plt


import numpy as np

from utilities import islicer, link_islicer, psnr, plotter2D
from utilities.show_utilities import show

%matplotlib inline<jupyter_output><empty_output><jupyter_text>Why PDHG?In the previous demo, we presented the __Tikhonov regularisation__ for tomography reconstruction, i.e.,$$ \underset{u}{\operatorname{argmin}} \|\mathcal{A} u - g\|^{2}_{2} + \alpha^{2}\|L u\|^{2}_{2} $$where we set $L = \nabla $ or $L = \mathbb{I}$. Due to the $\|\cdot\|^{2}_{2}$ terms, one can observe that the above objective function is differentiable. However, this is not always the case. Consider for example an $L^{1}$ norm for the fidelity, i.e., $\|\mathcal{A} u - g\|_{1}$ or an $L^{1}$ norm of the regulariser i.e., $\|u\|_{1}$.Using the __Primal-Dual hybrid gradient algorithm__, we can solve minimisation problems where the objective is not differentiable, and only a convexity is required. The algorithm is described in the [Appendix](Appendix) and for every iteration, we solve two subproblems, i.e., __primal & dual problems__ where the __[proximal operators](Proximal)__ have a closed form solution or can be solved efficiently using an iterative solver.  L$^{1}$ regularisationLet $L=\mathbb{I}$ in [Tikhonov regularisation](Tikhonov) and replace the$$\alpha^{2}\|L u\|^{2}_{2}\mbox{  with  } \alpha\|u\|_{1}, $$ which results to a non-differentiable objective function. Hence, we have $$ \underset{u}{\operatorname{argmin}} \|\mathcal{A} u - g\|^{2}_{2} + \alpha\|u\|_{1}. $$ Now, we would like to setup and run PDHG algorithm for the above minimisation problem. First, we load a phantom, from the [Tomophantom](https://github.com/dkazanc/TomoPhantom) package. We canchoose different 2D,3D & 4D phantoms from this [library](https://github.com/dkazanc/TomoPhantom/tree/master/PhantomLibrary/models). For the well-known Shepp-Logan phantom, use **_model = 1_** in the code below.<jupyter_code>model = 12 # select a model number from the library
N = 256 # set dimension of the phantom
path = os.path.dirname(tomophantom.__file__)
path_library2D = os.path.join(path, "Phantom2DLibrary.dat")

phantom = TomoP2D.Model(model, N, path_library2D) 

# Define image geometry.
ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, 
                   voxel_size_x = 0.1,
                   voxel_size_y = 0.1)
im_data = ig.allocate()
im_data.fill(phantom)

show(im_data, title = 'TomoPhantom', cmap = 'inferno')<jupyter_output><empty_output><jupyter_text>The next step is to create an __AcquisitionGeometry__ and __AcquisitionData__ using a __Projection Operator__ $\mathcal{A}$.<jupyter_code># Create AcquisitionGeometry and AcquisitionData 
detectors = N
angles = np.linspace(0, np.pi, 180, dtype='float32')
ag = AcquisitionGeometry('parallel','2D', angles, detectors,
                        pixel_size_h = 0.1)

# Create projection operator using Astra-Toolbox. Available CPU/CPU
A = AstraProjectorSimple(ig, ag, device = 'gpu')

# Create an acqusition data (numerically)
sino_num = A.direct(im_data)

# Show numerical sinogram
show(sino_num, title = 'Numerical Sinogram', labels = ['Detectors','Angle'], cmap = 'inferno')

# Simulate Gaussian noise for the sinogram
gaussian_var = 0.2
gaussian_mean = 0 # or np.mean(sino_num.as_array())
n1 = np.random.normal(gaussian_mean, gaussian_var, size = ag.shape)
                      
sino_noisy = ag.allocate()
sino_noisy.fill(n1 + sino_num.as_array())
show(sino_noisy, title = 'Noisy Sinogram', labels = ['Detectors','Angle'], cmap = 'inferno')<jupyter_output><empty_output><jupyter_text>FBP and SIRT reconstuctionsLet's start with a Filtered Back Projection reconstruction and Simultaneous Iterative Reconstruction Technique.Recall, for FBP we type` fbp = FBP(ig, ag, filter_type = 'ram-lak', device = 'gpu') `` fbp.set_input(sino_noisy)`` fbp_recon = fbp.get_output()` and the available filters from Astra are:`‘ram-lak’, ‘shepp-logan’, ‘cosine’, ‘hamming’, ‘hann’, ‘none’, ‘tukey’, ‘lanczos’, ‘triangular’, ‘gaussian’, ‘barlett-hann’, ‘blackman’, ‘nuttall’, ‘blackman-harris’, ‘blackman-nuttall’, ‘flat-top’, ‘kaiser’, ‘parzen’, ‘projection’, ‘sinogram’, ‘rprojection’, ‘rsinogram’.`For SIRT, we type` x_init = ig.allocate()` ` sirt = SIRT(x_init = x_init, operator = A, data=sino_noisy, constraint = IndicatorBox(lower=0),max_iteration = 50)`` sirt.run(verbose=False) `` sirt_recon = sirt.get_output() ` __(Reminder)__ : If $\mathcal{A} u = g$, then:  $$u^{k+1} = \mathcal{P}_{C}(u^{k} + D \mathcal{A}^{T} M ( g - \mathcal{A} u^{k}))$$ where, $$    \begin{cases}    M = \frac{\mathbb{1}}{\mathcal{A}\mathbb{1}}, \quad m_{ii} = \frac{1}{\sum_{j} a_{ij}}, \quad\mbox{ sum over columns }\\    D = \frac{\mathbb{1}}{\mathcal{A}^{T}\mathbb{1}}, \quad d_{jj} = \frac{1}{\sum_{i} a_{ij}}\quad\mbox{ sum over rows }\\    \end{cases}    $$    and $\mathcal{P}_{C}$ is a projection onto a convex set C. For this demo, we will use a positivity constraint.<jupyter_code># Setup and run the FBP algorithm

fbp = FBP(ig, ag, filter_type = 'hann', device = 'gpu')
fbp.set_input(sino_noisy)
fbp_recon = fbp.get_output()
show(fbp_recon, title = 'FBP reconstruction', cmap = 'inferno')
# Setup and run the SIRT algorithm, with positivity constraint

x_init = ig.allocate()  

sirt = SIRT(x_init = x_init, operator = A, data=sino_noisy, constraint = IndicatorBox(lower=0),
           max_iteration = 50)
sirt.run(verbose=False)

sirt_recon = sirt.get_output()
show(sirt_recon, title = 'SIRT reconstruction', cmap = 'inferno')<jupyter_output><empty_output><jupyter_text>How to setup and run PDHGIn order to setup and run PDHG, we need to define the following:- The operator K.- The function $\mathcal{F}$ and $\mathcal{G}$.- Step-sizes $\sigma$ and $\tau$ such that $\sigma\tau\|K\|^{2}<1$, see [Appendix](Appendix).The setup and run of PDHG:` pdhg = PDHG(f = F, g = G, operator = K, tau = tau, sigma = sigma, max_iterations = maxiter)`` pdhg.run()` We can express our [problem](Lasso) into this [form](PDHG_form), if we let1. $K = \mathcal{A} \quad \Longleftrightarrow \quad $ `K = A`  1. $\mathcal{F}: Y \rightarrow \mathbb{R}, \mbox{ with } \mathcal{F}(z) := \frac{1}{2}\| z - g \|^{2}, \quad \Longleftrightarrow \quad$ ` F = 0.5 * L2NormSquared(sino_noisy)`1. $\mathcal{G}: X \rightarrow \mathbb{R}, \mbox{ with } \mathcal{G}(z) := \alpha\|z\|_{1}, \quad \Longleftrightarrow \quad$ ` G = alpha * L1Norm()`.**Note**: In order to use the PDHG algorithm, we do not need to setup the objective function, but rather the functions $\mathcal{F}, \mathcal{G}$ and operator $K$ that form the objective.Hence, we can verify that with the above setting we end up with this [form](PDHG_form) for $x=u$,  $$\underset{u}{\operatorname{argmin}} \frac{1}{2}\|\mathcal{A} u - g\|^{2}_{2} + \alpha\|u\|_{1} = \underset{u}{\operatorname{argmin}} \mathcal{F}(\mathcal{A}u) + \mathcal{G}(u) = \underset{x}{\operatorname{argmin}} \mathcal{F}(Kx) + \mathcal{G}(x) $$  Define operator K, functions $\mathcal{F}$  and $\mathcal{G}$<jupyter_code>K = A
F = 0.5 * L2NormSquared(b=sino_noisy)
alpha = 0.5
G = alpha * L1Norm()<jupyter_output><empty_output><jupyter_text>Compute operator norm and step-sizes $\sigma, \tau$, such that $\sigma\tau\|K\|^{2}<1$<jupyter_code>normK = K.norm()
sigma = 1/normK
tau = 1/normK<jupyter_output><empty_output><jupyter_text>Setup and run PDHG<jupyter_code># Setup and run PDHG
pdhg = PDHG(f = F, g = G, operator = K, sigma = sigma, tau = tau, 
            max_iteration = 500,
            update_objective_interval = 100)
pdhg.run()

pdhg_l1_recon = pdhg.get_output()<jupyter_output><empty_output><jupyter_text>Show reconstruction result<jupyter_code># Compute PSNR quality measure, and show reconstuction and ground truth
PSNR = psnr(im_data.as_array(), pdhg_l1_recon.as_array())

show(pdhg_l1_recon, title = 'Reconstruction: L1 regularisation: PSNR = {0:.2f}'.format(PSNR), cmap = 'inferno')
show(im_data, title = 'Ground Truth', cmap = 'inferno')

# Plot middle line profile
plt.figure(figsize=(30,15))
plt.rcParams.update({'font.size': 22})
plt.rcParams.update({'lines.linewidth': 5})
plt.plot(np.linspace(1,N,N),im_data.subset(horizontal_y = 128).as_array(), label = 'Ground Truth')
plt.plot(np.linspace(1,N,N),pdhg_l1_recon.subset(horizontal_y = 128).as_array(), label = 'L1 Reconstruction')
plt.legend()
plt.title('Middle Line Profiles')
plt.show()<jupyter_output><empty_output><jupyter_text>Exercise 1: Regularising parameterRun the code above, using different regularisation parameters:1. $\alpha = 1e-6$, $\alpha = 1e6$. What do you observe for these asympotic values of $\alpha$.1. Add more [noise](noise), e.g., __gaussian_var = 5__ and try to find a _good_ regularisation parameter.  Exercise 2: PDHG Convergence We say that the PDHG converges if the __Primal-Dual Gap = Primal objective + Dual objective__ $\rightarrow 0 $,  see  [here](Appendix) . Moreover, the speed of convergence of the PDHG algorithm depends heavily on the choise of step-sizes $\sigma$ and $\tau$. Run the [code](sigma_tau) above with $\sigma=1$, $\tau = \frac{1}{\sigma \|K\|^2}$ and 500 iterations. What do you observe?  Total Variation RegularisationNow, we can continue with the setup of the PDHG algorithm using the Total variation regulariser defined [here](TomoTV). Compared to the [Tikhonov](Tikhonov) problem, we use $L = \nabla$ but with a non-smooth term which is the __MixedL21Norm__ define [here](TomoTV). In addition, we enforce a positivity constraint for our solution using $\mathbb{I}_{\{u>0\}}(u)$. BlockFunctionIn order to setup the above [problem](TomoTV) we need to use the __BlockFramework__ and in particular the concepts of __BlockOperator__, __BlockDataContainer__ and __BlockFunction__. In the previous demos, we have already discussed BlockOperator and BlockDataContainer. BlockFunction behaves similarly to a BlockOperator/BlockDataContainer but instead of _stacking_ operators or DataContainers we _append_ functions as below:$$ F = [f_{1}, f_{2}]$$For example, let $$\begin{align}f_{1}: Y \rightarrow \mathbb{R}, \quad f_{1}(z_{1}) = \alpha\,\|z_{1}\|_{2,1}, \mbox{ ( the TV term ) }\\f_{2}: X \rightarrow \mathbb{R}, \quad f_{2}(z_{2}) = \frac{1}{2}\|z_{2} - g\|_{2}^{2}, \mbox{ ( the data-fitting term ) }\end{align}$$and consider $z = (z_{1}, z_{2})\in Y\times X$, then $$F(z) : = F((z_{1},z_{2}) = f_{1}(z_{1}) + f_{2}(z_{2})$$ Why BlockFunction?With the above form, $F(z)$ is a __separable sum__ of decoupled functions. __So, why do we need to write in this form__? It turns out, that the __proximal operator of a BlockFunction__ can be separated, i.e., $$\mathrm{prox}_{\tau F}(z) = \begin{bmatrix}\mathrm{prox}_{\tau f_1}(z_1)\\\mathrm{prox}_{\tau f_2}(z_2)\end{bmatrix} = \begin{bmatrix}\underset{w_{1}\in Y}{\operatorname{argmin}} \tau f_{1}(w_{1}) + \frac{1}{2}\|w_{1} - z_{1}\|^{2}\\\underset{w_{2}\in X}{\operatorname{argmin}} \tau f_{2}(w_{2}) + \frac{1}{2}\|w_{2} - z_{2}\|^{2}\end{bmatrix}$$Hence, in order to solve the __dual problem__ of the [PDHG](Appendix) algorithm, we need to solve the two decoupled problems as shown above.  PDHG for TV regularisationAt this stage, we have all the necessary ingredients in order to setup our PDHG algorithm for Total variation tomography reconstruction. We define K as a BlockOperator, with the Gradient and Projection operator:$$ K = \begin{bmatrix}\nabla\\A\end{bmatrix}$$` K = BlockOperator(Grad, A)`The function $\mathcal{F}$, is a [BlockFunction](BlockFunction) with 1. a function $\alpha\|\cdot\|_{2,1}\quad\Longleftrightarrow\quad$ ` alpha * MixedL21Norm() ` term that represents the Total variation regularisation ,1. a function $\frac{1}{2}\|\cdot -g \|_{2}^{2}\quad\Longleftrightarrow\quad$ `0.5 * L2NormSquared(sino_noisy)` term that represents the data fitting.Hence, $F = [f1, f2] \quad\Longleftrightarrow\quad $ ` F = BlockFunction( MixedL21Norm(), L2NormSquared(sino_noisy))`Finally, we have the function $\mathcal{G} = \mathbb{I}_{\{u>0\}}(u) \quad\Longleftrightarrow\quad$ ` G = IndicatorBox(lower=0)`Again, we can verify that with the above setting we can express our problem into this [form](PDHG_form), for $x=u$$$\begin{align}\underset{u}{\operatorname{argmin}}\alpha\|\nabla u\|_{2,1} + \frac{1}{2}\|\mathcal{A} u - g\|^{2}_{2} + \mathbb{I}_{\{u>0\}}(u) =  \underset{u}{\operatorname{argmin}} f_{1}(\nabla u) + f_{2}(\mathcal{A}u) + \mathbb{I}_{\{u>0\}}(u) \\ = \underset{u}{\operatorname{argmin}} F(\begin{bmatrix}\nabla \\\mathcal{A}\end{bmatrix}u) + \mathbb{I}_{\{u>0\}}(u) = \underset{u}{\operatorname{argmin}} \mathcal{F}(Ku) + \mathcal{G}(u) = \underset{x}{\operatorname{argmin}} \mathcal{F}(Kx) + \mathcal{G}(x) \end{align}$$  Exercise 2: Setup and run PDHG algorithm with the Total variation regulariserFollow the steps below:1. Define Gradient operator and BlockOperator K. The ImageGeometry, AcquisitionGeometry and Projection Operator and are already defined.1. Define BlockFunction F ( do not forget the regularisation parameter ).1. Define Function G.1. Computer operator norm of K, sigma and tau.1. Setup and run PDHG.1. Show reconstruction result.The reconstruction will look like the figure below after running PDHG for 1000 iterations, with $\sigma=1$, $\tau = \frac{1}{\sigma\|K\|^{2}}$ and $\alpha=2$.<jupyter_code># Define Gradient operator and BlockOperator K

### START CODE HERE (2 lines) ### 
    
    # Grad = Gradient(...)
    # K = BlockOperator(... , ...  )
    
### END CODE HERE ###

# Define BlockFunction F using the MixedL21Norm() and the L2NormSquared() 

### START CODE HERE (4 lines) ### 
    
    # alpha = 2
    # f1 = MixedL21Norm(...)
    # f2 = L2NormSquared(...)
    # F = BlockFunction(...  , ... )
    
### END CODE HERE ###

# Define BlockFunction G, as a positivity constraint

### START CODE HERE (1 lines) ### 
    
    # G = IndicatorBox(...)

### END CODE HERE ###

# Compute operator norm and choose step-size sigma and tau such that sigma*tau||K||^{2}<1 (3 lines)  

### START CODE HERE (3 lines) ### 

    # normK = ...
    # sigma = ...
    # tau = ...

### END CODE HERE ###

# Setup and run PDHG

### START CODE HERE (1 lines) ### 

    # pdhg = PDHG(f = ..., g = ..., operator = ...., tau = ..., sigma = ..., 
    #        max_iteration = 1000, update_objective_interval = 500)
    # pdhg.run(verbose = True)
    # pdhg_tv_recon = pdhg.get_output()

### END CODE HERE ###


# Show reconstruction, compared middle line profiles

# show(pdhg_tv_recon, title = 'TV reconstruction: PSNR = {0:.2f}'.format(PSNR), cmap = 'inferno')
# show(im_data, title = 'Ground Truth', cmap = 'inferno')

# # Plot middle line profiles
# plt.figure(figsize=(30,15))
# plt.rcParams.update({'font.size': 22})
# plt.rcParams.update({'lines.linewidth': 5})
# plt.plot(np.linspace(1,N,N),im_data.subset(horizontal_y = 128).as_array(), label = 'Ground Truth')
# plt.plot(np.linspace(1,N,N),pdhg_tv_recon.subset(horizontal_y = 128).as_array(), label = 'TV Reconstruction')
# # plt.plot(np.linspace(1,N,N),pdhg_l1_recon.subset(horizontal_y = 128).as_array(), label = 'L1 Reconstruction')
# plt.legend()
# plt.title('Middle Line Profiles')
# plt.show()
# Define Gradient Operator and BlockOperator 
Grad = Gradient(...)
K = BlockOperator(...,...)

# Define BlockFunction F using the MixedL21Norm() and the L2NormSquared()
alpha = 2
f1 =  MixedL21Norm(...)
f2 = 0.5 * L2NormSquared(...)
F = BlockFunction(...,...)

# Define BlockFunction G, as a positivity constraint
G = IndicatorBox(...)

# Compute operator norm and choose step-size sigma and tau such that sigma*tau||K||^{2}<1
normK =  ...
sigma = ...
tau = ...

pdhg = PDHG(f = ..., g = ..., operator = ...., tau = ..., sigma = ..., 
            max_iteration = 1000, update_objective_interval = 500)
pdhg.run(verbose = True)

pdhg_tv_recon = pdhg.get_output()

# Compute PSNR quality measure, and show reconstuction and ground truth
PSNR = psnr(im_data.as_array(), pdhg_tv_recon.as_array(), data_range=np.max(im_data.as_array()))

show(pdhg_tv_recon, title = 'TV reconstruction: PSNR = {0:.2f}'.format(PSNR), cmap = 'inferno')
show(im_data, title = 'Ground Truth', cmap = 'inferno')

# Plot middle line profiles
plt.figure(figsize=(30,15))
plt.rcParams.update({'font.size': 22})
plt.rcParams.update({'lines.linewidth': 5})
plt.plot(np.linspace(1,N,N),im_data.subset(horizontal_y = 128).as_array(), label = 'Ground Truth')
plt.plot(np.linspace(1,N,N),pdhg_tv_recon.subset(horizontal_y = 128).as_array(), label = 'TV Reconstruction')
# plt.plot(np.linspace(1,N,N),pdhg_l1_recon.subset(horizontal_y = 128).as_array(), label = 'L1 Reconstruction')
plt.legend()
plt.title('Middle Line Profiles')
plt.show()<jupyter_output><empty_output><jupyter_text>Exercise 3: Find the best regulariser for TV reconstruction and compare with FBP and SIRT algorithmsFor this exercise, we ask the user to1. Run TV reconstruction with 5 different regularising parameters, e.g., $\alpha = [0.05, 0.5, 2, 5, 10]$.1. Choose the best reconstruction according to the highest PSNR.<jupyter_code># Create a list of possible regularisation parameters

### START CODE HERE (1 line) ### 
    
    # alpha = [ , , , , , ]
    
### END CODE HERE ###

# Compute operator norm and choose step-size sigma and tau such that sigma*tau||K||^{2}<1

### START CODE HERE (3 lines) ### 
    
    # normK = ...
    # sigma = ...
    # tau = ...
    
### END CODE HERE ###

# Allocate space for reconstruction results
    # recon_all_geometry = ImageGeometry(channels=len(alpha), voxel_num_x = N, voxel_num_y = N)
    # recon_all = recon_all_geometry.allocate()
    # PSNR = np.zeros(len(alpha))
    
# Define function G, and data-fitting term that is an function-element of the BlockFunction

### START CODE HERE (2 lines) ### 
    
    # G = ...
    # f2 = ...
    
### END CODE HERE ###

# Run PDHG for different regularising parameters

### START CODE HERE (3 lines) ### 

    #for i in range(len(alpha)):

    #    f1 = ...
    #    F = BlockFunction(..., ...)

    #    pdhg = PDHG(f = ..., g = ..., operator = ..., tau = ..., sigma = ..., 
    #            max_iteration = 1000, update_objective_interval = 500)
    #    pdhg.run(verbose = True)

    #    recon_all[i] = pdhg.get_output().as_array()

        # Compute PSNR quality measure, and show reconstuction and ground truth
    #    PSNR[i] =  compare_psnr(im_data.as_array(), pdhg.get_output().as_array(), data_range = 1 )

    #    show(pdhg.get_output(), title = 'TV reconstruction: alpha = {0:.2f}'.format(alpha[i]))
    #    plt.pause(0.5)
    
### END CODE HERE ###

# Print Parameters and PSNR values

    # print('The parameters are {}  '.format(alpha))
    # print('The PSNR values {}  '.format(PSNR))

# Use islicer to scroll over the different reconstructions
    
# from utilities import islicer
# stack_recon = ImageData(np.array([fbp_recon.as_array(), 
#                                   sirt_recon.as_array(),
#                                   pdhg_l1_recon.as_array(),
#                                   recon_all.subset(channel=0).as_array(),
#                                   recon_all.subset(channel=1).as_array(),
#                                   recon_all.subset(channel=2).as_array(),
#                                   recon_all.subset(channel=3).as_array(),
#                                   recon_all.subset(channel=4).as_array(),
#                                   im_data.as_array()]))
# # from ccpi.framework import BlockDataContainer
# # stack_recon = BlockDataContainer(fbp_recon, sirt_recon)
# islicer(stack_recon, 0, title = ['FBP','SIRT','L1 reconstruction',
#                                  'TV a = {}'.format(alpha[0]),
#                                  'TV a = {}'.format(alpha[1]),
#                                  'TV a = {}'.format(alpha[2]),
#                                  'TV a = {}'.format(alpha[3]),
#                                  'TV a = {}'.format(alpha[4]),                                 
#                                  'Ground Truth'], 
#         minmax=(0,1), size=(10,10), slice_number = 0)    
# Setup a list of regularising parameters

# alpha = [ , , , , , ]

# Compute operator norm and choose step-size sigma and tau such that sigma*tau||K||^{2}<1

# K = BlockOperator(..., ...)
# normK = ...
# sigma = ...
# tau = ...

# Allocate space for reconstruction results

# recon_all_geometry = ImageGeometry(channels=len(alpha), voxel_num_x = N, voxel_num_y = N)
# recon_all = recon_all_geometry.allocate()
# PSNR = np.zeros(len(alpha))

# Define BlockFunction G, as a positivity constraint

# G = ...
# f2 = ...

# for i in range(len(alpha)):
      
#     f1 = ...
#     F = BlockFunction(..., ..)

#     pdhg = PDHG(f = ..., g = ..., operator = ..., tau = ..., sigma = ..., 
#             max_iteration = 1000, update_objective_interval = 500)
#     pdhg.run(verbose = True)

#     np.copyto(recon_all.as_array()[i],pdhg.get_output().as_array())

#     # Compute PSNR quality measure, and show reconstuction and ground truth
#     PSNR[i] =  psnr(im_data.as_array(), pdhg.get_output().as_array(), data_range=1)
    
#     show(pdhg.get_output(), title = 'TV reconstruction: alpha = {0:.2f}'.format(alpha[i]), cmap = 'inferno')
#     plt.pause(0.5)
 
print('The parameters are {}  '.format(alpha))
print('The PSNR values {}  '.format(PSNR))
from utilities import islicer
stack_recon = ImageData(np.array([fbp_recon.as_array(), 
                                  sirt_recon.as_array(),
                                  pdhg_l1_recon.as_array(),
                                  recon_all.subset(channel=0).as_array(),
                                  recon_all.subset(channel=1).as_array(),
                                  recon_all.subset(channel=2).as_array(),
                                  recon_all.subset(channel=3).as_array(),
                                  recon_all.subset(channel=4).as_array(),
                                  im_data.as_array()]))
# from ccpi.framework import BlockDataContainer
# stack_recon = BlockDataContainer(fbp_recon, sirt_recon)
islicer(stack_recon, 0, title = ['FBP','SIRT','L1 reconstruction',
                                 'TV a = {}'.format(alpha[0]),
                                 'TV a = {}'.format(alpha[1]),
                                 'TV a = {}'.format(alpha[2]),
                                 'TV a = {}'.format(alpha[3]),
                                 'TV a = {}'.format(alpha[4]),                                 
                                 'Ground Truth'], 
        minmax=(0,1), size=(10,10), slice_number = 0)
plt.figure(figsize=(30,15))
plt.rcParams.update({'font.size': 22})
plt.rcParams.update({'lines.linewidth': 5})
plt.plot(np.linspace(1,N,N),fbp_recon.subset(horizontal_y = 128).as_array(), label = 'FBP', color = 'brown')
plt.plot(np.linspace(1,N,N),fbp_recon.subset(horizontal_y = 128).as_array(), label = 'SIRT', color = 'orange')
plt.plot(np.linspace(1,N,N),pdhg_l1_recon.subset(horizontal_y = 128).as_array(), label = 'L1 Reconstruction', color = 'blue')
plt.plot(np.linspace(1,N,N),recon_all.subset(channel=1, horizontal_y=128).as_array(), label = 'TV Reconstruction', color = 'red')
plt.plot(np.linspace(1,N,N),im_data.subset(horizontal_y = 128).as_array(), label = 'Ground Truth', color = 'black')
plt.legend()
plt.title('Middle Line Profiles')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.153634299750642, 'fraction_numerical': 0.017790018236629575, 'mean_word_length': 3.735636235459993, 'pattern_counts': {'":': 0, '<': 48, '<?xml version=': 0, '>': 50, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 28, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '321864', 'license': 'Apache-2.0', 'n_tokens_mistral': 9284, 'n_tokens_neox': 8641, 'n_words': 2451, 'path': 'training/2019_SynergisticSymposium/04_PDHG_Tomography.ipynb', 'repo_name': 'Sophilyplum/CIL-Demos', 'size': '26879'}