# Detection and estimation of rigid patient motion during SPECT acquisition with exponential data consistency conditions

My Hoang Hoa Bui, Antoine Robert, Ane Etxebeste, Simon Rit

► To cite this version:

My Hoang Hoa Bui, Antoine Robert, Ane Etxebeste, Simon Rit. Detection and estimation of rigid patient motion during SPECT acquisition with exponential data consistency conditions. Fully 3D Image Reconstruction in Radiology and Nuclear Medicine, May 2025, Shanghai, China. hal-05375603

HAL Id: hal-05375603

https://hal.science/hal-05375603v1

Submitted on 21 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Detection and estimation of rigid patient motion during SPECT acquisition with exponential data consistency conditions

My Hoang Hoa Bui¹, Antoine Robert¹, Ane Etxebeste¹, and Simon Rit¹

¹INSA-Lyon, Université Claude Bernard Lyon 1, CNRS, Inserm, CREATIS UMR 5220, U1294, F-69373, Lyon, France

**Abstract** Rigid patient motion during single photon emission computed tomography (SPECT) acquisition can introduce artifacts in the reconstructed SPECT images. While translational motion can be estimated by data-driven methods with subpixel accuracy, rotations are more challenging. With a known convex region K of constant attenuation coefficient encompassing all injected activity, we can derive the exponential projections from SPECT measurements and evaluate their pair-wise exponential data consistency conditions (eDCCs). This study aimed to propose a fully data-driven method using eDCCs to detect and estimate rigid patient motion. The proposed method was evaluated on simulated data of a liver patient and allowed a great improvement in the reconstructed image quality with an activity recovery coefficient above 95% with respect to the reconstructed images from motion-free projections.

## 1 Introduction

Rigid patient motion during single photon emission computed tomography (SPECT) acquisition is likely to occur due to the long acquisition time, up to 40 minutes, which can introduce artifacts in the reconstructed SPECT images. Exponential data consistency conditions (eDCCs) are a set of mathematical equations characterizing the range of the exponential Radon transform. In our previous work, we showed the ability of a data-driven method based on eDCCs to detect and estimate translational motion in SPECT with sub-pixel accuracy [1].

Rigid motion is likely to also have a rotational component which can be corrected mainly with reconstruction-based methods [2, 3] or with external hardware devices [4]. However, the former may be computationally expensive or have limited accuracy while the latter can be costly and time-consuming.

eDCCs may offer an appealing alternative solution to detect and estimate rigid patient motion as they operate in the projection domain and they do not require any additional hardware. This work aimed to propose and evaluate a fully data-driven solution for rigid motion estimation in SPECT using eDCCs.

## 2 Materials and methods

### 2.1 Exponential data consistency conditions

In three-dimensional (3D) parallel SPECT, the measurement of an activity distribution f and an attenuating medium μ, can be modeled by the 2D attenuated Radon transform in each detector line v:

$$ p_{\lambda}^{\mu}(u, v) = \int_{-\infty}^{+\infty} f(u\alpha_{\lambda} + v\beta_{\lambda} + t\gamma_{\lambda}) \exp\left(-\int_{t}^{+\infty} \mu(u\alpha_{\lambda} + v\beta_{\lambda} + t'\gamma_{\lambda}) dt'\right) dt \quad (1) $$

where ($\alpha_{\lambda}, \beta_{\lambda}, \gamma_{\lambda}$) defines the detector orientation as in Fig. 1 with $\alpha_{\lambda} = (\cos \lambda, \sin \lambda, 0)$, $\beta_{\lambda} = (0, 0, 1)$ and $\gamma_{\lambda} = (-\sin \lambda, \cos \lambda, 0)$.

If there is a known convex region K where the attenuation is constant, i.e. $\mu(x) = \mu_0$ if $x \in K$, and all radioactivity is encompassed inside K, i.e. $f(x) = 0$ if $x \notin K$, exponential projections can be derived

$$ p_{\lambda}(u, v) = p_{\lambda}^{\mu}(u, v)C_{\lambda}(u, v). \quad (2) $$

The conversion factor $C_{\lambda}(u, v)$ was introduced in [5]

$$ C_{\lambda}(u, v) = \exp \left( \tau_{\lambda,u,v} \mu_0 + \int_{\tau_{\lambda,u,v}}^{+\infty} \mu(u\alpha_{\lambda} + v\beta_{\lambda} + t'\gamma_{\lambda}) dt' \right) \quad (3) $$

with $\tau_{\lambda,u,v}$ the coordinate where the photon leaves the K region along the integral line towards the $(u, v)$ detector pixel. eDCCs state that for $\gamma_{\lambda} \cdot \gamma_{\lambda'} \neq -1$, $p_{\lambda}(\cdot, v)$ and $p_{\lambda'}(\cdot, v)$ are consistent if and only if

$$ P_{\lambda}(\sigma_{\lambda,\lambda'}, v) = P_{\lambda'}(\sigma_{\lambda',\lambda}, v) \quad (4) $$

with $P_{\lambda}(\sigma_{\lambda,\lambda'}, v) = \int_{-\infty}^{+\infty} p_{\lambda}(u, v) \exp(\sigma_{\lambda,\lambda'}u) du$ the two-sided Laplace transform of projection line $p_{\lambda}(\cdot, v)$ and

$$ \sigma_{\lambda,\lambda'} = -\sigma_{\lambda',\lambda} \quad (5) $$

$$ = \mu_0 \cdot \operatorname{sgn}((\alpha_{\lambda} + \alpha_{\lambda'}) \cdot (\gamma_{\lambda'} - \gamma_{\lambda})) \frac{\|\gamma_{\lambda} \times \gamma_{\lambda'}\|}{\gamma_{\lambda} \cdot \gamma_{\lambda'} + 1} \quad (6) $$

When there is no rigid motion between $\lambda$ and $\lambda'$, $\sigma_{\lambda,\lambda'} = -\sigma_{\lambda',\lambda} = \mu_0 \tan((\lambda - \lambda')/2)$ as previously defined in [1, 6].

<figure><img src="image_2.png" /><figcaption>Figure 1: Schematic illustration of the pair-wise geometry for acquired projections $p_{\lambda}^{\mu}$ and $p_{\lambda'}^{\mu}$, when there is no rigid motion.</figcaption></figure>

### 2.2 eDCCs with a known rigid motion

Suppose that the object underwent a rigid motion between the acquisition of projections $\lambda$ and $\lambda'$, characterized by a 3D rotation matrix $A$ (parameterized by its Euler angles defining
sequential rotations about x, y and z) and a 3D translation **a**. The emission and attenuation maps after the rigid motion are $\tilde{f}(\mathbf{x}) = f(\mathbf{A}\mathbf{x} + \mathbf{a})$ and $\tilde{\mu}(\mathbf{x}) = \mu(\mathbf{A}\mathbf{x} + \mathbf{a})$, respectively.

Since $p_{\lambda}^{\mu}$ and $p_{\lambda'}^{\mu}$—acquired at $(\alpha_{\lambda'}, \beta_{\lambda'}, \gamma_{\lambda'})$ with the rigid motion—are now measurements of two different objects, they are likely to be inconsistent [6]. Such motion is equivalent to the acquisition of a projection $\tilde{p}_{\lambda'}^{\mu}$ of the object before motion with a different detector geometry which can be defined if the rigid motion is known (Fig. 2). Specifically, the corresponding geometry of the detector is an orientation $(\tilde{\alpha}_{\lambda'}, \tilde{\beta}_{\lambda'}, \tilde{\gamma}_{\lambda'}) = (\mathbf{A}\alpha_{\lambda'}, \mathbf{A}\beta_{\lambda'}, \mathbf{A}\gamma_{\lambda'})$, followed by a translation $(\Delta u, \Delta v) = (-\mathbf{a} \cdot \mathbf{A}\alpha_{\lambda'}, -\mathbf{a} \cdot \mathbf{A}\beta_{\lambda'})$. The transformed detector orientation accounts for the object's rotation, while the detector shift accounts for the object's translation [1].

Although $p_{\lambda}^{\mu}$ and $\tilde{p}_{\lambda'}^{\mu}$ are measurements of the same object $f$, their corresponding rows $p_{\lambda}^{\mu}(\cdot, v)$ and $\tilde{p}_{\lambda'}^{\mu}(\cdot, v)$ are no longer imaging the same part of $f$. However, there is a set of planes, parallel to $\gamma_{\lambda}$ and $\tilde{\gamma}_{\lambda'}$ (i.e. orthogonal to $\gamma_{\lambda} \times \tilde{\gamma}_{\lambda'}$) which image the same 2D slice of $f$, e.g. the plane $\mathcal{P}$ in green in Fig. 2.

<figure><img src="image_3.png" /><figcaption>**Figure 2:** Schematic illustration of the pair-wise eDCC with the new detector geometry for $\lambda'$ (in dark gray) accounting for the rigid motion. $\mathcal{P}$ is the plane defined by the origin and the two unit vectors $\gamma_{\lambda}$ and $\tilde{\gamma}_{\lambda'}$.</figcaption></figure>

If the motion is known, pairwise eDCCs can therefore be calculated along the two lines at the intersection of these planes and the two detectors. These pairs of lines can be described in two new detector coordinate systems $(\tilde{\alpha}_{\lambda,\lambda'}, \tilde{\beta}_{\lambda,\lambda'}, \gamma_{\lambda})$ and $(\tilde{\alpha}_{\lambda',\lambda'}, \tilde{\beta}_{\lambda',\lambda'}, \tilde{\gamma}_{\lambda'})$ with

$$ \tilde{\beta}_{\lambda,\lambda'} = \tilde{\beta}_{\lambda',\lambda} = \frac{\gamma_{\lambda} \times \tilde{\gamma}_{\lambda'}}{\|\gamma_{\lambda} \times \tilde{\gamma}_{\lambda'}\|} \quad (7) $$

and

$$ \tilde{\alpha}_{\lambda,\lambda'} = \gamma_{\lambda} \times \tilde{\beta}_{\lambda,\lambda'}; \quad \tilde{\alpha}_{\lambda',\lambda} = \tilde{\gamma}_{\lambda'} \times \tilde{\beta}_{\lambda,\lambda'} \quad (8) $$

The rotated and aligned projections according to these new coordinate systems are

$$ \tilde{p}_{\lambda,\lambda'}^{\mu}(u,v) = p_{\lambda}^{\mu} \left( \tilde{\mathbf{M}}_{\lambda,\lambda'} \begin{bmatrix} u \\ v \end{bmatrix} \right) \quad (9) $$

$$ \tilde{p}_{\lambda',\lambda}^{\mu}(u,v) = p_{\lambda'}^{\mu} \left( \tilde{\mathbf{M}}_{\lambda',\lambda} \begin{bmatrix} u - \mathbf{a} \cdot \mathbf{A}\alpha_{\lambda'} \\ v - \mathbf{a} \cdot \mathbf{A}\beta_{\lambda'} \end{bmatrix} \right) \quad (10) $$

with the 2D rotation matrices

$$ \tilde{\mathbf{M}}_{\lambda,\lambda'} = \begin{bmatrix} -\tilde{\alpha}_{\lambda,\lambda'} & - \\ -\tilde{\beta}_{\lambda,\lambda'} & - \end{bmatrix} \begin{bmatrix} \alpha_{\lambda} & \beta_{\lambda} \\ - & - \end{bmatrix} \quad \text{and} \quad (11) $$

$$ \tilde{\mathbf{M}}_{\lambda',\lambda} = \begin{bmatrix} -\tilde{\alpha}_{\lambda',\lambda} & - \\ -\tilde{\beta}_{\lambda',\lambda} & - \end{bmatrix} \begin{bmatrix} \mathbf{A}\alpha'_{\lambda} & \mathbf{A}\beta'_{\lambda} \\ - & - \end{bmatrix}. \quad (12) $$

The eDCCs (4) can be evaluated for each row of the rotated projections $\tilde{p}_{\lambda,\lambda'}^{\mu}$ and $\tilde{p}_{\lambda',\lambda}^{\mu}$ after the conversion to exponential projections using (3) with the corresponding detector geometry.

## 2.3 Rigid motion estimation

For the sake of simplicity, we only considered a single inter-frame rigid patient motion right before projection $\lambda_0$, i.e. $\lambda_0$ divides the acquisition into no-motion projections $p_{\lambda}^{\mu}$ for $\lambda < \lambda_0$ and projections with motion $p_{\lambda'}^{\mu}$ for $\lambda' \ge \lambda_0$. To estimate $\lambda_0$ and the rigid patient motion parameters **A** and **a**, we minimized the inconsistency between projections by solving the minimization problem

$$ (\hat{\lambda}, \hat{\mathbf{A}}, \hat{\mathbf{a}}) = \underset{\lambda \in \Omega_{\lambda}}{\operatorname{argmin}} F(\lambda, \mathbf{A}, \mathbf{a}) \quad (13) $$

with $\Omega_{\lambda}$ the set of gantry positions of the acquisition, the cost function

$$ F(\lambda, \mathbf{A}, \mathbf{a}) = \frac{1}{K} \sum_{\lambda' < \lambda \le \lambda'} \mathrm{NE}(\lambda, \lambda'), \quad (14) $$

$K$ the number of evaluated projection pairs $\tilde{p}_{\lambda,\lambda'}^{\mu}$ and $\tilde{p}_{\lambda',\lambda}^{\mu}$, (fulfilling $\lambda < \lambda' \le \lambda$ and $\gamma_{\lambda} \cdot \tilde{\gamma}_{\lambda'} \ne -1$) and the noise-aware metric, $\mathrm{NE}(\lambda, \lambda')$ described in our previous works [1, 6]. The problem was solved by first calculating

$$ (\hat{\mathbf{A}}_{\lambda}, \hat{\mathbf{a}}_{\lambda}) = \underset{\mathbf{A} \in \mathrm{SO}(3), \mathbf{a} \in \mathbb{R}^3}{\operatorname{argmin}} F(\lambda, \mathbf{A}, \mathbf{a}) \quad (15) $$

for every $\lambda \in \Omega_{\lambda}$ with the Nelder-Mead downhill simplex algorithm [7] and then taking the $\lambda$ value which minimizes $F$, i.e.

$$ (\hat{\lambda}, \hat{\mathbf{A}}_{\lambda}, \hat{\mathbf{a}}_{\lambda}) = \underset{\lambda \in \Omega_{\lambda}}{\operatorname{argmin}} F(\lambda, \hat{\mathbf{A}}_{\lambda}, \hat{\mathbf{a}}_{\lambda}). \quad (16) $$

In favor of minimizing the computation time, we defined a constraint search range of the optimizer, i.e. $[-5, 5]$ degrees for the **A** rotation parameters and $[-20, 20]$ mm for the **a** translation components. Moreover, all six motion parameters were normalized to the range $[0, 1]$ to homogenize the optimization step size across all optimized variables.

## 2.4 Data acquisition and motion simulation

We simulated SPECT projections with a voxelized CT image of a liver patient and an ellipsoid source of size (20, 40, 30) mm located inside the liver. The $^{99m}$Tc activity was 120 MBq. With RTK [8] and GATE [9], different sets of projections were generated to model different physical effects in the SPECT system: RTK attenuated projections modeling (1); RTK attenuated with Poisson noise; RTK attenuated, Poisson noise and PSF (detector spatial resolution);
GATE simulated projections modeling all physical effects. One SPECT acquisition consists of 60 step-and-shoot projections covering 360° circular orbit around the z axis of radius 38 cm. Each projection image was generated for 20 s with 256 × 256 pixels of 2.39 mm isotropic spacing. Data were acquired within a primary energy window (15% centered around 140.5 keV) and a lower scatter window (15% centered around 120 keV) and scatter corrected with the dual-energy window method [10].

One set of projections, considered as a reference, was acquired with RTK or GATE simulations without motion. Another set was acquired after simultaneously rotating the attenuation and emission maps 2° around each axis, then translating by (10, 10, 10) mm. From these two sets, we can generate various acquisitions with inter-frame rigid motion occurring at different projections $\lambda_0$.

The convex region K was the whole liver segmented in the original attenuation map.

## 2.5 Quantitative assessment of motion correction

The estimated motion was defined at the minimum of 20 optimized cost functions with different initializations. This process was then repeated 10 times to evaluate the robustness of the proposed method. The six parameters of the estimated motion and the shifted source center, defined as $\hat{\mathbf{s}}_0 + \hat{\mathbf{a}}$ with $\mathbf{s}_0$ the center of the simulated ellipsoid, were compared to their respective ground truth.

The motion-free, non-corrected and motion-corrected projections (accounting for the rigid motion in the acquisition geometry) were reconstructed by the ordered subset expectation maximization (OSEM) implementation of RTK with attenuation correction, using 5 iterations and 4 subsets. The image quality was evaluated by the relative activity recovery coefficient

$$rARC(\%) = \frac{A}{A_{\text{motion-free}}} \times 100\%. \quad (17)$$

A and $A_{\text{motion-free}}$ are the activities in the source region obtained from evaluated projections (non-corrected or motion-corrected) and from motion-free projections, respectively.

## 3 Results

Fig. 3 illustrates the estimation error of the rigid motion occurring right before different projections $\lambda_0$ during the acquisition for RTK attenuated projections (left column) and GATE simulated projections (right column). Without simulating physical effects in the SPECT system, the performance of motion correction with eDCCs was quite robust across all studied values of $\lambda_0$. The proposed method returned sub-pixel accuracy for translational components (top panel) and sub-degree accuracy for rotational components (middle panel). The variation of motion error across different repetitions remained stable and relatively small for all components (narrow and overlayed shaded areas) with RTK simulations. Compared to RTK attenuated projections, the bias in translation and rotation estimation of GATE simulated data was more pronounced, except for z-rotation. Although the proposed method with eDCCs was less accurate in estimating

<figure><img src="image_4.png" /><figcaption>**Figure 3:** Average motion estimation error over 10 repetitions of the proposed method as a function of different motion projection $\lambda_0$ for respectively RTK attenuated and GATE simulated projections. The colored background represents one standard deviation of 10 repetitions.</figcaption></figure>

motion for GATE simulated data, the offset in the center of the emitting source as compared to the reference position was still recovered with 2-pixel accuracy (black curves) while this value was below 1 mm for RTK attenuated projections. To further evaluate the influence of various physical effects on the performance of eDCCs motion correction, the averages of the optimized cost function values and the offsets of the source center across 10 repetitions are presented in Table 1. From top to bottom, the optimized value of the cost function increases with the addition of physical effects which are not modeled by (1). For RTK-simulated projections with Poisson noise or even with detector spatial resolution, sub-pixel accuracy was achieved with minimal fluctuations. In contrast, for GATE-simulated projections, despite the application of scatter correction, the bias was more pronounced. Nevertheless, the variation remained within sub-pixel range.

<table><caption>**Table 1:** Average optimized cost function and misplacement of the emission map center over 10 repetitions of motion estimation, for mid-acquisition rigid motion on different projection sets.</caption><thead><tr><th>Projections</th><th>F(\hat{m}, \tilde{\mathcal{R}}_m)<br>[unitless]<br>mean ± std</th><th>Source center offset<br>[mm]<br>mean ± std</th></tr></thead><tbody><tr><td>RTK attenuated</td><td>0.050 ± 0.005</td><td>0.3 ± 0.3</td></tr><tr><td>RTK attenuated, noise</td><td>0.157 ± 0.002</td><td>0.2 ± 0.1</td></tr><tr><td>RTK attenuated, noise, PSF</td><td>0.289 ± 0.001</td><td>0.8 ± 0.2</td></tr><tr><td>GATE scatter-corrected</td><td>0.405 ± 0.002</td><td>4.1 ± 0.7</td></tr><tr><td>GATE simulated</td><td>0.437 ± 0.001</td><td>3.7 ± 0.5</td></tr></tbody></table>

Qualitative assessment of the reconstructed SPECT image is shown in Fig. 4 for GATE simulated projections. Motion correction shows an enhancement in image quality as compared to the reference case without motion with a strong blur
reduction and restoration of the orientation of the ellipsoid.

<figure><img src="image_5.png" /><figcaption>**Figure 4:** SPECT images obtained from motion-free (left), uncorrected (middle) and motion-corrected (right) projections simulated with GATE for a mid-acquisition motion.</figcaption></figure>

To quantitatively assess the reconstructed SPECT image after correction for rigid motion with eDCCs, the rARC before and after motion correction of GATE simulated data is provided side-by-side in Fig. 5. Without correction, the values of rARC computed on reconstructed images from GATE simulated projections spread in a wide range from 82% to nearly 100%. In contrast, motion correction brought a comparable ARC to motion-free images, above 95% for all studied cases.

<figure><img src="image_6.png" /><figcaption>**Figure 5:** Relative activity recovery coefficient (rARC) of GATE-simulated projections as a function of motion projection $\lambda_0$.</figcaption></figure>

## 4 Discussion

In this study, eDCCs were employed for the first time to detect and estimate 3D rigid patient motion, including both the 3D translation and rotation. The method demonstrated the robustness of detecting sudden inter-frame rigid motions in idealized SPECT measurements simulated with RTK. For motion occurring at different projections $\lambda_0$, motion correction applied to RTK-attenuated projections showed a good agreement with the ground truth. Translational components achieved sub-pixel accuracy, while rotational errors were below one degree. In contrast, motion estimated from GATE-simulated projections was not fully accurate. In both cases, the method obtained a small variance among different repetitions, shown as thin-colored regions in Fig. 3. Interestingly, analysis of the emission map displacement revealed that many of these estimates were spatially close to the true emission map position, typically within two pixels of the source center.

Consistent with findings reported in [6], the performance of eDCC-based motion correction was found to be slightly degraded in the presence of Poisson noise and point spread function (PSF) effects, but significantly degraded by scatter, as shown in Table 1.

The current framework with Nelder-Mead downhill simplex optimization is computationally expensive. We expect a significant computation time reduction with gradient-based optimization.

Despite these limitations, eDCC-based rigid motion correction was highly effective in restoring reconstructed image quality. Specifically, the proposed method demonstrated the ability to restore image quality to nearly 100% activity of the reference SPECT image without motion.

## 5 Conclusions

In this study, eDCCs were used to detect and correct for rigid patient motion, including both the 3D translation and rotation. While the optimizer struggled to accurately estimate the true motion parameters due to the presence of scatter, it demonstrated its ability to reconstruct SPECT images with a quality comparable to that obtained from motion-free projections.

## Acknowledgment

This work was supported by grant ANR-21-CE45-0026 (SPECT-Motion-eDCC) from the French National Research Agency (ANR), performed within the framework of the SIRIC LYRICAN INCANSERM-DGOS-12563 and the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d'Avenir" (ANR- 11-IDEX-0007) operated by the ANR.

## References

[1] M. H. H. Bui et al. "Detection and correction of translational motion in SPECT with exponential data consistency conditions". *Physics in Medicine & Biology* (2025). DOI: 10.1088/1361-6560/adb09a.

[2] B. Hutton et al. "A hybrid 3D reconstruction/registration algorithm for correction of head motion in emission tomography". Vol. 2. 2000, 15/36–15/40 vol.2. DOI: 10.1109/NSSMIC.2000.950043.

[3] R. R. Fulton et al. "Use of 3D reconstruction to correct for patient motion in SPECT". *Physics in Medicine & Biology* 39.3 (1994), p. 563. DOI: 10.1088/0031-9155/39/3/018.

[4] J. E. McNamara et al. "A flexible multicamera visual-tracking system for detecting and correcting motion-induced artifacts in cardiac SPECT slices". *Medical Physics* 36.5 (2009), pp. 1913–1923. DOI: https://doi.org/10.1118/1.3117592.

[5] F. Natterer. *The Mathematics of Computerized Tomography*. NY: John Wiley & Sons, 1986.

[6] A. Robert et al. "Ability of exponential data consistency conditions to detect motion in SPECT despite other physical effects". *17th International Meeting on Fully 3D Image Reconstruction in Radiology and Nuclear Medicine*. Stony Brook, NY, United States, 2023, pp. 337–340.

[7] J. A. Nelder et al. "A Simplex Method for Function Minimization". *The Computer Journal* 7.4 (1965), pp. 308–313. DOI: 10.1093/comjnl/7.4.308.

[8] S. Rit et al. "The Reconstruction Toolkit (RTK), an open-source cone-beam CT reconstruction toolkit based on the Insight Toolkit (ITK)". *Journal of Physics: Conference Series* 489 (2014). DOI: 10.1088/1742-6596/489/1/012079.

[9] D. Sarrut et al. "Advanced Monte Carlo simulations of emission tomography imaging systems with GATE". *Phys Med Biol* 66.10 (2021), 10TR03. DOI: 10.1088/1361-6560/abf276.

[10] K. F. Koral et al. "SPECT dual-energy-window Compton correction: scatter multiplier required for quantification". *Journal of Nuclear Medicine* 31.1 (1990), pp. 90–98.