Singular value decomposition analysis of back projection operator of maximum likelihood expectation maximization PET image reconstruction

In emission tomography maximum likelihood expectation maximization (ML-EM) image reconstruction technique1,2 has replaced the analytical approaches (e.g. the widely used filtered back projection) in several applications, since ML-EM offers improvements in spatial resolution and stability due to the more accurate modelling of the system and to the ability of accounting for noise structure.3 In exchange ML-EM has only a linear rate of convergence2 and its computational cost is still tedious even with the rapidly increasing processing capacity of current computers. Thus a significant part of recent research activities aims at accelerating the algorithm.3,4,5 Another important property partly connected to the low convergence rate is the maximal resolution achievable for a given reconstruction method given a certain noise level towards which most of the developments are directed.6,7,8,9

In order to achieve improvement in both convergence rate and spatial resolution an ML-EM positron emission tomography (PET) reconstruction code has been developed with graphical processing unit (GPU) based Monte Carlo engine.10 GPUs parallel threads allow for running the inherently parallel neutral particle Monte Carlo transport simulations approximately hundred times faster than on a comparably priced CPU thus significantly reducing the time required for the reconstruction. As increased computational capacity allows for better physics modelling the main novelty of this code is the ability of full particle transport modelling as accurate as it is worthwhile in hope of improving image quality.11

Contrary to expectations such faithful physics modelling in the back projection step of the algorithm causes strong artefacts: modelling positron range leads to tissue dependent inhomogeneity artefacts in the reconstructed image. Furthermore, these inhomogeneities disappear when simplified Monte Carlo simulations are used without. All the differences between the two cases occur in the system matrix (derived from the Monte Carlo simulations) of the back projection step. These differences were analysed with respect to the convergence properties and stability to noise in a smaller test system by means of singular value decomposition (SVD) which is a powerful when analysing rectangular matrices. We found a significant advantage of the matrix belonging to the simplified simulations in terms of both singular values and vectors that characterized the convergence properties and stability of the algorithm. In other words more accurate physical modelling is less efficient in terms of convergence and these differences explained the perceived artefacts. However, after numerous iterations accurate modelling gives better reconstruction for low noise cases. Taking advantage of these results we created an a posteriori filtering matrix applied in each iteration after the back projection step with which we could further amplify these differences for speeding up the convergence, but without spoiling the stability to noise.

This paper is organised as follows: in the first subsection of Materials and methods the details of our reconstruction code and a simplified system are described. Second subsection contains the notations. Third subsection gives a short tutorial for SVD including the closely related Picard condition and the corresponding convergence and noise analysis. Results section is divided into four subsections, which present the perceived artefact in detail, the use of faithful modelling and our newly developed SVD filtering method with comparison to the original algorithm. Finally, possible solutions for the implementation are offered with theoretical considerations to further research. Discussion summarises the impact of the SVD analysis and filter and presents our connecting further research goals.

A 3D Monte Carlo based ML-EM image reconstruction code named PANNI (PET Aimed Novel Nuclear Imager) has been developed in the framework of the TeraTomo project.12,13

PANNI is a Monte Carlo based image reconstruction software written for GPUs using C and CUDA environments surmising roughly 40 000 lines.11 The key feature of our software is the possibility of faithful Monte Carlo modelling which accounts for positron range, gamma photon-matter interaction and detector response supported by advanced variance reduction methods. Detectors around the object are positioned on a quasi-cylindrical surface with dodecagon cross section. Detector response is either simulated or a pre-generated tabulated response function may be used. Positron range modelling simplifies to the following probability density function.14

℘(r)=aA2re−Ar+bB2re−Br$$\begin{array}{} \wp(r)=aA^2re^{-Ar}+bB^2re^{-Br} \end{array} $$

with r being the positron range distance, a, A, b and B material dependent constants.

As sampling each of the terms is equivalent to sampling the sum of two exponentially distributed random variables x can be obtained by using double exponential sampling.15 Advanced variance reduction methods are implemented for source angular sampling outgoing direction and energy biasing and for free flight sampling. The Monte Carlo engine has been validated against MCNP5.16 The code is capable of simulating 10⁸ photon pairs per second on a commercially available GPU (NVidia GeForce 690).

Both the forward projection and the back projection steps are carried out via the Monte Carlo method. In the back projection step some of the physics modelling may be turned off.

The code has been tested with two geometries, a sophisticated scanner geometry (“full system”) and a simplified smaller system (“1D model”). Acquisition geometry for the full system can be set as wished, in our current setup it consists of a dodecagon with inscribed circle radius of 8.7 cm, packed on each side with an array of 39 × 81 LYSO detector pixels of 1.17 mm sided squares comparable to a small animal PET scanner similar to the Mediso nanoScan PET/CT scanner. Coincidence counting is accepted between detector pixels on opposite and next to opposite dodecagon sides (1:3 coincidence). The voxel space of the full system is divided into 128×128×128 voxels (0.3 mm sided) and contains a water-cylinder (light grey area in Figure 1 and Figure 2) except for a smaller cylindrical area containing bone material (dark grey area in Figure 1 and Figure 2) Activity phantom for the evaluation is a cylindrical ring of ¹⁵O partially located in bone material (the more commonly used ¹⁸F gives less conspicuous results). From now on simplified modelling means the neglect of the positron range effect in the back projection Monte Carlo simulations in contrast with faithful modelling which accounts for positron range.

Figure 1

Figure 2

The voxel space of the 1D model is a 38.4 mm long interval containing 256 voxels half of which is located in bone material the other half in water. Detector pixels are assigned in two parallel sections each containing 81 crystals of a size 1.17.×.1.17 mm. Every pixel is in coincidence with every pixel on the opposite side. Roughly speaking the 1D model is a cross-section of the full system geometry of PANNI (Figure 3).

Figure 3

The 1D model contains only positron range modelling, neither gamma photon-matter interaction nor detector response modelling is included. Detection is based on the angle of view of the detector from a given voxel. The two analysed settings are: positron range neglected (back projection posrange OFF) and positron range modelled (back projection posrange ON) in the back projection. Forward projection always accounts for positron range.

We have compared two reconstruction results for the sophisticated scanner geometry: one with full physics modelling in the back projection (Figure 1) and one omitting positron range in the back projection (Figure 2). After 80 iterations faithful modelling produced the reconstruction in Figure 1. The cross section of the cylindrical ring phantom in radial direction is originally a box function which is blurred due to gridding and averaging in a given voxel (similarly to partial volume effect). Therefore, it can be well approximated by a gaussian with a Full Width at Half Maximum of 3 voxels which is indicated by the red line in Figure 1 and Figure 2. The Full Width at Half Maximum/Full Width at Tenth Maximum is calculated by fitting a gaussian in radial direction along the ring separately for each angular position with resolution of 1 degree.

Accounting for positron range in the back projection caused systematic inhomogeneity in the reconstructed image in Figure 1. The activity estimation in the bone material is appropriate (FWHM = 3.5 voxel) in contrast with the activity of the voxels located in water which is underestimated (FWHM = 5 voxels). The inhomogeneity reduces with simplified back projection, when positron range is neglected in the Monte Carlo simulation. Also the FWHM is reduced in the water area (Figure 2).

The effect of modelling any physical phenomenon appears in the system matrix, thus our analysis aims at finding the differences in the back projection system matrix caused by positron range. Due to its size the system matrix of the full system cannot be stored thus the 1D model was used for further calculations. For real and valid results a comparison of the test system with the full system is needed. Simulations confirmed the analogous behaviour as the perceived artefact appeared in the case of positron range modelling in the back projection while neglecting positron range abolished the problem. Thus, the analysis focuses on the positron range effect, back projection posrange OFF and back projection posrange ON settings are compared. (The positron range free path is sampled with two random variables similarly to the full system). As the corresponding system matrix is of a size 6561x256 it can be directly stored and also the numerical SVD calculation may be carried out.

Comparing the back projection posrange OFF and back projection posrange ON case in terms of singular values of the system matrix, Figure 4 shows the positive difference for the first 133 index belonging to the former setting.

Figure 4

In the light of the convergence analysis of Ref.17 smaller singular values mean that the corresponding frequency components of the solution are later reconstructed with positron range modelling compared to the positron range neglecting back projection. This space frequency characterises the singular vectors of the voxel space (Figure 5), which is a second but not less significant difference between the two types of system matrices.

Figure 5

Figure 5 on the left accounts only for the symmetries of the system while on the right reflects also the tissue map of the volume. As the average positron free path is much longer in water than in the bone material space frequency of every basis vector is smaller in the area located in the water. Thus the reconstruction of the activity of these voxels is significantly slower and this property is the reason of the obtained artefact resulted from faithful modelling in the back projection and the solution to the perceived anomalous behaviour.

The third and final difference occurs in the sinogram space basis vectors with which the measurement-forward projection Hadamard ratio can be unfolded in a given iteration. The absolute value of the obtained spectral coefficients can be seen in Figure 6 after 15 iterations for the positron range neglecting and modelling case respectively (the spectrum varies slowly through iterations).

Figure 6

Faster decay in the spectral coefficients equals to heavier blurring in the back projection.18 This means that the positron range modelling gathers less information from the Hadamard ratio in a given iteration than the positron range neglecting back projection. This also implies the faster convergence and higher stability to noise.

Advantage of faithful modelling

ML-EM with faithful modelling converges to the exact solution.1,2 Algorithm using simplified back projection converges to different fix point due to unmatched forward-backward projector pair and only approximates the solution.2,3,5 These statements were verified by test simulations with varying noise level (Figure 7). In presence of high noise semi-convergence property of the algorithm3,19 is dominant, reaching the optimum estimate as fast as possible is desired which is the advantageous property of the simplified back projection. In low noise case after numerous iterations the faithful modelling outperforms the simplified one reaching much better activity estimate (Figure 7). SVD analysis showed the disadvantage of the former in the term of convergence speed but due to matched forward-backward projector pair it converges to the exact solution1 in contrast with the simplified modelling.

Figure 7

This result resolves the contradiction as additional information indeed leads to better image reconstruction thus the perceived anomaly was only apparent. Simplification just luckily affects the behaviour of the algorithm from a mathematical point of view. However, this form is not the ideal back projection operator but the one that is easy to implement without much modification to the original algorithm. To obtain a better back projection operator the previously listed advantages can be amplified with a posteriori manipulation and a close to ideal form can be reached. As a possible degree of freedom, singular values of the simplified back projection operator can be further increased similarly to the accompanying effect of the simplified modelling. In this case U and V matrices are unchanged. The possible fix points of the algorithm can be obtained from the next equations (ratio is in a Hadamard sense) as the update process multiplies (also in Hadamard sense) the current estimate by 1 when the following is true:

ATyr|A|=ATyrAT∗1⋮1Lor−space=1⋮1voxelspace$$\begin{array}{} \displaystyle \frac{{\boldsymbol A}^{\boldsymbol T} y_r}{|{\boldsymbol A}|}=\frac{{\boldsymbol A}^{\boldsymbol T}y_r}{{\boldsymbol A}^{\boldsymbol T}*\begin{pmatrix} {1}\\{\vdots}\\1\end{pmatrix}_{Lor-space}}=\begin{pmatrix} {1}\\{\vdots}\\1\end{pmatrix}_{voxel \,space} \end{array} $$

Rearranging:

ATyr−1⋮1ATy~=0$$\begin{array}{} {\boldsymbol A}^{\boldsymbol T}\left(y_r-\begin{pmatrix} {1}\\{\vdots}\\1\end{pmatrix}\right){\boldsymbol A}^{\boldsymbol T}\tilde{y}=0 \end{array} $$

Using the dyadic definition of SVD:

ATy~∑r=1Rank(A)σrvrurT∗y~=0$$\begin{array}{} \displaystyle {\boldsymbol A}^{\boldsymbol T}\tilde{y}\left(\sum\limits_{r=1}^{Rank({\boldsymbol A})}\sigma_rv_ru_r^T\right)*\tilde{y}=0 \end{array} $$

v_r and u_r are the columns of matrix V and U respectively. As V is an orthogonal matrix the linear combination above equals to zero precisely when every σrurTy~$\begin{array}{} \sigma_ru_r^T\tilde{y} \end{array} $

coefficient equals to zero. The singular values are all nonzero thus the urTy~$\begin{array}{} u_r^T\tilde{y} \end{array} $

dot product equals to zero for ∀r. This implies that UTy~=0.$\begin{array}{} {\boldsymbol U}^{\boldsymbol T}\tilde{y}=0. \end{array} $ Matrix U remains the same with singular value modification so the possible fix points are unchanged.

SVD filter

The easiest way to modify the singular value spectrum of the back projection system matrix is the application of a matrix of the form:

B=VD∗VT$$\begin{array}{} {\boldsymbol B}={\boldsymbol V}{\boldsymbol D}^*{\boldsymbol V}^{\boldsymbol T} \end{array} $$

D^∗ is the matrix with which D^∗D has the desired form. In noiseless (test) case, i.e. when there is no noise added to simulations, this form is (the scalar multiple of) the identity matrix, in agreement with the convergence analysis5,17,20, as the singular values are clustered together as far as possible. The SVD filter fastens the convergence with two orders of magnitude. However, it cannot be applied straightforward for the real, noisy case.

The measurement process equals to Ax = UDV^Tx where x stands for the activity distribution. Thus, the measurement attenuates its frequency components according to the singular value spectrum (multiplication with matrix D) and adds some noise to the result. As so only those components which fit to the discrete Picard condition (Figure 6) can be amplified.

We performed several reconstructions with such an SVD filter. The best result was obtained when D^∗ diagonal matrix contains elements of the form (Figure 8):

1σipifi<cut−off0ifi≥cut−offfor some 0<p=1−ε and small ε.$$\begin{array}{} \left\{\begin{array}{} \frac{1}{\sigma_i^p}&if\,i\,\lt cut-off\\ 0&if\,i\geq cut-off \end{array}\right.\text{for some 0}\lt p=1-\varepsilon \text{ and small }\varepsilon. \end{array} $$

for some and small.

Then D^∗D is also diagonal with the following entries: σi1−pifi<cut−off0ifi≥cut−off$\begin{array}{} \left\{\begin{array}{} \sigma_i^{1-p}&if\,i\,\lt cut-off\\ 0&if\,i\geq cut-off \end{array}\right. \end{array} $ which are closer to 1, so the singular value spectrum of the resulted back projection system matrix is contracted.

The conventional ML-EM formula looks like as follows (y_r is the pointwise, i.e. Hadamard ratio vector of measured and forward-projected data, yr=ymAx.$\begin{array}{} y_r=\frac{y_m}{{\boldsymbol A}x}. \end{array} $

1 is a sinogram space vector containing ones in every coordinate):

xn+1=xnATyrAT1=xnVDUTyrVDTUT1$$\begin{array}{} \displaystyle x^{n+1}=x^n\frac{{\boldsymbol A}^{\boldsymbol T}y_r}{\boldsymbol A^{\boldsymbol T }{\bf 1}}=x^n\frac{{\boldsymbol V}{{\boldsymbol D}}{\boldsymbol U}^{\boldsymbol T}y_r}{{\boldsymbol V}{\boldsymbol D}^{{\boldsymbol T}}{\boldsymbol U}^{{\boldsymbol T}}{\bf 1}} \end{array} $$

Instead, we use the modified iteration formula which looks like as follows (B^T = B as being symmetric and V^TV = I as V is orthogonal. I is the identity matrix. The ratio and the multiplication in the update process of xⁿ is in Hadamard sense):

xn+1=xnBATyrATB1=xnVD∗VTVDUTyrVD∗TVTVDTUT1=xnVD∗DUTyrVD∗TDTUT1$$\begin{array}{} \displaystyle x^{n+1}=x^n\frac{{\boldsymbol B}{\boldsymbol A}^{{\boldsymbol T}}y_r}{{\boldsymbol A}^{\boldsymbol T}{\boldsymbol B}{\bf 1}}=x^n\frac{{\boldsymbol V}{\boldsymbol D}^*{\boldsymbol V}^{\boldsymbol T}{\boldsymbol V}{\boldsymbol D}{\boldsymbol U}^{\boldsymbol T}y_r}{{\boldsymbol V}{\boldsymbol D}^{*{\boldsymbol T}}{\boldsymbol V}^{\boldsymbol T}{\boldsymbol V}{\boldsymbol D}^{\boldsymbol T}{\boldsymbol U}^{\boldsymbol T}{\bf 1}}=x^n\frac{{\boldsymbol V}{\boldsymbol D}^*{\boldsymbol D}{\boldsymbol U}^{\boldsymbol T}y_r}{{\boldsymbol V}{\boldsymbol D}^{*{\boldsymbol T}}{\boldsymbol D}^{\boldsymbol T}{\boldsymbol U}^{\boldsymbol T}{\bf 1}} \end{array} $$

Figure 8 shows the result of the reconstruction compared to regular ML-EM algorithm with both back projection posrange ON and back projection posrange OFF settings (6 × 10⁵ positron used). The L₂-norm of the difference of the activity distribution and the current activity estimate is presented for each setting after a given number of iterations. Smaller L₂-norm means better agreement.

Figure 8

SVD filter outperforms the best setting so far as the initial convergence is faster and better agreement is reached in every iteration. Furthermore, the rise of the discrepancy due to semi-convergence occurs later. Additionally, the faster initial convergence of the positron range neglecting back projection (back projection posrange OFF) can be seen compared to back projection posrange ON.

In our paper all of the three SVD matrices from the factorisation of the system matrix were analysed. The results explained the perceived artefact as the convergence speed of the scheme with positron range modelling back projection was material dependent and the further advantage of the simplified back projection in terms of overall convergence speed and stability to noise. The presented SVD filter further amplified these favourable properties so as to fasten the algorithm but preserving its robustness. Additionally the use of faithful modelling was pointed out when high computational capacity is available.

A posteriori filtering is in close relation with such a priori methods when a regularizing term is added to the likelihood function in the problem formulation for decreasing noise sensitivity and accelerate the convergence. The resulted filtering term is then present usually in the nominator of the backprojection in additive form and arises from known constraints about the imaging process. In general, some kind of regularisation is almost always required due to the ill-posed nature of the reconstruction problem. Our SVD filter approaches from a bit different point of view: it does not require any a priori knowledge. The effect of B matrix is not strictly regularisation but rather deconvolution which accelerates the convergence and the deconvolution process itself has to be regularised due to the presence of noise. Thus, the process has a filtering effect as well.

Monte Carlo simulations result slightly different system matrix elements through iterations which imply slightly different SVDs and B matrices. So as to make the most of the presented SVD filtering technique our further aim is to find the connection with conventional deconvolution methods which obtain exactly the same effect but the filtering factors can be recalculated in each iteration tailored to the given back projection system matrix. The special form (PSWF which is strongly connected to Fourier-transform22,23) of the singular vectors makes this direction promising.