Frozen Hybrid Kernelised Expectation Maximisation For Bremsstrahlung SPECT Reconstruction

Background : Selective internal radiation therapy with Yttrium-90 microspheres is an e ﬀ ective therapy for liver cancer and liver metastases. Yttrium-90 is mainly a high-energy beta particle emitter. These beta particles emit Bremsstrahlung radiation during their interaction with tissue making post-therapy imaging of the radioactivity distribution feasible. Nevertheless, image quality and quantiﬁcation is di ﬃ cult due to the continuous energy spectrum which makes resolution modelling, and attenuation and scatter estimation challenging. Methods : In this study, a modiﬁed hybrid kernelised expectation maximisation is used to improve resolution and contrast and reduce noise. The iterative part of the kernel was frozen at the 72 nd sub-iteration to avoid over-ﬁtting of noise and background. A NEMA phantom with spherical inserts was used for the optimisation and validation of the algorithm, and data from 5 patients treated with Selective internal radiation therapy were used as proof of clinical relevance of the method. Results : The results suggest a maximum improvement of 56% for region of interest mean recovery coe ﬃ cient at ﬁxed coe ﬃ cient of variation and better identiﬁcation of the hot volumes in the NEMA phantom. Similar improvements were achieved with patient data, showing 47% mean value improvement over the gold standard used in hospitals. Conclusions : Such quantitative improvements could facilitate improved dosimetry calculations with SPECT when treating patients with Selective internal radiation therapy, as well as provide a more visible position of the cancerous lesions in the liver.


Introduction
Selective internal radiation therapy (SIRT) with 90 Y microspheres, also known as radioembolisation, is a recommended treatment [1][2][3] for metastatic colorectal cancer and hepatocellular carcinoma. It involves the injection of radioactive microspheres through the hepatic arteries, which flow preferentially to tumours due to their increased vasculature. The microspheres lodge in the blood capillaries of the tumour, delivering a high localised absorbed dose while minimising the dose absorbed by the healthy liver parenchyma. [4,5].
As beta particles emitted by 90 Y are deflected by atomic nuclei in tissue, Bremsstrahlung photons are produced making imaging following the therapy feasi-ble. Imaging with single photon emission computed tomography (SPECT) is used after the 90 Y treatment to qualitatively assess the activity distribution of the radionuclide [6][7][8]. Nevertheless, there are a number of challenges that degrade image quality and limit quantification. Firstly, traditional single photon imaging is based on the detection of mono-energetic photons and the images of broad-spectrum Bremsstrahlung photons have degraded spatial resolution (up to 20 mm) and artefacts [9,10]. In fact, attenuation is estimated assuming mono energetic photons which is not true for 90 Y given the wide energy spectrum [11]. Scatter correction methods utilising multiple energy windows cannot be used due to the continuous Bremsstrahlung spectrum and there is no peak that can be distinguished from the scattered photons and some studies have attempted the use of monte carlo (MC) techniques [6]. In addition, high energy photons may not be fully attenuated by collimator septa and low energy photons undergo scattering interactions with the collimator. Clinical protocols are based on previous studies which have shown that image quality can be improved by limiting the energy window to 50-250 keV and using a medium energy collimator [12][13][14][15]. Even though the optimisation of the protocol provides improvements, the partial volume effect (PVE) still makes quantification difficult. For this reason we propose and investigate a solution that reduces PVE and noise by using a synergistic reconstruction taking advantage from anatomical and functional information. The Hybrid Kernelised Expectation Maximisation (HKEM) method [16][17][18] was used to reconstruct data from a NEMA phantom and patient data with 90 Y. HKEM was derived from KEM proposed by Wang and Qi [19] to improve reconstruction of short frames by using a selection of reconstructed frames as side information, and Hutchcroft et al. [20] where the side information was from magnetic resonance imaging (MRI). In the following years different studies have shown the potential of this technique [21][22][23][24][25][26][27][28][29] and it has been used in clinical applications such as cardiovascular imaging [30,31] and cancer theranostics [32].
One potential disadvantage of HKEM is instability in terms of noise propagation of the hybrid kernel over the iterations, which can be further amplified by the lack of scatter correction. We therefore extend the HKEM algorithm to allow the iterative part of the kernel to be frozen at a chosen iteration number so that the kernel is based on an image with lower noise. We refer to the proposed technique as frozen HKEM (FHKEM). This paper is organised as follows. Section 2 describes the mathematical aspects of the hybrid kernelized reconstruction algorithm and the proposed extension. It also describes the experimental methodology, reconstruction and the image analysis. Section 3 presents results and a comparison of the different standard algorithms with the proposed one, and a discussion of these results is provided. In Section 3.3 further remarks and the limitations are described and conclusions are drawn in   3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 2M e t h o d s

Algorithm description
In the kernel method, we consider the image, λ, written as a linear combination (1) where k gj ,i st h egj th element of the kernel matrix, k, α g is the g th element of the coefficient vector that we need to estimate, and N j is the number of feature vectors used to estimate the kernel element relative to the image voxel j. Note that the usual additive term used for scatter correction is not present in this formulation since no scatter correction was performed. The FHKEM algorithm can then be described as: where 8 < : In the original HKEM, f n is always equal to the iteration number n, while for FHKEM, f n stops being updated when it reaches the selected iteration F n . α n g is the g th kernel coefficient estimated at iteration n, y i is the i th sinogram bin, L is the total number of bins, c ij is the probability that an event occurring in voxel j is detected in the i th sinogram bin.
The kernel matrix calculation is described in detail in [16], However a short description is as follows: the k gj element of the kernel can be written as follows: where is the kernel coming from the CT image and  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 is the part coming from the SPECT iterative update. The quantity x j is the coordinate of the j th voxel, z (fn) j is the feature vector that is calculated from the f th SPECT update image, and σ c , σ s , σ ds and σ ds are scaling parameters for the distances in (3) and (4). σ c and σ s control the edge preservation from the anatomical image and functional image, σ ds and σ ds control the weight of the neighbouring voxels based on the distance from the central voxel in the neighbourhood, σ ds and σ ds are in this case redundant as the images have the same dimension and voxel size.

Phantom Data
The Phantom data were acquired at the National Physical Laboratory (NPL), UK, using the Mediso AnyScan SCP. Reconstruction for this scanner has previously been implemented [33,34] in the Software for Tomographic Image Reconstruction (STIR) [35]. A NEMA phantom was scanned containing 6 spherical inserts of different volume and the same 90 Y activity concentration. The diameter of each sphere was 10 mm, 13 mm, 17 mm, 22 mm, 28 mm and 37mm and the activity at the time of scanning was 0.255± 0.001 MBq, 0.511± 0.003 MBq, 1.19± 0.01 MBq, 2.58± 0.01 MBq, 5.34± 0.03 MBq, 12.58± 0.07 MBq and the background was filled with water. The data were acquired for 2 hours, with 120 60 s projections. The energy window was set between 75 and 225 keV. A parallel-hole medium energy general purpose (MEGP) collimator was used. The CT image was acquired for attenuation estimation and was used as anatomical image in the HKEM and FHKEM.

Clinical Data
Clinical data were acquired at the Royal Surrey NHS Foundation Trust in Guildford, UK, using the GE Optima 640 SPECT/CT system for 5 patients who were treated with SIRT with 90 Y resin microspheres (SIR-Spheres). The patient data involved in this study were anonymised. The injected activity was in the range between 1.5-2.2 GBq and the SPECT acquisitions lasted 40 minutes, with 120 20 s projections. The energy window was set between 75 and 225 keV. A medium energy general purpose parallel-hole MEGP collimator was used. The CT image was acquired for attenuation estimation and was used as anatomical image in the HKEM and FHKEM.

Reconstruction Setup and Analysis
The point spread function (PSF) for the NPL Mediso gamma camera with collimator was obtained with data of multiple linear 90 Y sources at different distances from the detector. A linear fit was used to study the dependence of the sigma of the Gaussian curve on the distance from the detector. Since no PSF measurements were available for the GE Optima used for the patient data a theoretical PSF model was estimated from the collimator properties [36].
For optimisation purposes, many reconstructions were carried out with different combinations of the parameters, such as the number of subsets for all algorithms, set to 12, all the sigma parameters for the HKEM (which are described in table 1), the number of neighbours, and the iteration at which the iterative kernel is frozen. The anatomical image used as side information is a sequentially acquired 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65  Table 1 FHKEM selected parameter values Neighbours 5×5×5 functional edge σs 1 anatomical edge σc 0.3 spatial distance σ ds = σ dc 3 kernel frozen at Sub-iteration Fn 72 CT image which was manipulated to create spatial inconsistencies between CT and SPECT. In particular some of the spheres were removed from the image. The data was reconstructed using ordered subsets expectation maximisation (OSEM) with no PSF (OSEM-noPSF), OSEM, HKEM and FHKEM. For the patient data the image reconstructed with OSEM, 2 iteration using GE Xeleris was also added into the comparison. All the (F)HKEM results reported in this manuscript are obtained with the use of the same resolution modelling as the OSEM. The SPECT image size was 128⇥128⇥128, while the voxel size was 4⇥4⇥4m m 3 . The CT image size for the phantom was 512⇥512⇥82, while the voxel size was 0.98⇥0.98⇥5mm 3 . For the patient the CT image size was 512⇥512⇥161, while the voxel size was 0.98⇥0.98⇥2.5 mm 3 . The CT images were re-sampled to match the SPECT image properties. The choice of the parameter settings was based on the results in Figure 2 that is the "best" parameter value is chosen so that it provides the highest uptake while suppressing more noise.
To estimate the recovery coefficient, the mean value of the inner voxels in the biggest sphere of the OSEM image was divided by the input activity to determine a "calibration factor". In this way we can have a measure of the degradation due to PVE for each sphere. Recovery curve are usually used to correct activity values in patient images, however this is not done here as the patient data is acquired with a different scanner. The ROI analysis was carried out in terms of mean value, CoV, and contrast to noise ratio (CNR) using the ROIs shown in Figure 1(a) and 1(d). For the patient data the chosen ROI was extracted from the hottest lesion in the liver and the background (bgr) ROI from the surrounding cold liver. Preliminary investigation on the optimisation of (F)HKEM parameters in terms of ROI values and CoV was carried out. An example of this optimisation is given in Figure 2(a), which reports the trade-off between ROI mean (kBq/ml) and CoV at different values of σ ds while the other parameters are fixed, and Figure 2(b) for the choice of the sub-iteration where the kernel is frozen, F n . In Figure 2 it is possible to see the criteria of such optimisation, the "best" parameter value is chosen so that it provides the highest uptake while suppressing more noise. For example, it can be seen how σ ds = 1 voxels provides higher CoV than σ ds > 1 voxels but without increasing the ROI mean significantly. The same can be seen for F n = 600 (HKEM) and F n < 600 where the relative increase in ROI value is smaller than the relative increase of CoV. In fact, there is a difference of 12% between the CoV obtained with HKEM and FHKEM with F n = 72 while there is only 3% between the mean ROI   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63 64 65 values. Note that this plot refers to patient 4 as the differences were more visible than the phantom.

Phantom data
The comparison between the reconstruction algorithms for the NEMA phantom is shown in Figure 3, which reports the comparison in terms of the mean value in each ROI versus the CoV in the cold background using OSEM-noPSF, OSEM, HKEM and FHKEM with iterative kernel frozen at the 72 nd sub-iteration. It can be seen that for all the spheres the PSF provides a big improvement in quantification and the proposed FHKEM outperforms the other methods for all regions while providing better noise suppression with a recovery improvement of up to 56% compared to OSEM-noPSF. The recovery improvement gradually decreases with the size of the lesion. This could be related to the resolution limitation for lesions with size comparable to the kernel neighbourhood. Smaller voxel size could allow further improvement but more investigation is needed. Figure 4, shows the reconstructed images of the NEMA phantom using OSEM-noPSF, OSEM, HKEM and FHKEM with functional kernel frozen at the 72 nd sub-iteration, (comparison performed at the 50th iteration). A clear improvement can be seen when PSF is used while brighter lesions and better defined edges are evident when using HKEM. When looking at the difference between HKEM and FHKEM one can observe that FHKEM can obtain higher contrast spheres, less deformation for the smallest spheres, as well as smoother background .   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63 64 65

Patient data
When considering the application to patient data Figure 5 highlights the different behaviour of HKEM when the iterative part of the kernel continues to be updated or it is frozen at the 72 nd sub-iteration using the patient data. Figure 5(a) shows the mean ROI value (counts per voxel) of the lesion as a function of the background CoV, whereas Figure 5(b) plots the CoV against the number of iterations. The comparison between OSEM-noPSF, OSEM, HKEM and FHKEM is reported in Figure 6. In this plot the lesion mean value is plotted against the CoV in the cold liver for 4 different patients. In addition, the vertical and horizontal lines represent the values of ROI mean and CoV for the image reconstructed with the vendor software. The patient data have proven to be consistent with the phantom data. In fact, it can be seen from Figure 5 that for all the patients the ROI values obtained with HKEM and FHKEM is very close while the CoV becomes more different with higher number of iterations. For example, for patient 3 the improvement in terms of CoV, at iteration 50, of FHKEM over HKEM is 28%. Furthermore, freezing the kernel at early iterations makes the iterative algorithm more stable against noise propagation. Similarly to Figure 3, Figure 6 confirms that FHKEM outperforms all the other algorithms when it comes to reduce noise while providing activity concentration recovery higher or comparable to the other techniques. In this figure, one can also notice the quantitative improvement over the vendor reconstruction (vertical and horizontal lines). In fact, even considering very early iterations (or similar CoV) the mean ROI value can be improved up to 44% when using FHKEM.
In this work the CNR was also considered as a metric of comparison. Figure  7 compares the CNR values obtained using OSEM-noPSF, OSEM, HKEM and FHKEM frozen at the 72 nd sub-iteration for 4 different patients.The horizontal line indicates the CNR value obtained with the vendor reconstruction. This provides a confirmation of the benefits of FHKEM over the other algorithms, and in particular it shows a stabilisation of the CNR due to the "deceleration" of the noise propagation. Note that because the quantitative results for the patients are consistent  we are only showing the results for patients 1-4 for style purposes. Nonetheless, to demonstrates the consistency of the results, the images are shown for patient 5 as well.
Finally, Figure 8 shows the images reconstructed with the vendor, OSEM-noPSF, OSEM, HKEM and FHKEM frozen at the 72 nd sub-iteration for all patients. The images reconstructed with STIR are all at the 50 th full iteration (600 sub-iterations), while the vendor image was reconstructed using 2 full iterations. From a qualitative point of view severe PVE can be seen in the vendor image and OSEM-noPSF compared to the reconstructions using PSF modelling. The noise, visible in OSEM has propagated, even though reduced, through HKEM, whereas for FHKEM it is visibly reduced while showing similar contrast.
It could be argued that the comparison with the vendor software is not fair as the study uses only a single image (the one used in the clinical examination) at 2 iterations. However, the quantification comparison reports the results of the other algorithms at each iteration. To make the study more complete and fair, in Figure  9, all reconstructions images are illustrated at the iteration providing the maximum 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 CNR. This was done because comparing at the same CoV was not always possible for every patient. The figure shows still higher contrast for the proposed method, however there is no difference between HKEM and FHKEM since the maximum CNR is reached when the kernel is frozen.

Further Remarks and Limitations
This study demonstrates qualitative and quantitative improvement for Bremsstrahlung SPECT which could lead to more accurate dosimetry after SIRT with 90 Y.
Further improvement could be achieved by introducing scatter correction. Several studies have investigated MC simulation techniques for the estimation of the scattered Bremsstrahlung photons [37][38][39][40].