Simultaneous estimation of the emission and attenuation maps in TOF-PET

1. Simultaneous estimation of the attenuation and emission in TOF-PETJoint work VUB (M Defrise), KULeuven (J. Nuyts, A. Rezaei) Herhaling titel van presentatie

2. Simultaneous estimation of the emission and attenuation maps in TOF-PET PET-CT PET data : gg coincidences CT data: x-ray attenuation Reconstruction of the organ density Reconstruction with correction for the g attenuation "m map" "Emission map": tracer's biodistribution 2

3. Simultaneous estimation of the emission and attenuation maps in TOF-PET PET-MRI PET data : gg coincidences MR data: H+ density, T1,T2,.. segmentation + a priori m assignments No direct link with attenuation coefficients ? Reconstruction with correction for the g attenuation "m map" ? "Emission map": tracer's biodistribution e.g.: V. Keereman, Y. Fierens, C. Vanhove, T. Lahoutte, S. Vandenberghe, 2011 MRI-based attenuation correction for micro-SPECT, Mol. Imag. Hofmann M et al, Eur J Nucl Mol Imag 2009, 3

4. Simultaneous estimation of the emission and attenuation maps in TOF-PET PET without additional data ? PET data : gg coincidences ? "m map" "Emission map" - Welch A, Campbell C, Clackdoyle R, Natterer F, Hudson M, Bromiley A, Mikecz P, Chillcot F, Dodd M, Hopwood P, Craib S, Gullberg G T and Sharp P 1998 Attenuation correction in PET using consistency information, IEEE Trans. Nuclear Science -Nuyts J, Dupont P, Stroobants S, Benninck R, Mortelmans L, Suetens P 1999 Simultaneous maximum a posteriori reconstruction of attenuation and activity distributions from emission sinograms, IEEE Trans Med Imag: MLAA - Papers by A. Bronnikov 4

5. h(l - c /2) Time-of-flight PET in a nutshell Anger 1966. 1980’s: Allemand et al, Ter Pogossian et al, Mullani et al, Lewellen et al detector d2, time t2 Measure the arrival time difference t= c (t2 – t1 )/2 localize decay along LOR with uncertainty l = c /2  FWHM = 500 ps l FWHM = 75 mm o e+ l detector d1, time t1 “Each coincident event yields more information than without time-of-flight, by a factor ≈ patient diameter/lFWHM ” Snyder et al 1980’s, Tomitani 1981, Watson 2007, Surti et al 2006/8, Popescu & Lewitt 2006, …

6. h(l - c /2) Time-of-flight PET: analytical data model detector d2, time t2 o e+ l detector d1, time t1 corrected data attenuation factor activity image TOF profile measured data for all lines L and all t attenuation image

7. Simultaneous estimation of the emission and attenuation maps in TOF-PET TOF-PET without additional data ? TOF-PET data : gg coincidences "m map" "Emission map" M. Defrise, A. Rezaei, J. Nuyts, TOF-PET data determine the attenuation sinogram up to a constant, Phys Med Biol. 2012 - Proof based on the consistency conditions for TOF-PET data (Continuous and Discrete Data Rebinning in TOF-PET,IEEE Trans Med Imag 2008) - Analytic algorithm to determine the m map - Not as stable as a CT scan but benefit: m map is matched to emission map - Iterative approaches: MLAA and MLACF 7

8. - Analytic solution for the 2D problem- Some results- Analytic solution for the 3D problem- Discrete methods- Some results Herhaling titel van presentatie

9. Analytic solution for the 2D problem Herhaling titel van presentatie

10. The non-attenuated 2D TOF PET data for a gaussian TOF profile satisfy the local consistency condition Proof: insert the definition of p and use

11. Some insight of the consistency condition The most likely annihilation point (MLA) of a TOF-LOR is the point x* at the maximum of the TOF profile, i.e at l=t t' t In the limit s=0 the equation means that all TOF-LORs sharing the same MLA x* must be equal (or: the characteristic curves of this PDE in the 3D space (f,s,t) are loci of constant MLA).

12. The attenuated 2D TOF PET data are Vocabulary: a is the attenuation factor m is the attenuation coefficient The proof of uniqueness and the algorithm exploit: p=m/a must be consistent a is independent of the TOF variable t

13. Since a > 0, we have For each LOR (f,s) for which m(f,s,t)>0, applying a least-square fit w.r.t. the TOF variable t and with only two unknown parameters yields

14. The solution to the least-square fit is with quantities directly calculated from the measured data: When does this work ? Answer: For all LORs such that - m > 0 - along the LOR, the activity is not restricted to a point source

15. Procedure: - apply above method to estimate in the interior of supp(Rf) - integrate the gradient within supp(Rf) to find Rm up to a global additive constant, hence a=exp(-Rm) up to a global factor. (Least-square solution found by solving a Poisson equation, we use a Landweber iteration) - estimate the constant (various methods, e.g. add small external object with known activity or attenuation) - either directly use the estimated a - or reconstruct m

16. Some details on the implementation (slide J. Nuyts) The noisy TOF sinogram is first smoothed in 3D. That seems essential, without this the noise propagation is not acceptable. A moderate smoothing seems enough: 2.5 pixels radially and angularly. In TOF-direction I smoothed with a Gaussian with a width of 70% of the TOF-kernel. This addition smoothing in the TOF-direction was taken into account by increasing s in the analytical expressoin. In a real PET system, there will be such smoothing along the TOF-direction anyway by binning into a relatively small number of bins. In this test program, I stored the TOF in 128 bins, avoiding problems due to binning. Straightforward implementation of the analytical inversion then yields derivative sinograms with extreme values near the edges, but with only moderate noise propagation inside the object. To find a region with reliable values, a sinogram was computed by summing the original noisy sinogram over all TOF-bins (nonTOFsino). This sinogram was thresholded, the resulting binary image was strongly eroded to yield a (overly conservative) mask. Values inside that mask should definitely be subject to moderate noise only. For the radial derivative sinogram, the maximum value inside that mask is computed. All pixels in the derivative sinogram exceeding this maximum are labeled. Also their horizontal neighbors are labeled. In the subsequent Landweber-like iterations, those pixels will be excluded from the computations. Exactly the same thing was done for the angular derivative sinogram. Then the Landweber-like iterations are applied and a calculated sinogram is obtained. As can be seen in the images, it looks surprisingly good considering the noise in the derivative sinograms. I assume this is because this iterative procedure basically integrates in two directions. The resulting calculated sinogram should be OK except for a constant. The calculated sinogram values are typically lower than the true ones, because the iterations were started from a zero image. Therefore, I expect that a positive constant must be added. To estimate a lower limit of that constant, I compute a temporary sinogram: I take the median using a 7x7 mask to exclude extreme values, and I use nonTOFsino to set the background to zero. Then I compute the minimum of that temporary sinogram. If that minimum is negative (which was the case for the few simulations I did) , I subtract it (thus adding a positive offset) from the calculated sinogram. Then the resulting calculated sinogram is reconstructed with FBP. It has terrible streaks in the background, but the part inside the object looks rather good. Because my offset was conservative, the reconstructed attenuation is still a bit low. But if this image can be segmented and if the attenuation is known in a region, the values can be easily corrected (see below). For comparison, I computed a noise-free blank scan from the known activity image (forward projection without attenuation), and using nonTOFsino as the transmission scan, an MLTR and FBP reconstruction were calculated. They look better of course, but considering that using a noise-free blank scan is an obvious advantage, I think this indicates that the analytical inversion is rather stable. In this simulation, the mean TOF-integrated count was 23.7. For a 3D PET system with 5 or 7 segments, that would correspond to a TOF-integrated count of about 4...5 per segment.

17. true activity true attenuation Slide J. Nuyts. Attenuated emission sinogram, summed over all 128 TOF-bins. Maximum = 106, mean = 23.7. TOF-resolution was 7.5 cm = 24 pixels FWHM. The simulation was done oversampling the image pixels with 3x3 samples/pixel, and with 3 subsamples per detector. 128 x 128 pixels, pixel size = 3.125 mm To suppress the noise, the TOF-sinogram was smoothed with a Gaussian in 3D, using the following FWHM: 2.5 pixels radially, 2.5 pixels along the angle, and 0.7 x TOF-FWHM = 16.8 pixels in TOF-direction. In the reconstruction, s was set to the combined effect of the TOF-kernel and the additional smoothing along the TOF-dimension. Computed radial derivative (left), true radial derivative (right). Computed angular derivative (left), true angular derivative (right). Computed attenuation sinogram (left), true attenuation sinogram (right). Defrise-recon zeroed outside boundary MLTR from true activity image and noisy sinogram FBP from true activity image and noisy sinogram true attenuation Defrise-recon.

18. Simple procedure to find the offset. As mentioned above, the reconstruction from the calculated sinogram is not quantitatively accurate, because the analytical inversion is only accurate except for a constant in the sinogram. Quantification can be restored, if the attenuation coefficient in a part of the image is known. Here is a possible way to do this, based on the linearity of FBP. Let R be the original non-quantitative reconstruction. Make a non-TOF sinogram that is unity inside the object. Here, this was done by thresholding the TOF-integrated sinogram nonTOFsino, setting all non-zero values to 1. Reconstruct this sinogram, generating a reconstruction image corresponding to a sinogram offset of 1. We call it U. Segment the original reconstruction R, identifying a region where the attenuation coefficient is known a priori. We call this known attenuation value m. Compute for that region the mean value in R, called MR, and the mean value in U, called MU. Correct the image R as follows: Rnew = R + (m – MR) / MU * U For the 2D simulation, the segmentation was done by thresholding R with a value of mtissue / 2. Because the quantification of R is not that bad, this separates tissue and bone from the lungs. We also discard the background by assuming zero attenuation where the sinogram count is zero (convex hull). Thus, the segmentation contains bone and tissue. We compute MR as the median value of the selected pixels in R. We compute MU as the median of the corresponding pixels in U. The image was then corrected using the known attenuation of tissue. The result is shown below. Result of the correction applied to the simulation shown above. True image Reconstruction R Corrected image Rnew True image Reconstruction R Corrected image Rnew Result of the correction applied to a simulation with three times less counts (mean of 7 counts/non-TOF pixel), obtained with the same parameters.

19. True image Reconstruction R Corrected image Rnew Result of the correction applied to a noise-free simulation, using the same parameters. This illustrates the impact of the smoothing, which is needed to deal with noisy data.

20. Analytic solution for the 3D problem

21. The non-attenuated 3D TOF PET data for a gaussian TOF profile satisfy two consistency conditions The algorithm exploit: D(m/a)=K(m/a)=0 a is independent of t estimation of the 4D gradient PDEs: Defrise et al TMI 2008. Attenuation problem: submitted to IEEE MIC 2012, A. Rezaei et al

22. A discrete approach

23. Notations

24. Two possible approaches: I. MLAA: Poisson maximum-likelihood for l,m (Nuyts et al, Trans Med Imag 1999; for TOF: Rezaei et al IEEE Med Imag Conf 2011) II. MLACF: Poisson maximum-likelihood for l,a (Nuyts et al, submitted IEEE MIC 2012) + simple, fast, and works well on simulated data (see data below) - in 3D there are more unknowns a than m (more LORs than voxels)

25. is easily maximized w.r.t.attenuation coefficients a at fixed emission image l: we neglect the non-negativity constraint on m and use always y/p

26. We are left with the problem of maximizing An iterative algorithm based on a surrogate function: MLACF The MLACF algorithm is monotone Convergence is not guaranteed (since L is not concave) Generalizations for scatter background OK and implemented (Nuyts et al submitted MIC 2012, Panin et al submitted MIC 2012)

27. Examples of results with simulated 2D data

28. Estimating the sinogram of the attenuation map from TOF data. Results with the MLACF method. Large ellipse: 300 mm x 480 mm. Background 1, lungs 0.2, tumors 2, 3 and 4 Attenuation coef. Background 0.00966/mm, lungs 0.00266/mm, spine 0.01866/mm Most attenuated LOR: X 0.03. 32 TOF bins, FWHM = 50 mm, TOF bin spacing is 16 mm 256x256 emission sinogram (non attenuated, non-TOF) 256x256 mu map sinogram Simulated 256x256 mu map

29. Noise-free data MLACF 4 subsets 16 iterations TOF FBP w/o (top) and with (bottom) attenuation correction using exact mumap Hamming window

32. Noise-free data True mu map (transmission/blank) Grey scale 0 to 0.25 MLACF 4 subsets 16 iterations Grey scale 0 to 0.25

33. Noise-free data MLACF uses 4 subsets 16 iterations Note that the scaling of the estimated acfs is obtained by requiring the maximum estimated acf to be 1 at each iteration of the algorithm.

34. Noisy data Maximum count is 100, total count is 15 millions MLACF uses 4 subsets. FBP is the usual TOF FBP with a hamming window TOF FBP with attenuation corrections assuming known mumap MLACF 4 subsets 8 iterations (top), 16 iterations (bottom)

35. Noisy data Maximum count is 100, total count is 15 millions MLACF uses 4 subsets. True mu map (transmission/blank) Grey scale 0 to 0.25 MLACF 4 subsets 12 iterations Grey scale 0 to 0.25

36. Noisy data Maximum count is 100, total count is 15 millions MLACF uses 4 subsets 12 iterations Note that the scaling of the estimated acfs is obtained by requiring the maximum estimated acf to be 1 at each iteration of the algorithm.

37. Noisy data Maximum count is 20, total count is 3 millions MLACF uses 4 subsets. FBP is the usual TOF FBP with a hamming window TOF FBP and with attenuation correction assuming exact mumap MLACF 4 subsets 8 iterations

39. Estimating the sinogram of the attenuation map from TOF data. Results with attenuating medium outside the emission Large ellipse: 300 mm x 480 mm. Background 1, lungs 0.2, tumors 2, 3 and 4 Attenuation coef. Background 0.00966/mm, lungs 0.00266/mm, spine 0.01866/mm Bed: 0.01/mm 32 TOF bins, FWHM = 50 mm, TOF bin spacing is 16 mm Simulated 256x256 emission map Simulated 256x256 mu map includes an (uncomfortable) bed

40. Noise-free data (with bed attenuation) True mu map (transmission/blank) Grey scale 0 to 0.35 MLACF 4 subsets 16 iterations Grey scale 0 to 0.35 No rescaling of the mu map !!

41. Noisy data (with bed attenuation) max LOR count = 10, total count 1.26 million TOF FBP and with attenuation correction assuming exact mumap Grey scale [0,2] MLACF 4 subsets 6 iterations Grey scale [0, 2.26] MLACF 4 subsets 6 iterations Grey scale [0, 0.35]