# TV reconstruction with high quality CB XCT data

## Background
Data from QMUL (Graham Davis) of tooth. Data acquisition is 4 lengthy high quality CT scans that are reconstructed with (CIL) FDK and averaged to remove noise.

Current acquisition has a slight misalignment and averaging 4 FDK results in blurring.

## Normalisation of the operators

We want to solve the following:
$$\mathrm{argmin}_{u} \frac{1}{2}\|A u - g \|^{2} + \alpha \|\nabla u\|_{2,1}$$

Which we can rewrite as follows:

$$\begin{align} \mathrm{argmin}_{u} \frac{1}{2}\left\| \|A\| \left( \frac{A}{\|A\|} u -\frac{g}{\|A\|} \right) \right\|^{2} + \alpha \left\|\|\nabla\| \frac{\nabla}{\|\nabla\|} u\right\|_{2,1} \\
\|A\|^2 \mathrm{argmin}_{u} \frac{1}{2}\left\| \left( \frac{A}{\|A\|} u -\frac{g}{\|A\|} \right) \right\|^{2} + \alpha \frac{\|\nabla\|}{\|A\|^2}\left\|\ \frac{\nabla}{\|\nabla\|} u\right\|_{2,1}\\
\mathrm{argmin}_{u} \frac{1}{2}\left\| \left( \frac{A}{\|A\|} u -\frac{g}{\|A\|} \right) \right\|^{2} + \alpha \frac{\|\nabla\|}{\|A\|^2}\left\|\ \frac{\nabla}{\|\nabla\|} u\right\|_{2,1}
\end{align}$$

Therefore we have now an L2NormSquared with the argument scaled by $\|A\|$ and a mixed L21 norm scaled by $\|\nabla\|$. In SPDHG we deal with subsets, and in the case each subset is the same we can rewrite the equation above as

$$\mathrm{argmin}_{u} \frac{1}{2} \sum_i^{n-1}\left\| \left( \frac{A_i}{\|A_i\|} u -\frac{g_i}{\|A_i\|} \right) \right\|^{2} + \alpha \frac{\|\nabla\|}{\|A_i\|^2}\left\|\ \frac{\nabla}{\|\nabla\|} u\right\|_{2,1}+ \frac{1}{2} \frac{\|A_n\|}{\|A_i\|^2} \left\| \left( \frac{A_n}{\|A_n\|} u -\frac{g_n}{\|A_n\|} \right) \right\|^{2}$$

The last term will be equal to each member of the sum, if the $n$-th physical subset has the same geometry of the others, and the sum could be extended from 1 to $n$. In case the $n$-th subset has a different geometry from the other subsets the last term takes care of the difference.

## Subsets

### Number of subsets

Purely based on computational resources requirements, I would select the number of subsets so that the GPU is utilised for longer. This means that one has to monitor the GPU memory/computation usage and find the good spot where computation is longer and most of the memory is used. **I think**.

The TV subset will not depend on the physical subsets and it will take a long time.

**Are there mathematical considerations that would suggest how many physical subsets to choose?**

Matthias in https://doi.org/10.1117/12.2272946 states that, by visual investigation, one can observe that ['the artifacts get reduced with more subsets indicating that with more subsets the algorithm converges faster'](https://www.spiedigitallibrary.org/conference-proceedings-of-spie/10394/2272946/Faster-PET-reconstruction-with-a-stochastic-primal-dual-hybrid-gradient/10.1117/12.2272946.full#ID0EUFAE). Of course this is _faster_ to converge wrt to number of epochs, not clock time. What's not clear is the time it takes an epoch in the case of different amount of subsets.

### Selection
SPDHG in CIL is configured to have equal selection probability for each subset, aka uniform sampling. This means that, doing explicit TV with 25 subsets, results in a probability of 1/26 of selection of each subset. In this case, the TV term is going to be evaluated very few times in many iterations.

In my experiments I used non uniform sampling, and I increased the probability of the TV term to 25%, as seen in the picture below where the histogram of the frequency of selection of each subset is plotted. Each physical subset has a 3% probability to be selected, 75%/25 = 3%


In [Ehrhardt et al](https://www.spiedigitallibrary.org/conference-proceedings-of-spie/10394/2272946/Faster-PET-reconstruction-with-a-stochastic-primal-dual-hybrid-gradient/10.1117/12.2272946.full#s4) cited above, with a similar set-up with TV regularisation, the sampling is non-uniform with 50% probability of selection of a physical subset.

# Reconstructions

## Dataset 1

First Dataset from Graham. It had problems with centering of the 4 scans resulting in blurring of the average scan. For this reason here the comparison is between the data of the second scan only.

CENTRE:                 652.942383
CENTRE:                 652.960449
CENTRE:                 652.932617
CENTRE:                 652.919556

## Alpha = 0.0003

Clearly over regularised. After 100 epochs the tiny details are lost.

'args': {'--alpha': '0.0003', '--normalise-operators': True, 
         '--num_phys_subs': '25', '--adaptive-georg': False, 
         '--output_dir': '/mnt/data/edo/Dev/recon/SPDHG/QM_tooth/normalised/'}



## Alpha 0.0001

# TODO

- do recon with 50 Z slices
- compare reconstruction with 50 slices with the one with full data and try to see what happens with $\alpha$, if we can estimate it in small scale and then run in full scale.
- maybe use simulated dataset in dataexample
- what is voxel size of tooth?
- how to link alpha parameter for 2D centre slice to 3D

Graham's reconstruction is in a ‘.siz’ file. This contains the dimensions of the bin file in x, y and z sizes. The data is 32 bit floating point, little endian, incrementing in x, y then z (C contiguous).

https://petpp.github.io/data/slides/2021_MIC_Delplancke.pdf

$$
\frac{I}{I_0} = \exp {\int_{L_i} -\mu (s) ds}
$$

$$
b_i = -\log \frac{I_i}{I_0} = \int_{L_i} \mu(s) ds
$$

$$ \text{CCP}i $$ which is different from CCP*i*

$I$ $I_0$

$$
b_i = \sum_j a_{ij} u_j
$$

$$
f =
\begin{pmatrix}
||\cdot - \texttt{b} ||_2^2\\
|| \cdot ||_{2,1}
\end{pmatrix}
$$

$$
\texttt{K} =
\begin{pmatrix}
A \\
\nabla
\end{pmatrix}
$$