In emission tomography maximum likelihood expectation maximization (ML-EM) image reconstruction technique1,2 has replaced the analytical approaches (
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
This paper is organised as follows: in the first subsection of
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
with
As sampling each of the terms is equivalent to sampling the sum of two exponentially distributed random variables
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 (“
Top view of the reconstruction of the cylinder-ring mathematical phantom of the Figure 1
Top view of the reconstruction of the cylinder-ring mathematical phantom of the Figure 2
The voxel space of the
Mathematical phantom and system geometry for Figure 3
The
The effect of positron range modelling and the average positron range is accounted for by the system matrix. As an obvious tool for the analysis of rectangular matrices, the algorithm and the back projection step is examined by means of singular value decomposition. SVD is a factorisation of any
SVD was used for the analysis of convergence speed of the reconstruction algorithm with respect to the applied back projection. According to Liu
Sinogram space singular vectors can measure the information content of a given measurement-forward projection Hadamard ratio in the corresponding back projection step by means of Picard condition formalism. For the existence of a square integrable solution to the problem
In case of matrices instead of integral operators, the discrete Picard condition requires the spectral coefficients
to decay faster in average than the singular values.18 Despite back projection is not a direct inversion from this point of view the faster is the decay of the spectral coefficients
The aforementioned difference in the sinogram space basis affects the
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
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.
Singular values of the system matrix for positron range neglecting and modelling case. Increased values for the former imply the faster convergence of the corresponding (first 133) basis components of the activity estimate.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.
One of the voxel space singular vectors of the system matrices corresponding to positron range neglect (left – back projection posrange OFF) and positron range modelling (right – back projection posrange ON). Back projection posrange OFF reflects only the symmetries of the geometry. Back projection posrange ON accounts for the material map as well, increased position uncertainty due to longer average positron free path implies smaller space-frequency in water area.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).
Absolute value of the spectral coefficients of the measurement-forward projection Hadamard ratio in the sinogram basis corresponding to positron range neglect (left – back projection posrange OFF) and positron range modelling (right – back projection posrange ON). Faster decay means less information gathered as the coefficients of the horizontal plateau are corrupted by noise thus it represents an error level estimate. Due to one to one correspondence property of SVD between sinogram and voxel space singular (basis) vectors these basis coefficients of the activity are not hoped to be correctly estimatedFigure 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.
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.
The L2-norm of the difference between the activity distribution and the current estimate after a given number of iterations. Smaller value means better agreement. Subfigure on the left shows the result of the noiseless test case where the convergence of back projection posrange ON setting to the exact solution and the convergence of back projection posrange OFF setting to an other fix point is presented. Subfigure on the right shows the result of a simulated reconstruction with 106 positron used for measurement generation and in both forward and back projection Monte Carlo simulations. After slower initial convergence back projection posrange ON reaches much better activity estimate. Back projection posrange OFF converges to a different fix point similarly to noiseless case.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
Rearranging:
Using the dyadic definition of SVD:
coefficient equals to zero. The singular values are all nonzero thus the
dot product equals to zero for
The easiest way to modify the singular value spectrum of the back projection system matrix is the application of a matrix of the form:
The measurement process equals to
We performed several reconstructions with such an SVD filter. The best result was obtained when
for some and small.
Then
The conventional ML-EM formula looks like as follows (
Instead, we use the modified iteration formula which looks like as follows (
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 × 105 positron used). The L2-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 L2-norm means better agreement.
The L2-norm of the difference between the activity distribution and the current estimate after a given number of iterations. Smaller value means better agreement. Reconstruction with SVD filter outperforms the best setting so far in terms of faster initial convergence and the farther starting point of increasing discrepancy due to semi-convergence. Also the faster initial convergence of positron range neglecting back projection can be seen compared to positron range modelling back projection.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.
SVD filter requires the calculation of
Observing (the square of) the singular value spectra of the system matrix, low-pass filtering is not ideal but PSWFs are good approximation for the singular vectors. Owing to this favourable property, matrix B can be calculated if singular value spectrum has been obtained (
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.