# Week 2

You can access the Zoom meeting [here.](https://ukri.zoom.us/j/99401085708)

Click [here](https://hackmd.io/@SIRF-CIL-Fully3D/r1AxJKNou) to return to the main note.

**Mon 5th of July 14:00 – 16:30 GMT + 01:00**
[Local time and calendar invite](https://dateful.com/eventlink/2183545608)

14:00 – 15:00 Joint session

1. [Introduction to iterative reconstruction and regularisation](https://www.ccpsynerbi.ac.uk/sites/www.ccppetmr.ac.uk/files/20210705_Basics_Iterative.pdf)
2. [Implementation in SIRF/CIL](https://www.ccpsynerbi.ac.uk/sites/www.ccppetmr.ac.uk/files/SIRFCIL%20Course%20Monday%20Week%202%20intro.pdf)

15:00 – 16:30 GMT + 01:00 Sequential sessions for different modalities

1. 15:00-15:30 [CT](https://www.ccpsynerbi.ac.uk/sites/www.ccppetmr.ac.uk/files/CIL_training_week2.pdf)
2. 15:30-16:00 [MR](https://www.ccpsynerbi.ac.uk/sites/www.ccppetmr.ac.uk/files/Week2_Monday_MR_CKolbitsch.pdf)
3. 16:00-16:30 [PET](https://www.ccpsynerbi.ac.uk/sites/www.ccppetmr.ac.uk/files/SIRFCIL%20Course%20recon%20implementation%20-%20Copy.pdf)

**To access the [recording](https://ukri.zoom.us/rec/share/HixAiPpFn78gFdpXk8zptIE84FxRWxsKTcO-28Cy3k-cdql7FuDj8RJijA8OxFBW.5v92NjZG7FcCvMst) of the Zoom meeting use passcode - 3x9x3@3u**

---

**Tue 6th of July 14:00 - 15:00 GMT + 01:00**

Technical support session to get started with the tools for the course. https://ukri.zoom.us/j/99401085708 (Passcode: 1957351)

---

**Wed 7th of July 14:00 - 16:00 GMT + 01:00**
[Local time and calendar invite](https://dateful.com/eventlink/1804088546)
(Partially overlapping) sessions for support

1. 14:00 - 15:00 PET
2. 14:30 - 15:30 MR
3. 15:00 - 16:00 CT

**To access the [recording](https://ukri.zoom.us/rec/share/K73uFR29jmyW5VPmiRDMyg78wICinIOEDgfXMXNA7KlJyUnDF4hr2Q9FzijquIU_.S1bPU3gmLm82dwy0) of the Zoom meeting use passcode - N3g95?=N**

---

**Fri 9th of July 14:00 – 15.30 GMT + 01:00**
[Local time and calendar invite](https://dateful.com/eventlink/1870738523)
Joint session

1. [Summary of main learning objectives](https://www.ccpsynerbi.ac.uk/sites/www.ccppetmr.ac.uk/files/Week2_Friday%20intro.pdf)
2. [PET additions](https://www.ccpsynerbi.ac.uk/sites/www.ccppetmr.ac.uk/files/Week2_Friday_PET_Kris.pdf)
3. Example solutions
4. Overview of next week

**To access the [recording](https://ukri.zoom.us/rec/share/992zxqC0Wx_rxT_8V0g3wuKVRPAxE7eddTwKh8eSXQZWfPdQXosYvDEA-ZCPYfI.ov5A1vWZYlSPCAHU) of the Zoom meeting use passcode - +Nw^L6es**

## Supporting Material
**General**
* [Convex Optimization Short Course](https://web.stanford.edu/~boyd/papers/cvx_short_course.html)
**S. Boyd, S. Diamond, J. Park, A. Agrawal, and J. Zhang** contains slides, jupyter notebooks and three 60-90 minute lectures.
* [General image reconstruction lecture](https://www.youtube.com/watch?v=9d6miyrBvKw)
**Carola Schoenlieb**. This lecture shows how image reconstruction can be described as an inverse problem. Different approaches to solve this problem are presented starting with variational regularisation and leading to deep learning.
* [An introduction to continuous optimization for imaging (Champolle_Pock)](https://hal.archives-ouvertes.fr/hal-01346507/document)
**Antonin Chambolle, Thomas Pock**. Acta Numerica, Cambridge University Press (CUP), 2016, Acta Numerica, 25, pp.161-319. 10.1017/S096249291600009X.
Review paper on state-of-the-art methods.
* [A General Framework for a Class of First Order Primal-Dual Algorithms for Convex Optimization in Imaging Science](https://epubs.siam.org/doi/abs/10.1137/09076934X)
**Ernie Esser, Xiaoqun Zhang, and Tony F. Chan**, SIAM J. Imaging Sci., 3(4), 1015–1046
10.1137/09076934X
PGHG paper

**MR**
- [Basics of MRI and MR image reconstruction](https://www.youtube.com/watch?v=xCv38thzljw)
**Gastao Cruz** (CCP PETMR lecture at PSMR 2018).
Overview of the basics of MRI, from nuclear spins and magnetic fields to the sampling requirements for k-space. Also advanced topics such as reconstruction from undersampled data (GRAPPA, SENSE) are discussed.
- [MR image reconstruction using SIRF](https://www.youtube.com/watch?v=vC66LgPfRNM)
**Christoph Kolbitsch**. MR image reconstruction (k-space sampling, coil sensitivity maps, GRAPPA) and how it can be done in SIRF.
- [Advanced MR image reconstruction ](https://www.youtube.com/watch?v=C5C4juQ5-7M)
**Florian Knoll**. Advanced methods for MR image reconstruction mainly using machine learning.

**CT**
* [CIL paper 1](https://arxiv.org/abs/2102.04560)
* [CIL paper 2](https://arxiv.org/abs/2102.06126)

**PET**
* [Maximum Poisson likelihood PET reconstruction](https://youtu.be/pMx5brQe5yE)
Andrew Reader
* [Maximum likelihood - Expectation Maximisation](https://youtu.be/F7R1PC7Vmt0)
Andrew Reader
* [MAP image reconstruction](https://youtu.be/Uf-I4yaEIx8)
Andrew Reader

## Notebooks
Please check the [Software update instructions](/4pUuyPH1TpKBnbUmrgAgAA) first.

**General**
* [gradient_descent_mr_pet_ct.ipynb](https://github.com/SyneRBI/SIRF-Exercises/blob/master/notebooks/Synergistic/gradient_descent_mr_pet_ct.ipynb)

**MR**
* [d_undersampled_reconstructions.ipynb](https://github.com/SyneRBI/SIRF-Exercises/blob/master/notebooks/MR/d_undersampled_reconstructions.ipynb)
* [e_advanced_recon.ipynb](https://github.com/SyneRBI/SIRF-Exercises/blob/master/notebooks/MR/e_advanced_recon.ipynb)
* [f_create_undersampled_kspace.ipynb](https://github.com/SyneRBI/SIRF-Exercises/blob/master/notebooks/MR/f_create_undersampled_kspace.ipynb)

**CT**
* [01_optimisation_gd_fista.ipynb](https://github.com/TomographicImaging/CIL-Demos/blob/main/training/2021_Fully3D/Week2/01_optimisation_gd_fista.ipynb)
* [02_tikhonov_block_framework.ipynb](https://github.com/TomographicImaging/CIL-Demos/blob/main/training/2021_Fully3D/Week2/02_tikhonov_block_framework.ipynb)
* [03_PDHG.ipynb](https://github.com/TomographicImaging/CIL-Demos/blob/main/training/2021_Fully3D/Week2/03_PDHG.ipynb)
* [04_SPDHG.ipynb](https://github.com/TomographicImaging/CIL-Demos/blob/main/training/2021_Fully3D/Week2/04_SPDHG.ipynb)
* [05_Laminography_with_TV.ipynb](https://github.com/TomographicImaging/CIL-Demos/blob/main/training/2021_Fully3D/Week2/05_Laminography_with_TV.ipynb)
* [additional_exercises_data_resources.ipynb](https://github.com/TomographicImaging/CIL-Demos/blob/main/training/2021_Fully3D/Week2/additional_exercises_data_resources.ipynb)

If **not** running on the STFC cloud you will need to download the data see [CIL README](https://github.com/TomographicImaging/CIL-Demos/blob/c2bae4649b51ecc3f3d0bb77eacc9faeeb5f7503/training/2021_Fully3D/Week2/README.MD)

**PET**
See also the [PET README](https://github.com/SyneRBI/SIRF-Exercises/tree/master/notebooks/PET#week-2)
* [ML_reconstruction.ipynb](https://github.com/SyneRBI/SIRF-Exercises/blob/master/notebooks/PET/ML_reconstruction.ipynb)
* [DIY_OSEM.ipynb](https://github.com/SyneRBI/SIRF-Exercises/blob/master/notebooks/PET/DIY_OSEM.ipynb)
* [MAPEM.ipynb](https://github.com/SyneRBI/SIRF-Exercises/blob/master/notebooks/PET/MAPEM.ipynb)

## How to update SIRF-Exercises

Please refer to [this page](https://hackmd.io/4pUuyPH1TpKBnbUmrgAgAA?view).

## Frequently Asked Questions
Here we will collect questions that are common to many participants. Please only edit this section if you are part of the organisational team.

### I am running out of disk space. What to do?

Some of the exercises generate a fair amount of data, which could be as temporary files. Those are not always cleaned-up properly, especially on the cloud (as the instance can be killed when idle). You will have to delete files that you no longer need. This can be done in via the Jupyter notebook interface, for instance


## User questions
In this section anybody can ask a question. **Please mention the relevant notebook** in your question to make it easier for us to answer. It will also help others to check any previous questions/answers for a specific notebook they are working on. If you think your question does not fit to a specific notebook, please ask it in **General questions.**

## General questions
### Your question
Our answer

### Q: Kullback-Leibler divergence vs Poisson log-likelihood

In Andrew Reader's lecture it was mentioned that the data fidelity term in PET reconstruction (Poisson log-likelihood) is related to the Kullback-Leibler (KL) divergence. The discrete KL divergence is defined as
$$D_{KL}(p,q) = \sum_i p_i \log \frac{p_i}{q_i} $$
Assuming that $p$ is the measured data and $q$ is the estimated data, the Poisson logL is given as
$$-logL(p,q) = \sum_i q_i - p_i \log q_i$$
which can be re-written as
$$-logL(p,q) = \sum_i q_i - p_i \log q_i + \underbrace{(p_i \log p_i - p_i \log p_i)}_{0}$$
$$-logL(p,q) = \sum_i q_i + p_i \log \frac{p_i}{q_i} - p_i \log p_i $$
When optimizing over $q(x)$, the last term can be ignored, so we get
$$-logL^*(p,q) = \sum_i q_i + p_i \log \frac{p_i}{q_i}$$
which is similar to the KL divergence, but not the same.

**Is the negative Poisson logL the same or only similar to the KL divergence?**

(Edo) The Kullback-Leibler implementation in CIL is described [here](https://tomographicimaging.github.io/CIL/nightly/optimisation.html#kullbackleibler)
(Kris) Thanks for the nice formulas! The definition that you gave for the KL divergence as a (generalised) distance between probability distributions, e.g. appropriate for histograms [Wikipedia link](https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence). In optimisation theory, the definition in Edo's link is used, sometimes called the "generalised KL divergence", which is a (generalised) distance between vectors
$${KL}(p,q) = \sum_i p_i \log \frac{p_i}{q_i} - p_i + q_i$$
The latter is indeed the same as negative-Poisson up to terms in $p_i$ (which can be dropped as not dependent on the image $x$). Note that if $p$ and $q$ are (normalised) histograms, $\sum_i p_i=\sum_i q_i = 1$ and the 2 definitions of the KL divergence are the same. It is therefore safest to use the general one!

### Q: Step size ratio in (S)PDHG

(S)PDHG converges given that the product of the primal and dual step size is smaller than equal $1/L^2$ where $L$ is the norm of the forward operator. This in turn mean the ratio of the step sizes is a free parameter (we can multiple one by a factor and divide the other by the same factor).
From my experience, this step size ratio strongly affects the speed of convergence (especailly when using the Poisson log likelihood as data fidelity term).
**Is there any theory on how to choose the step size ratio "wisely" given the data, the forward model and the data fidelity term?**

**A** (Vaggelis): We can accelerate (S)PDHG algorithms using the so-called **strong convexity** property of either the regulariser term or fidelity term or both. This is discussed [here](https://hal.archives-ouvertes.fr/hal-00490826/document#page=14). Under this assumption, we can update the stepsizes sigma, tau analytically per iteration.

In terms of the Kullback Leibler (KL) divergence, we can use a modified version of the KL that is strongly convex and presented [here](https://epubs.siam.org/doi/pdf/10.1137/17M1134834/#page=18).

At the moment, there is no direct implementation in CIL that can handle these two concepts:
a) PDHG under strongly convex functions
b) Modified version of the KL divergence

However, there are two pull requests in CIL that implement the above and can be used in a CIL+SIRF environment:

a) https://github.com/TomographicImaging/CIL/pull/921
b) https://github.com/TomographicImaging/CIL/pull/917

An alternative for now, which is quite easy to implement, is to use a diagonal preconditioning for the sigma, tau stepsizes, presented [here](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.381.6056&rep=rep1&type=pdf). Using **Lemma 2**, we can compute $\sigma$ and $\tau$, not as scalar values, but as diagonal matrices, e.g., $\Sigma, T$

**Note** that, in most real cases we do not have access to a matrix form of the $K$ operator. However, assuming that our operator is **positive** ($Kx\geq0$, for any $x\geq0$) , we can have the following:

$\Sigma = \frac{\mathbb{1}_{Y}}{K * \mathbb{1}_{X}}$

$T = \frac{\mathbb{1}_{X}}{K^{T} * \mathbb{1}_{Y}}$

**Q** Thanks for the help! I don't that any of the accelerations in PDHG can be used, since the gradient of the Poisson logL is not Lipschitz (at least when there are no contaminations.) I also tried the "preconditioned" step sizes as proposed [here](https://epubs.siam.org/doi/pdf/10.1137/17M1134834/#page=18), but the problem is the same. In the paper they have $\gamma$ and $1/\gamma$ in the step sizes, but it is not mentioned how to choose $\gamma$ wisely. And according to my experience $\gamma$ has a huge impact on the speed of convergece. From my heuristic experience so far, a good intial estimated for $\gamma$ is $||x||_\infty$.

**A:** (Vaggelis) Do you use the modified KL divergence?

**Q** Yes, I do. I implemented the approach (the prox) shown by [Ehrhadt et al.](https://iopscience.iop.org/article/10.1088/1361-6560/ab3d07) (follow up paper of the general SPDHG paper).

I added a short demo that compares MLEM, L-BFGS-B and PDHG for optimizing the negative Poisson logL for a 2 voxel system [here](https://github.com/gschramm/teaching/blob/main/nucmed_imaging/logL_optimize.py). This demo nicely shows that effect of $\gamma$ on the speed of convergence of PDHG. Note, that the "toy" demo only requires numpy, scipy and matplotlib.

## Questions about the notebooks
### Your question
Our answer

### Q: (sirf.STIR) Difference between recon.update_current_estimate() and recon.update(arg)

I'm not sure exactly, but I thought the following lines:

recon.update(initial_image)

for i in range(1, num_subiters):
	recon.update_current_estimate()

would do the same as

current_image = initial_image.clone()

for i in range(1, num_subiters + 1)
	current_image = recon.update(current_image)

When I execute the first option though, I run into an error. Could you kindly advise on this? Thank you!

**A** (Kris): You can check the implementation of a function (in iPython) by using `??`. Try

??recon.update

which will tell you that it is defined in terms of `update_current_estimate`. (If you try to see the implementation of the latter, it won't be very instructive as it calls the underlying C/C++ layer. To see that source code, you'd have to go hunting in the C/C++ code).

So I agree with you that these should be equivalent, however the first should start with

recon.set_current_estimate(initial_image)

as otherwise you will be changing `initial_image` (which might be surprising).

Or did you have another error?

### Q: ImageData.add() vs. ImageData.fill()

This may be from last week, but since it still comes up this week I thought I'd still ask.
After generating a noisy array with

` noisy_array = numpy.random.poisson(image.as_array())`

I tried to add it to the image via

	noi_image = image.clone()
	noi_image.add(noisy_array)

While I realize that this is wrong in terms of statistics (the mean would double by adding the noiseless data), I thought it
would still create a noisy image in some sense, but using image.show and noi_image.show and evaluating their maxima tells me that no noise was added. Does this have to do with data types? Thank you in advance!

Note: noi_image.fill(noisy_array) works fine.

**A**

The [`fill`](https://github.com/SyneRBI/SIRF/blob/699b12ceb2e8d74d804c2a3620875d8a969a7843/src/xSTIR/pSTIR/STIR.py#L356-L385) method modifies the object you call the method on: it does copy the values that are in the argument. The argument can be a number, a numpy array or a suitable `DataContainer`. `fill` does not return any value.

The [`add`](https://github.com/SyneRBI/SIRF/blob/699b12ceb2e8d74d804c2a3620875d8a969a7843/src/common/SIRF.py#L179-L206) method adds a number or a `DataContainer` to a `DataContainer`. The behaviour of `add` is different depending on whether the argument `out` is passed, i.e.

- `image.add(other_image, out=yet_another_image)` will store the result of `image + other_image` in a pre allocated `DataContainer` named `yet_another_image`. In this case it will not return anything.
- `image.add(other_image)` will create a new `DataContainer` to store the result and return it.

In your code

noi_image.add(noisy_array)

tries to add a numpy array to a `DataContainer`. There are 2 problems with this:

1. the addition is not defined (did it not [raise](https://github.com/SyneRBI/SIRF/blob/699b12ceb2e8d74d804c2a3620875d8a969a7843/src/common/SIRF.py#L189) an `Error`?)
2. in the case you would have added, the result would have been stored in a new `DataContainer`, but the result is not stored into any variable so after the line executes, the result gets disposed of. Your `noi_image` is not changed

Because the `add` method can be used inline, you could use

noisy_data = noi_image.clone()
noisy_data.fill(noisy_array)
noi_image.add(noisy_data, out=noisy_data)

`help(image.add)` will tell you that it returns the result, so you'd want to do

res = noi_image.add(noisy_array)

However there are some tricky things here. If `noisy_array` is a `numpy.ndarray` and `noi_image` is a `sirf` object, then I believe that `res` will be a `numpy.ndarray` (I haven't tried...)

### ***Question for Prof. Kolbitsch*** #
I am trying to answer your questions...

* If the images from different coils are combined by SOS why do you still need the phase info? Such info was used in reconstructing the image from each single coil. Why do you still need such info when you add the images together?

A: You don't need the phase. A SOS combination removes the phase.

* The disadvantage of SOS, that I can think of, comes from adding the contribution of each coil indiscriminately without taking into account the coil position.
I am afraid I am missing something important here.

A: This is true, but for many practical coils, it is near optimal.

* Simply summing the images shouldn't get rid of random noise?
Once an image is recontructed the value of its voxels is a real number, isn't it?

A: Noise is not removed by adding (if noise is uncorrealted, then adding complex noise can reduce the overall amount of noise). Noise is not completely uncorrelated in parallel imaging.

* I think noise in PI MRI comes from undersampling the K-space and small FOV so aliased replicates for each point.

A: The fundamental source of noise is unchanged. The reonstruction can amplify noise (where the g-factor is high, or similarly when the condition number of the system is high)

* The missing values in the K-space are filled out by interpolation before the Image reconstruction so they will not affect the image quality

* Organ motion during the acquisition will result in aliasing. Tipically the lungs.

### Questions on UnderSampled MRI #
* **Could you please explain the shift of the DC frequency when you perfotm the 2D FFT**?

A: DC is at the first array element in an FFT. This is left over from the days when swapping things around cost time and for some reason having the data arranged this way is more efficient. Its just an annoying convention.

* Q: **Why is the undersampled reconstruction squeezed, but covers the whole FOV?**
A: Because each point in the K-space contains information of the whole scanned object...?
A: The image FOV is 1/(k-space coordinate separation). The pixel size is FOV/N.

* Q: **What artifacts appear in the zero-filled reconstruction?**
A: Artifacts due to truncation of the FFT signal representation...?

* Q: **Why are they artifacts all fine-lined**?
A: Because some phase encoding is missing...?

* Q: **How come they only appear in one direction?**
A: Because some phase encoding is missing while no frequency encoding is missing...?

* Q: **Why are there are artifacts in the reconstruction but not in the coilmaps?**
A: Each coilmap reconstruction comes only from the central part of the K-space that is not undersampled...?

* Q: **In what respect did a GRAPPA reconstruction:
improve the resulting image?
deterioate the resulting image?**
A: It has got rid of the fine-lined artifact. Since it used interpolation to fill out the missing values the artifacts are smeared all over the image...?

### Questions about cost function minimization #
1. The loop exits because the max number of iterations is reached.
I increased the delay in "animation.ArtistAnimation" to observe the evolution of the reconstructed image. I cannot see any relevant change from one step of the iteration to the next one. Why?

Q: **Where is the noise coming from?**
A: The algorithm combined images from each coil that are reconstructed from an under-sampled K-space...?

Q: **Why has not every high-frequency artifact vanished?**
A: The artifacts come from combining together a finite number of aliased images...? I looked up the explanation of SENSE on [http://www.mriquestions.com/senseasset.html](https://)

### What is the e_advanced_recon Preprocessed Data?
* The e_advanced_recon loads 112 K-spce lines from each coil (Yes/No)?
**A:** (Christoph) Yes, that k-space dimensions are 112 phase encoding lines, 4 coils and 512 samples along each readoutline.
* Is the following process sort of a Machine-Learning one? You seem to use "preprocessed_data" as a training set whose result is "grappa_images". The resulting operator "E" can then be used with another acquisition set (Yes/No)?
E= pMR.AcquisitionModel(preprocessed_data, grappa_images)
**A:** (Christoph) No, this has nothing to to with machine learning. In order to set up an MR AcquisitionModel we need the information about the acquired k-space which is provided by `preprocessed_data`and we need the information about which image we want to reconstruct (mainly the dimensions). This is defined by `grappa_images`.

### Question about CIL AcquisitionGeometry
Can I update the rotation axis in the AG without creating a new AG? When I read in data via eg. the txrm reader, the AG is defined with it, and I don't have immediate access to the parameters.

**A** The `rotation_axis` is available in the instance of the `AcquisitionGeometry`, `ag` under `ag.config.system.rotation_axis`. You can set the `position` and `direction` (the second only for 3D) as:

reader = TXRMReader(...)
data = reader.read()
ag = data.geometry

# to change the geometry:
ag.config.system.rotation_axis.position = (X,Y,Z)
ag.config.system.rotation_axis.direction = (v_X,v_Y,v_Z) # a vector

You need to be careful in changing the geometry as some changes might result in an invalid geometry, e.g. if you change the size of the acquisition panel. But in this case it should be OK.

I'm not entirely sure if the update of the geometry triggers a full compatibility check of the geometry and the data. If you want to check that, you could create a new `AcquisitionData` as

ag = old_data.geometry.copy()

ag.config.system.rotation_axis.position = (X,Y,Z)
ag.config.system.rotation_axis.direction = (v_X,v_Y,v_Z) # a vector
new_data = AcquisitionData(old_data.as_array(), geometry=ag)

Notice that the data is passed by reference so there is no duplication of the memory occupation. If this works you can dispose of `old_data` without problems (again Python just removes a reference to the data).

### Question about notebook "**f_create_undersampled_kspace**" ##
I am confused about what kind of data is contained in "preprocessed_data". Can you please clarify?
You reconstruct the image from a fully-sampled K-Space using the content of "preprocessed_data". Then you interrogate such dataset to know which phase-encoding points have been acquired. I would say all the Ky points have been acquired from its usage in the previous step.
**A:** (Christoph) The aim of this notebook is to take a fully sampled data set and then take a subset to create a "new" dataset which now has a certain undersampling pattern.

### How do multiple coils acquire k-space
You have 4 coils.
Does each coil generate a K-space or just a line in a shared K-space?

**A:** (David) All coils receive at the same time. This is why it is often called "parallel imaging". Each coil receives its own k-space. All coils receive at the same k-space coordinates (the coordinates are determined by the imaging gradients which are independent of the coils). If the data is fully sampled (data eventually acquired at all k-space coordinates), then the k-space for each coil is the Fourier Transform of the product of the object with the sensitivity of that coil.

### Which Python development environment to use?
I wonder whether Python has a development workspace like MatLab has. After the school is over I would like to try something on my Mac where Python is installed but I can only call it from a terminal window.
I am not keen on using the notebooks which are valuable tools for teaching.
For instance, C# uses Visual Studio.
Is Python intepreter integrated with some kind of IDE or does it have its own development environment as MatLab has?

**A** There are a range of environments you could use - its personal preference. [This guide](https://www.analyticsvidhya.com/blog/2020/06/5-python-ide-analytics-data-science-programming/) looks like a reasonable summary to me. The before the SIRF training, the STIR training used to use Spyder for exercises, but we have moved away from this. You can try out Jupyter Lab on the cloud instances using [these instructions](https://github.com/SyneRBI/SIRF-Exercises/blob/a38d31188fb2dbd13709977a6dcdc8de41d24a50/DocForParticipants.md#jupyter-lab) (just replace everything from `/tree/` in the URL to `/lab` - `https://training.jupyter.stfc.ac.uk/user/<your user>/lab`). I (Ash) personally use VS Code, however VS Code, PyCharm and PyDev are better suited to library development compared to scientific use.

**A** (Kris) Indeed, there are many IDE environments out there. Spyder still works fine (better than it used to). VS Code has nice features for running the interface locally while the python kernel is actually running remote on your univ system. I'm fairly sure all of them can read Jupyter notebooks as well and act as clients (with possibly a developer- friendlier interface than `jupyter`). I don't do enough Python development myself anymore to have a preference.

**A** (Edo) After some changes, I started using [Visual Studio Code](https://code.visualstudio.com/) which is language agnostic and has a system of extensions that allow you to do many things. Personally I use the Microsoft Python extension, the SSH extension and the CMake extension. With these I am able to run, debug all python code I need. I know that Gemma has used Visual Studio Code to debug a C++ code running remotely.

### Question about PDHG

How is the dual gap calculated?

**A** The primal-dual gap is the difference between the objective of the primal and dual problems. Remember [`PDHG`](https://github.com/TomographicImaging/CIL-Demos/blob/main/training/2021_Fully3D/Week2/03_PDHG.ipynb) solves the following optimisation problem

$$
x = \text{argmin} \ F (K u) + G(u)
$$

The primal dual gap is then defined by $g = p-d$ where

$$
p = F(Kx) + G(x) \\
d = - F^*(y) - G^*(- K^T(y))
$$

where $p$ is the primal objective, $d$ is the dual objective and $F^*$ and $G^*$ the convex conjugates of $F$ and $G$, $x$ is the primal and $y$ the dual (current) solution respectively.

**A** (Georg) Remember that the above only holds if the functions $F$ and $G$ are **proper, convex, lower-semicontinous** functions (which is true for most data fidelity terms and popular priors, but not for all) such that they are equal to their convex biconjugate $(F^{**} = F, G^{**} = G)$ which means that

$$
\text{min}_x \ F (K x) + G(x) \ \ \ \text{(primal problem)}
$$

is equivalent to

$$
\text{min}_x \ \underbrace{\text{max}_y \ \langle Kx, y \rangle - F^* (y)}_{F^{**}(Kx) = F(Kx)} + G(x) \ \ \ \text{(primal-dual problem)}
$$

is also equivalent to (by swapping min and max)

$$
\text{max}_y \ -F^* (y) - G^* (-K^* y) \ \ \ \text{(dual problem)}
$$

#### **follow-up Q** (Kris). Does CIL have a list of priors/data-fidelities that are "proper, convex, lower-semicontinous" and therefore can be used with pdhg? If not, maybe we can make one here...

**A:** Currently, in CIL we have the following functions described in [Table 3](https://arxiv.org/pdf/2102.04560.pdf#page=10).

Most of the functions have an easy `proximal` method, e.g., `L2NormSquared` that can be used in _proximal-type algorithms_. In case we need the `proximal` method of its `convex_conjugate`, i.e., $f^{*}$ we use the **Moreau Decomposition Formula** to compute the `proximal_conjugate`:

$$\mathrm{prox}_{\tau f}(x) + \tau\,\mathrm{prox}_{\tau^{-1}f^{*}}(x/\tau) = x.$$

**Recall:** The _proximal operator of function $f$_ is

$$\mathrm{prox}_{\tau f}(x) := \underset{u}{\operatorname{argmin}}\frac{1}{2}\|x-u\|_{2}^{2} + \tau f(u), \quad\mbox{for any } x.$$

A nice table with a lot of **proximal computations** are reported [here](https://archive.siam.org/books/mo25/mo25_ch6.pdf#page=49).

**Example:** In CIL, we do not have an implementation of the **Affine Function**, see [Section 6.2.2](https://archive.siam.org/books/mo25/mo25_ch6.pdf#page=4), (probably we should) defined as

$$ f(x) = \,<\alpha, x> + b, \,\alpha\in \mathbb{R}^{n}, b\in\mathbb{R} $$

Then,

$$\mathrm{prox}_{\tau f}(x) = x - \tau\alpha$$

In python code, using the CIL design of the `Function Class`, the **AffineFunction** can be implemented as:

from cil.optimisation.functions import Function

class AffineFunction(Function):
    
    def __init__(self, alpha, beta):
        
        self.alpha = alpha
        self.beta = beta
        
    def __call__(self, x):
        
        return self.alpha.dot(x) + self.beta
        
    def convex_conjugate(self, x):
    
        return -self.beta    
    
    def proximal(self, x, tau, out=None):
        
        if out is None:
            
            return x - tau*self.alpha
        
        else:
            
            x.subtract(tau*self.alpha, out=out)

### Is convergence of PDHG always monotonic in terms of the primal and dual, even with correctly set step size can the objectives oscillate?

**A:** No the primal dual gap is not monotonic. To monitor it per iteration, it is better to use `update_objective_interval=1` for the `PDHG` algorithm as:

pdhg= PDHG(f = f, g = g, operator = K, sigma = sigma,
           tau = tau, max_iteration = 1000,
           update_objective_interval = 1)
pdhg.run() 

These lines will run 1000 iterations printing to screen the objective value at each iteration. This may be too much. So you can tweak the amount of prints during `run`

# print every N iterations
N = 10
pdhg.run(print_interval=N)

Then you can retrieve all the objective function values and the iteration number at which that had been evaluated from the algorithm as `list`s with

obj_values = pdhg.objective
iters = pdhg.iterations

Notice that the evaluation of the objective function value is not required by PDHG and it is only performed to monitor the convergence, but it is a costly operation so you may spend a lot of time evaluating the objective value. Therefore we advise to set the `update_objective_interval` to something larger than 1, i.e. evaluate the objective at each iteration. However, there is no way to disable the evaluation of the objective value except using a very high number (higher than `max_iteration` will practically disable it).

Notice that you can still use `print_interval=1` without evaluating the objective function. In this way you can have an idea of the current iteration of the algorithm without getting any information about convergence.

## Question about notebook "b_kspace_filter"
What does the plot entitled "Image Value"represent?
Is it a profile of the image along an axis?
**A:** (Christoph) Yes this plot shows a line of the image. Th term "image value" refers to the pixel value in the image.

## Can TIGRE not (yet?) reconstruct 2D cone beam data?
I tried to extract the centre slice from a 3D cone beam dataset and reconstruct that, but I get `TIGRE main loop fail`.

**A** TIGRE + CIL can handle 2D data. We achieve this by creating a 3D dataset with one slice only [under the hood](https://github.com/TomographicImaging/CIL/blob/39b6f7a722afb6d5f0e2d47a99ce8266378c2a65/Wrappers/Python/cil/plugins/tigre/ProjectionOperator.py#L47-L50).

I modified the [Walnut notebook](https://github.com/TomographicImaging/CIL-Demos/blob/main/training/2021_Fully3D/Week1/01_intro_walnut_conebeam.ipynb) to do work in 2D (on my laptop with Windows10 + CIL from conda)

import numpy as np
import os

from cil.io import TXRMDataReader, TIFFWriter
from cil.processors import TransmissionAbsorptionConverter, \
                           CentreOfRotationCorrector, Slicer
from cil.plugins.tigre import FBP
from cil.utilities.display import show2D, show_geometry

base_dir = os.path.abspath(r"C:\Users\ofn77899\Data\walnut")
data_name = "valnut"
filename = os.path.join(base_dir, data_name, "valnut_2014-03-21_643_28/tomo-A/valnut_tomo-A.txrm")

data3D = TXRMDataReader(file_name=filename).read()

# extract the centre slice
data = data3D.get_slice(vertical='centre') 

show_geometry(data.geometry)

data = TransmissionAbsorptionConverter()(data)

data.reorder(order='tigre')

ig = data.geometry.get_ImageGeometry()

fbp =  FBP(ig, data.geometry)
recon = fbp(data)

show2D(recon, fix_range=(-0.01,0.06))


Interestingly my laptop fails on the reconstruction of the whole dataset (3D) with the following error

tigre/Source/TIGRE_common.cpp (7): Main loop fail 
tigre/Source/TIGRE_common.cpp (14): CBCT:CUDA:Atb unspecified launch failure

In my case binning the 3D dataset made the whole 3D recon work, so I presume is a memory error.

binner = Binner(roi={'horizontal': (None, None, 4),'vertical': (None, None, 4)})
data = binner(data)

### Q follow up...
I did the same, but with a different (Zeiss) dataset. It reconstructs in 3D, but after extracting the centre slice, I get the same two TIGRE errors. I also can't use the CentreOfRotationCorrection (for the same reason I think).

## Question about Coil Sensitivity
* I am afraid I cannot answer the posted question "**Why is there noise in some regions of the coilmaps and outside the object?**"
I get the signal each coil measures is stronger in certain areas and weaker in other ones.
* I would like to go over the other two posted questions:
"**Is this noise in the coilmap going to have a strong negative impact on the image quality in this region of the combined image?**"
Why "Yes" or why "No"?
"**In which organ in the human anatomy would you expect a coilmap to look similarly noisy?**" I daresay any organ where there is motion, like Lungs...?

**A:** (Christoph) If there is no signal in the MR image (e.g. air or lung tissue) then the coil maps also cannot be estimated there and show noise. If there is a good match between coil maps and image (i.e. no motion) then this is usually not a problem.

## Question about DIY_OSEM subset sensitivity
The notebook states, "for a strict implementation of OSEM, you need to compute 'subset sensitivities' for the division." Are we supposed to recalculate the sensitivity image for a given subset? When I explicitely calculate subset_sensitivity in the function below, I end up with a thatched artifact in my output image shown below. The OSEM image was calculated with 4 subsets and one iteration. When I provide the initial calculation of ```sensitivity```, the artifact is not present.

def OSEM(acquired_data, acq_model, initial_image, num_iterations):
    estimated_image=initial_image.clone()

    for subset_num in range(acq_model.num_subsets):
        acq_model.subset_num = subset_num
        
        # Unclear if this is what is meant by compute "subset sensitivities" for the division, it causes artifacts
        subset_sensitivity = acq_model.backward(acquired_data.get_uniform_copy(1))

        for i in range(num_iterations):
            quotient = acquired_data/acq_model.forward(estimated_image)  # y / (Ax + b)
            mult_update = acq_model.backward(quotient)/subset_sensitivity       # A^t * quotient / A^t1
            estimated_image *= mult_update                               # update (in place)
    
    return estimated_image


**A** (Kris) See my answer in the Friday session of this week. In addition, I believe your code is not really OSEM: you are doing `num_iterations` for every subset (i.e. keeping the subset fixed). This generally doesn't work very well. The "trick" with subsets is to let every update use a different subset, such that on average all of the data does get used. In your scheme, any problems of 1 subset will be "baked in", and they then have to be removed by the next subset. I guess that's why you see the "hatch" pattern.
Of course, one can consider various hybrid algorithms (it might be ok to iterate a subset a few times, and you can even think about doing a few subsets in parallel and then combining updates. There's a growing literature on that.

**Followup Q:** I have updated my function to the form:

def OSEM(acquired_data, acq_model, initial_image, num_iterations):
    estimated_image=initial_image.clone()
    
    for i in range(num_iterations):
        for n in range(acq_model.num_subsets):
            acq_model.subset_num = n
            subset_sensitivity = acq_model.backward(acquired_data.get_uniform_copy(1),subset_num=n) # A^t_s * 1
            
            quotient = acquired_data/acq_model.forward(estimated_image,subset_num=n)   # y / (A_s * x + b); there is no way to restrict acquired data to a particular subset
            mult_update = acq_model.backward(quotient,subset_num=n)/subset_sensitivity # A^t_s * quotient / A^t_s*1
            estimated_image *= mult_update                                             # update (in place)
    
    return estimated_image

which still produces the thatched artifact in reconstructed images when I set `num_subsets > 2`. Is there an error in my function? Could someone please provide a working example of their function if they do not see the thatched artifact? Also, if I am updating `acq_model.subset_num = n`, do I still need to explicitely set `subset_num=n` in my function calls?

**A** (Kris) I haven't been able to try your code yet (I hope to get round to it by Wednesday). This looks ok though so I'm not sure. Of course, you're now recomputing subset-sensitivities multiple times, but that's an optimisation to worry about later.
You should not have to add `subset_num=n` in the function call anymore after setting `acq_model.subset_num`. You can of course do a quick experiment to confirm, or check the relevant code with `??`.

**A** (Ash) This appears to be because there are some points at the edge of the FOV that are not captured in all projections.

Check this out:

acq_model.num_subsets = 1
sensitivity = acq_model.backward(acquired_data.get_uniform_copy(1))

acq_model.num_subsets = 4
plt.figure()
for subset in range(acq_model.num_subsets):
    acq_model.subset_num = subset
    subset_sensitivity = acq_model.backward(acquired_data.get_uniform_copy(1))
    # get a boolean image of where total sensitivity is nonzero
    # but the subset sensitivity is zero
    subset_missing = (sensitivity.as_array() != 0) & (subset_sensitivity.as_array() == 0)
    plt.subplot(2, 2, subset+1)
    plt.title(f'subset {subset}')
    plt.imshow(subset_missing[0, :, :])


If you plot the estimated sinogram data at each iteration, you find that there are now undefined (NaN) regions in subsets 2 and 3 (if you call them 0-3):


So we have an issue where some subsets aren't sensitive to these voxels - the multiplicative update then sets them as NaN. However, then the other subsets are sensitive to these voxels so the NaNs are ending up in the sinogram.
Fortunately, these undefined values don't make it into our update (to Kris: why? Because the backprojection treats NaNs in `quotient` as 0?), but they mean that any LORs that run through these regions are ignored, but this wasn't accounted for in the sensitivity calculations. So we have too small of a multiplicative update on these LORs.

Options to fix this would be:
- Set all estimate NaNs to 0 before forward projection
- This is easiest, and it makes sense to assume anything outside the FOV is zero. I assume this might introduce artefacts if you had activity in these voxels.
- Calculate sensitivity using an image of ones inside the FOV of all subsets - then the missing LOR in `backward(quotient)` should be normalised in the division with `subset_sensitivity`

Here is the results of the first solution:

def OSEM(acquired_data, acq_model, initial_image, num_iterations):
    estimated_image=initial_image.clone()

    for i in range(num_iterations):
        for n in range(acq_model.num_subsets):
            acq_model.subset_num = n
            subset_sensitivity = acq_model.backward(acquired_data.get_uniform_copy(1))
            
            # Zero-fill NaNs before forward projection
            estimated_image_arr = estimated_image.as_array()
            estimated_image_arr[numpy.isnan(estimated_image_arr)] = 0
            estimated_image.fill(estimated_image_arr)
            
            quotient = acquired_data/acq_model.forward(estimated_image)
            mult_update = acq_model.backward(quotient)/subset_sensitivity
            
            estimated_image *= mult_update

    return estimated_image


**A**: (Kris) Thanks Ash for the investigation. This is quite useful information. First a note that `sirf.STIR.OSMAPOSLReconstructor` will handle this differently. It does a check in any division and effectively will set NaN to 0 (relevant code is [here](https://github.com/UCL/STIR/blob/master/src/include/stir/numerics/divide.inl)).
Like Ash, I am puzzled that the NaN only creates the hatch pattern. I would have expected that it will generate a whole line of NaNs in the next update, and soon everything becomes unusable. Possibly our `fill` function handles NaN but I didn't think so (it'd be undesirable). Maybe someone can check by doing a `fill` and then `as_array` to see if NaNs have disappeared.
Of course, all of this is a clear illustration why OSEM is not a convergent algorithm, and why it requires balanced subsets to even have a fighting chance... This is why the notebook says "you will need to think about NaNs"...

By the way, there are other anomaly cases where you can see OSEM creating trouble: suppose that a subset is very noisy, and hence the multiplicative update would set some voxels to zero (extreme example: a subset with only 1 LOR with zero counts: all voxels on that LOR would be set to zero). However, future updates will not get rid of those zeroes, even if other subsets disagree. There are 2 (related) solutions to this:
- threshold the multiplicative update (`OSMAPOSLReconstruction.set_minimum_relative_change(.1)` will make sure each update doesn't reduce a voxel to more than 10%). Relevant code is [here](https://github.com/UCL/STIR/blob/877511b607a90cbc3c3faebf9fa4e2dc3ea72f1b/src/iterative/OSMAPOSL/OSMAPOSLReconstruction.cxx#L553-L567)
- use an additive update and threshold it to stay slightly "inside" the non-negative region

### QUESTION ABOUT DOWNLOADING CIL
Please, can someone advise how to download CIL from Mac and integrate it with the other already accessible material without spoiling Docker configuration? Thank you so much.

**A** (Kris) as said on Friday, the Docker image should create the CIL-demos in its `/devel` directory. But if it doesn't, you can do it on the "host" but cloning/copying it yourself into `../SIRF-SuperBuild/docker/devel`

(Edo) **The docker instance has both SIRF and CIL installed**. It should have the [CIL-Demos](https://github.com/TomographicImaging/CIL-Demos) installed in the same directory as `SIRF-Exercises`. If you don't have the `CIL-Demos` you can just do as Kris suggests: clone the repository in `/path_to/SIRF-SuperBuild/docker/devel` where you need to change `path_to` with the path to the SIRF-SuperBuild directory.