Maximum-likelihood joint image reconstruction and motion estimation with misaligned attenuation in TOF-PET/CT

This work is an extension of our recent work on joint activity reconstruction/motion estimation (JRM) from positron emission tomography (PET) data. We performed JRM by maximization of the penalized log-likelihood in which the probabilistic model assumes that the same motion field affects both the activity distribution and the attenuation map. Our previous results showed that JRM can successfully reconstruct the activity distribution when the attenuation map is misaligned with the PET data, but converges slowly due to the significant cross-talk in the likelihood. In this paper, we utilize time-of-flight PET for JRM and demonstrate that the convergence speed is significantly improved compared to JRM with conventional PET data.


Introduction
Patient motion in positron emission tomography (PET) is a source of error due to possible mismatches between the PET data and the computed tomography (CT) attenuation map (μ-map) (Nyflot et al 2015). In Bousse et al (2016) we developed a motion compensated reconstruction scheme for gated PET data, namely joint reconstruction/motion estimation (JRM), to jointly estimate both the activity distribution and the motion field, by penalized likelihood maximization. Unlike previous works in this field (Jacobson andFessler 2003, Blume et al 2010), our model assumed that both the activity distribution and the μ-map are affected by motion. This model led to the following result: the JRM-reconstructed PET gates are the same for any input μ-maps derived from deformations of a common μ-map. This is because the estimated motion automatically accounts for PET/μ-map misalignment. However, in case of large mismatches, our results in Bousse et al (2016) showed that JRM needs a high number of iterations due to the significant cross-talk in the joint-likelihood.
Recent work Rezaei et al (2012), (2014) demonstrated the potential of jointly reconstructing the activity distribution and the μ-map from time-of-flight (TOF) PET data. Since JRM with warped μ-map presents some similarities with joint activity/attenuation reconstruction from emission data, in the sense that the warped μ-map must match the PET projections at each gate, TOF PET can increase the JRM convergence rate-especially in situations where the input μ-map is misaligned with the PET data.
In this work we investigate the ability of JRM with TOF PET to deal with a misaligned μ-map. Section 2 presents the JRM optimization problem for TOF PET and non-TOF PET, and summarizes the algorithm. In section 3 we compare the convergence rate of JRM with TOF PET and non-TOF PET on a simulated end-expiration PET gate and a simulated endinhalation μ-map. Results are discussed in section 4.

TOF maximum-likelihood for joint image reconstruction and motion estimation
In this section we use notations similar to Bousse et al (2016, sections II and III).
2.1.1. Motion-free model. The activity distribution and the attenuation map respectively take the form of 2 functions ∈ + C f and µ ∈ + C , where + C denotes the set of non-negative continuous functions defined on R 3 . The attenuation map μ is reconstructed independently from separate measurements such as x-ray CT or segmented MRI. TOF measured counts are represented by a collection n 1, , t are respectively the detector and the time bin indices. In absence of motion, g i,t follows a Poisson distribution of expectation ¯( is the TOF PET system response function at detector/time bin (i, t), τ is the scanning time, L i is the segment connecting the detectors of bin i, s i, t is the expected background events (random/scatter) at bin (i, t) and ⊂ Ω R 3 is a compact set representing the field of view.

Model with motion.
In practice, f and μ are affected by patient motion. For quasi-cyclic motion (respiratory, cardiac), acquired emission data are regrouped into n g gates, each of which corresponds to a patient state. Each TOF-PET data vector ( ) = ∈ = × N g g l it l i t n n n n , , , 1 , b t b t , = … l n 1, , g , corresponds to a single gate in which the patient is assumed to be static. The motion at gate l is represented by a diffeomorphism → ϕ R R : l 3 3 that deforms both the activity distribution f and the attenuation map μ. The deformed activities and μ-maps are f is the warping operator defined as . The number of detected counts g i, t, l at bin (i, t), gate l is a Poisson random variable of expectation where s i, t, l is the expected number of background events at bin (i, t), gate l, τ l is the duration of gate l and a i and H i t , are defined in (1).

Optimization problem.
The log-likelihood of the TOF-PET data is Maximum-likelihood joint image reconstruction/motion estimation (JRM) in TOF PET and non-TOP PET consists of solving the following optimization problem where L is either L TOF or L PET , = D D ng and D denotes the set of diffeorphism on R 3 . Solving (4) returns a single image f reconstructed from all gates g l , = … l n 1, , g . The reconstructed activity image at gate l is f f l lˆˆφ = ϕ W . It should be noted that (2) and therefore (3) depend on f and μ only via ϕ W f l and µ ϕ W l . The activity distribution f corresponds to an unobserved state consistent with μ (but not necessarily with any g l ). This led to a theoretical result for JRM in non-TOF PET in Bousse et al (2016, proposition 1), that can be naturally extended to TOF PET: if µ 1 and µ 2 are 2 different attenuation maps such that 2 1 µ µ ψ = for some ψ ∈ D, then there exists a bijection between the maximizers (when they exist) of ( returns similar reconstructed gates, and more specifically, JRM can be performed with a μ-map misaligned with each gate, provided this misaligned results from some diffeomorphism. Moreover, the estimated motion realigns the μ-maps to the data at each gate l. When μ is aligned with one gate g l0 , f can be initialized by maximizing the log-likelihood at gate l 0 with ϕ = Id l0 . If μ is misaligned with every gate, there is no consistent data to obtain a good initial estimate of f, thus rendering the optimization problem (4) more difficult. For example, in non-TOF PET (see (Bousse et al (2016, section IV.B)), it takes 50 to 100 iterations to solve (4) with a misaligned μ-map. The objective of this work is to see if TOF PET can facilitate the optimization.

Joint reconstruction scheme
2.2.1. Discretization. We used the same discretization scheme as in Bousse et al (2016), initially introduced in Jacobson and Fessler (2003) and Jacobson (2006).
( ) R M n m , denotes the space of real × n m matrices. The activity distribution and attenuation map are respectively represented by ∈ + R f n v and µ ∈ + R n v , where n v denotes the total number of voxels. A deformation ϕ is parametrized by a B-spline coefficient vector , , t is the probability that an annihilation occurring at voxel j is detected in bin (i, t). The expected number of counts is redefined as where i j , is the length of the intersection of L i with voxel j.
In comparison, the non-TOF PET discrete log-likelihood is To enforce image and motion smoothness, 2 penalty terms f U( ) and ( ) θ V are added to L TOF and L PET (quadratic penalties, see Bousse et al (2016)). The discrete JRM formulation is formulated as We solve (8) and (9) by alternating maximization in f and θ. We proceed as in Bousse et al (2016): a complete iteration consists of a maximization in θ-performed with a limitedmemory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) quasi-Newton (QN) line-search algorithm (Nocedal and Wright 2006, chapter 7)-and a maximization in f -performed with a motion compensated (MC) modified maximum-likelihood expectation-maximization (MMLEM) (De Pierro 1995) reconstruction from the gated data ( ) = g l l n 1 g using the current estimated motion parameter θ.
The gradient of L TOF in α l is similar to the PET case described in Bousse et al (2016) with a summation over the time index t: , . The motion parameter θ reduces to a single α. The aim is to reconstruct a single activity image with a misaligned μ-map. The estimated deformation corresponds to the misalignment between the input μ-map and g.
Each MC-MMLEM sub-optimimization was performed with 100 iterations which in our experiments was sufficient for convergence (no stopping rule was used), using the current α estimate. Prior to each MC-MMLEM sub-optimimization, the image f was reinitialized to a blank image to avoid attenuation artifacts due to incomplete motion estimation in the previous iteration. We used the Fortran implementation proposed in Zhu et al (1997) for each L-BFGS sub-optimimization, with 80 iterations (sufficient for convergence in our experiments).

Simulation set-up
Our simulation set-up is similar to Bousse et al (2016, section IV.B) with Poisson noise. We used the XCAT phantom to generate ground truth end-expiration activity f and attenuation map µ ( × × 160 160 42 volumes with 3.125 mm edge cubic voxels, corresponding to a 500 mm field of view, see figures 1(a) and (b)). We used a = × × n 53 53 16 c B-spline grid to parametrize the motion. A misaligned end-inhalation μ-map, μ (figure 1(c)), was also generated.
The spatial resolution of the PET was set to 6 mm FWHM and the temporal TOF resolution set to 500 ps FWHM. We used 10 TOF bins (332 ps width). The line integral operator L was adjusted to match the spatial resolution of the PET. We generated TOF and non-TOF data Poisson random vectors from f and µ as:

Figures 2 and 3 show the JRM reconstructed activities α W f
PET PET and W f TOF TOFα at several iterations (iteration 0 corresponds to the first MMLEM reconstruction with no motion compensation). It can be observed that PET JRM requires 50 to 100 iterations to reduce μ-map misalignment artifacts significantly whereas TOF JRM requires only 5. The relative differences μ µ − α W PET and μ µ − α W TOF are shown in figures 4 and 5 6 . Results show that the warped μ-map using TOF JRM estimated motion, i.e., μ α W TOF , becomes similar to µ in 5 iterations, whereas it takes more than 50 iterations for the warped μ-map using non-TOF JRM estimated motion. Note that since TOF PET forward and back projections are computationally