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

The activity distribution and the attenuation map respectively take the form of 2 functions $f\in {{\mathcal{C}}^{+}}$ and $\mu \in {{\mathcal{C}}^{+}}$, where ${{\mathcal{C}}^{+}}$ denotes the set of non-negative continuous functions defined on ${{\mathbb{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 $\boldsymbol{g}=\left({{g}_{i,t}}\right)_{i,t=1}^{{{n}_{\text{b}}},{{n}_{\text{t}}}}\in {{\mathbb{N}}^{{{n}_{\text{b}}}\times {{n}_{\text{t}}}}}$. The subscripts $i\in \left\{1,\ldots,{{n}_{\text{b}}}\right\}$ and $t\in \left\{1,\ldots,{{n}_{\text{t}}}\right\}$ are respectively the detector and the time bin indices. In absence of motion, g_i,t follows a Poisson distribution of expectation ${{\bar{g}}_{i,t}}(\,f,\mu )$:

$$\begin{eqnarray}&&{{\bar{g}}_{i,t}}(\,f,\mu )=\tau {{a}_{i}}(\mu ){{\mathcal{H}}_{i,t}}\,f+{{s}_{i,t}},\end{eqnarray}$$

with

$$\begin{eqnarray}&&{{\mathcal{H}}_{i,t}}\,f\triangleq {{\mathop{\int}^{}}_{\Omega }}f\left(\boldsymbol{r}\right){{h}_{i,t}}\left(\boldsymbol{r}\right)\,\text{d}\boldsymbol{r}\quad \text{and}\quad {{a}_{i}}(\mu )\triangleq \exp \left(-{{\mathop{\int}^{}}_{{{L}_{i}}}}\mu \left(\boldsymbol{r}\right)\,\text{d}\boldsymbol{r}\right)\end{eqnarray} \tag{ 1 }$$

where ${{h}_{i,t}}:{{\mathbb{R}}^{3}}\to {{\mathbb{R}}^{+}}$ 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 $\Omega \subset {{\mathbb{R}}^{3}}$ is a compact set representing the field of view.

In practice, f and μ are affected by patient motion. For quasi-cyclic motion (respiratory, cardiac), acquired emission data are regrouped into ${{n}_{\text{g}}}$ gates, each of which corresponds to a patient state. Each TOF-PET data vector $\boldsymbol{g}_{{l}}=\left({{g}_{i,t,l}}\right)_{i,t=1}^{{{n}_{\text{b}}},{{n}_{\text{t}}}}\in {{\mathbb{N}}^{{{n}_{\text{b}}}\times {{n}_{\text{t}}}}}$, $l=1,\ldots,{{n}_{\text{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 ${{\varphi}_{l}}:{{\mathbb{R}}^{3}}\to {{\mathbb{R}}^{3}}$ that deforms both the activity distribution f and the attenuation map μ. The deformed activities and μ-maps are ${{\mathcal{W}}_{{{\varphi}_{l}}}}\,f$ and ${{\mathcal{W}}_{{{\varphi}_{l}}}}\mu$, where ${{\mathcal{W}}_{\varphi}}:{{\mathcal{C}}^{+}}\to {{\mathcal{C}}^{+}}$ is the warping operator defined as ${{\mathcal{W}}_{\varphi}}h=h\circ \varphi$ for all $h\in {{\mathcal{C}}^{+}}$. The number of detected counts g_{i, t, l} at bin (i, t), gate l is a Poisson random variable of expectation

$$\begin{eqnarray}&&{{\bar{g}}_{i,t,l}}\left(\,f,{{\varphi}_{l}},\mu \right)\triangleq {{\tau}_{l}}{{a}_{i}}\left({{\mathcal{W}}_{{{\varphi}_{l}}}}\mu \right){{\mathcal{H}}_{i,t}}{{\mathcal{W}}_{{{\varphi}_{l}}}}\,f+{{s}_{i,t,l}},\end{eqnarray} \tag{ 2 }$$

where s_{i, t, l} is the expected number of background events at bin (i, t), gate l, ${{\tau}_{l}}$ is the duration of gate l and a_i and ${{\mathcal{H}}_{i,t}}$ are defined in (1).

The $\log$-likelihood of the TOF-PET data is

$$\begin{eqnarray}&&{{L}^{\text{TOF}}}\left(\,f,\boldsymbol{\varphi },\mu \right)=\underset{l=1}{\overset{{{n}_{\text{g}}}}{\mathop \sum}}\,\underset{i=1}{\overset{{{n}_{\text{b}}}}{\mathop \sum}}\,\underset{t=1}{\overset{{{n}_{\text{t}}}}{\mathop \sum}}\,{{g}_{i,t,l}}\log {{\bar{g}}_{i,t,l}}\left(\,f,{{\varphi}_{l}},\mu \right)-{{\bar{g}}_{i,t,l}}\left(\,f,{{\varphi}_{l}},\mu \right),\end{eqnarray} \tag{ 3 }$$

where $\boldsymbol{\varphi }\triangleq \left({{\varphi}_{l}}\right)_{l=1}^{{{n}_{\text{g}}}}$. Similarly, the non-TOF PET $\log$-likelihood is

$$\begin{eqnarray}&&{{L}^{\text{PET}}}\left(\,f,\boldsymbol{\varphi },\mu \right)=\underset{l=1}{\overset{{{n}_{\text{g}}}}{\mathop \sum}}\,\underset{i=1}{\overset{{{n}_{\text{b}}}}{\mathop \sum}}\,{{g}_{i,l}}\log {{\bar{g}}_{i,l}}\left(\,f,{{\varphi}_{l}},\mu \right)-{{\bar{g}}_{i,l}}\left(\,f,{{\varphi}_{l}},\mu \right),\end{eqnarray}$$

with ${{g}_{i,l}}=\sum_{t=1}^{{{n}_{\text{t}}}}{{g}_{i,t,l}}$ and ${{\bar{g}}_{i,l}}\left(\,f,{{\varphi}_{l}},\mu \right)=\sum_{t=1}^{{{n}_{\text{t}}}}{{\bar{g}}_{i,t,l}}\left(\,f,{{\varphi}_{l}},\mu \right)$.

Maximum-likelihood joint image reconstruction/motion estimation (JRM) in TOF PET and non-TOP PET consists of solving the following optimization problem

$$\begin{eqnarray}&&\left(\,\hat{f},\hat{\boldsymbol{\varphi }}\right)=\underset{f\in {{\mathcal{C}}^{+}},\boldsymbol{\varphi }\in \mathcal{D}}{\mathop{\text{arg}\,\text{max}}}\,L\left(\,f,\boldsymbol{\varphi },\mu \right)\end{eqnarray} \tag{ 4 }$$

where L is either ${{L}^{\text{TOF}}}$ or ${{L}^{\text{PET}}}$, $\mathcal{D}={{\mathcal{D}}^{{{n}_{\text{g}}}}}$ and $\mathcal{D}$ denotes the set of diffeorphism on ${{\mathbb{R}}^{3}}$. Solving (4) returns a single image f reconstructed from all gates $\boldsymbol{g}_{{l}}$, $l=1,\ldots,{{n}_{\text{g}}}$. The reconstructed activity image at gate l is ${{\mathcal{W}}_{{{{\hat{\varphi}}}_{l}}}}\,\hat{f}=\hat{f}\circ {{\hat{\varphi}}_{l}}$.

It should be noted that (2) and therefore (3) depend on f and μ only via ${{\mathcal{W}}_{{{\varphi}_{l}}}}\,f$ and ${{\mathcal{W}}_{{{\varphi}_{l}}}}\mu$. The activity distribution f corresponds to an unobserved state consistent with μ (but not necessarily with any $\boldsymbol{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 ${{\mu}_{1}}$ and ${{\mu}_{2}}$ are 2 different attenuation maps such that ${{\mu}_{2}}={{\mu}_{1}}\circ \psi$ for some $\psi \in \mathcal{D}$, then there exists a bijection between the maximizers (when they exist) of $\left(\,f,\boldsymbol{\varphi }\right)\mapsto L\left(\,f,\boldsymbol{\varphi },{{\mu}_{1}}\right)$ and the maximizers of $\left(\,f,\boldsymbol{\varphi }\right)\mapsto L\left(\,f,\boldsymbol{\varphi },{{\mu}_{2}}\right)$ such that each pair of maximizers $(\,{{\hat{f}}_{1}},{{\hat{\boldsymbol{\varphi }}}^{1}})$–$(\,{{\hat{{{f}_{2}}}}_{{}}},{{\hat{\boldsymbol{\varphi }}}^{2}})$, defined by this bijection, satisfies ${{\mathcal{W}}_{\hat{\varphi}_{l}^{1}}}\,{{\hat{f}}_{1}}={{\mathcal{W}}_{\hat{\varphi}_{l}^{2}}}\,{{\hat{f}}_{2}}$ and ${{\mathcal{W}}_{\hat{\varphi}_{l}^{1}}}{{\mu}_{1}}={{\mathcal{W}}_{\hat{\varphi}_{l}^{2}}}{{\mu}_{2}}$.

Proposition 1 in Bousse et al (2016) does not claim the existence and uniqueness of a maximizer for $\left(\,f,\boldsymbol{\varphi }\right)\mapsto L\left(\,f,\boldsymbol{\varphi },\mu \right)$. However, if $(\,{{\hat{f}}_{1}},{{\hat{\boldsymbol{\varphi }}}^{1}})$ is a likely candidate to maximize $\left(\,f,\boldsymbol{\varphi }\right)\mapsto L\left(\,f,\boldsymbol{\varphi },{{\mu}_{1}}\right)$, then there exists $(\,{{\hat{f}}_{2}},{{\hat{\boldsymbol{\varphi }}}^{2}})\in {{\mathcal{C}}^{+}}\times \mathcal{D}$ that is an equally likely candidate to maximize $\left(\,f,\boldsymbol{\varphi }\right)\mapsto L\left(\,f,\boldsymbol{\varphi },{{\mu}_{2}}\right)$ such that ${{\mathcal{W}}_{\hat{\varphi}_{l}^{1}}}\,{{\hat{f}}_{1}}={{\mathcal{W}}_{\hat{\varphi}_{l}^{2}}}\,{{\hat{f}}_{2}}$ and ${{\mathcal{W}}_{\hat{\varphi}_{l}^{1}}}{{\mu}_{1}}={{\mathcal{W}}_{\hat{\varphi}_{l}^{2}}}{{\mu}_{2}}$. Results from Bousse et al (2016) demonstrated that JRM with ${{\mu}_{1}}$ and ${{\mu}_{2}}={{\mu}_{1}}\circ \psi$ 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 $\boldsymbol{g}_{{{{l}_{0}}}}$, f can be initialized by maximizing the $\log$-likelihood at gate l₀ with ${{\varphi}_{{{l}_{0}}}}=\text{Id}$. 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.

We used the same discretization scheme as in Bousse et al (2016), initially introduced in Jacobson and Fessler (2003) and Jacobson (2006). ${{M}_{n,m}}\left(\mathbb{R}\right)$ denotes the space of real $n\times m$ matrices. The activity distribution and attenuation map are respectively represented by $\boldsymbol{f}\in \mathbb{R}_{+}^{{{n}_{\text{v}}}}$ and $\boldsymbol{\mu }\in \mathbb{R}_{+}^{{{n}_{\text{v}}}}$, where ${{n}_{\text{v}}}$ denotes the total number of voxels. A deformation φ is parametrized by a B-spline coefficient vector $\boldsymbol{\alpha }=\left(\boldsymbol{\alpha }^{{X}},\boldsymbol{\alpha }^{{Y}},\boldsymbol{\alpha }^{{Z}}\right)\in {{\mathbb{R}}^{3\times {{n}_{\text{c}}}}}$, where ${{n}_{\text{c}}}$ denotes the number of B-spline control points. The warping operator ${{\mathcal{W}}_{\varphi}}$ becomes a square matrix $\boldsymbol{W}_{{\boldsymbol{\alpha }}}\in {{M}_{{{n}_{\text{v}}},{{n}_{\text{v}}}}}\left(\mathbb{R}\right)$, defined by the voxel interpolating functions, the B-spline coefficients $\boldsymbol{\alpha }$ and the B-spline basis. The B-spline coefficient corresponding to the motion ${{\varphi}_{l}}$ at gate l is denoted $\boldsymbol{\alpha }_{{l}}$. The entire collection of B-spline coefficients is denoted $\boldsymbol{\theta }=\left(\boldsymbol{\alpha }_{{l}}\right)_{l=1}^{{{n}_{\text{g}}}}$.

The TOF operator ${{\mathcal{H}}_{i,t}}$ takes the form of a system matrix $\boldsymbol{H}\in {{M}_{{{n}_{\text{b}}}\times {{n}_{\text{t}}},{{n}_{\text{v}}}}}\left(\mathbb{R}\right)$ where ${{\left[\boldsymbol{H}\right]}_{(i-1){{n}_{\text{t}}}+t,j}}\triangleq {{p}_{i,j,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

$$\begin{eqnarray}&&{{\bar{g}}_{i,t,l}}\left(\boldsymbol{f},\boldsymbol{\alpha }_{{l}},\boldsymbol{\mu }\right)\triangleq \tau {{a}_{i}}\left(\boldsymbol{W}_{{\boldsymbol{\alpha }_{{l}}}}\boldsymbol{\mu }\right){{\left[\boldsymbol{H}\boldsymbol{W}_{{\boldsymbol{\alpha }_{{l}}}}\,\boldsymbol{f}\right]}_{(i-1){{n}_{\text{t}}}+t}}+{{s}_{i,t,l}},\end{eqnarray} \tag{ 5 }$$

where ${{a}_{i}}\left(\boldsymbol{\mu }\right)\triangleq \exp \left(-{{\left[\boldsymbol{L}\boldsymbol{\mu }\right]}_{i}}\right)$ and the line integral matrix $\boldsymbol{L}\in {{M}_{{{n}_{\text{b}}},{{n}_{\text{v}}}}}\left(\mathbb{R}\right)$ is defined by ${{\left[\boldsymbol{L}\right]}_{i,j}}={{\ell}_{i,j}}$ where ${{\ell}_{i,j}}$ is the length of the intersection of L_i with voxel j.

The discrete TOF $\log$-likelihood is derived from (3):

$$\begin{eqnarray}&&{{L}^{\text{TOF}}}\left(\,\boldsymbol{f},\boldsymbol{\theta },\boldsymbol{\mu }\right)=\underset{l=1}{\overset{{{n}_{\text{g}}}}{\mathop \sum}}\,\underset{i=1}{\overset{{{n}_{\text{b}}}}{\mathop \sum}}\,\underset{t=1}{\overset{{{n}_{\text{t}}}}{\mathop \sum}}\,{{g}_{i,t,l}}\log {{\bar{g}}_{i,t,l}}\left(\,\boldsymbol{f},\boldsymbol{\alpha }_{{l}},\boldsymbol{\mu }\right)-{{\bar{g}}_{i,t,l}}\left(\,\boldsymbol{f},\boldsymbol{\alpha }_{{l}},\boldsymbol{\mu }\right).\end{eqnarray} \tag{ 6 }$$

In comparison, the non-TOF PET discrete $\log$-likelihood is

$$\begin{eqnarray}&&{{L}^{\text{PET}}}\left(\,\boldsymbol{f},\boldsymbol{\theta },\boldsymbol{\mu }\right)\triangleq \underset{l=1}{\overset{{{n}_{\text{g}}}}{\mathop \sum}}\,\underset{i=1}{\overset{{{n}_{\text{b}}}}{\mathop \sum}}\,{{g}_{i,l}}\log {{\bar{g}}_{i,l}}\left(\,\boldsymbol{f},\boldsymbol{\alpha }_{{l}},\boldsymbol{\mu }\right)-{{\bar{g}}_{i,l}}\left(\,\boldsymbol{f},\boldsymbol{\alpha }_{{l}},\boldsymbol{\mu }\right)\end{eqnarray} \tag{ 7 }$$

with ${{g}_{i,l}}=\sum_{t=1}^{{{n}_{\text{t}}}}{{g}_{i,t,l}}$ and ${{\bar{g}}_{i,l}}\left(\,\boldsymbol{f},\boldsymbol{\alpha }_{{l}},\boldsymbol{\mu }\right)=\sum_{t=1}^{{{n}_{\text{t}}}}{{\bar{g}}_{i,t,l}}\left(\,\boldsymbol{f},\boldsymbol{\alpha }_{{l}},\boldsymbol{\mu }\right)$.

To enforce image and motion smoothness, 2 penalty terms $U\left(\,\boldsymbol{f}\right)$ and $V\left(\boldsymbol{\theta }\right)$ are added to ${{L}^{\text{TOF}}}$ and ${{L}^{\text{PET}}}$ (quadratic penalties, see Bousse et al (2016)). The discrete JRM formulation is formulated as

$$\begin{eqnarray}&&\,\left({{\hat{\boldsymbol{f}}}^{\text{TOF}}},{{\hat{\boldsymbol{\theta }}}^{\text{TOF}}}\right)=\underset{\boldsymbol{f}\geqslant \mathbf{0},\boldsymbol{\theta }}{\mathop{\text{arg}\,\text{max}}}\,{{\Phi}^{\text{TOF}}}\left(\,\boldsymbol{f},\boldsymbol{\theta },\boldsymbol{\mu }\right)\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\left(\,{{\hat{\boldsymbol{f}}}^{\text{PET}}},{{\hat{\boldsymbol{\theta }}}^{\text{PET}}}\right)=\underset{\boldsymbol{f}\geqslant \mathbf{0},\boldsymbol{\theta }}{\mathop{\text{arg}\,\text{max}}}\,{{\Phi}^{\text{PET}}}\left(\,\boldsymbol{f},\boldsymbol{\theta },\boldsymbol{\mu }\right)\end{eqnarray} \tag{ 9 }$$

with ${{\Phi}^{\text{TOF}}}\left(\,\,\boldsymbol{f},\boldsymbol{\theta },\boldsymbol{\mu }\right)={{L}^{\text{TOF}}}\left(\,\,\boldsymbol{f},\boldsymbol{\theta },\boldsymbol{\mu }\right)-\beta U\left(\,\,\boldsymbol{f}\right)-\gamma V\left(\boldsymbol{\theta }\right)$, ${{\Phi}^{\text{PET}}}\left(\,\,\boldsymbol{f},\boldsymbol{\theta },\boldsymbol{\mu }\right)={{L}^{\text{PET}}}\left(\,\,\boldsymbol{f},\boldsymbol{\theta },\boldsymbol{\mu }\right)-\beta U\left(\,\,\boldsymbol{f}\,\right)-$$\gamma V\left(\boldsymbol{\theta }\right)$. The JRM-reconstructed activity images corresponding to gate l are $\boldsymbol{W}_{{\hat{\boldsymbol{\alpha }}_{l}^{\text{TOF}}}}{{\hat{\boldsymbol{f}}}^{\text{TOF}}}$ for TOF PET and $\boldsymbol{W}_{{\hat{\boldsymbol{\alpha }}_{l}^{\text{PET}}}}{{\hat{\boldsymbol{f}}}^{\text{PET}}}$ for non-TOF PET.

We solve (8) and (9) by alternating maximization in $\boldsymbol{f}$ and $\boldsymbol{\theta }$. We proceed as in Bousse et al (2016): a complete iteration consists of a maximization in $\boldsymbol{\theta }$—performed with a limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) quasi-Newton (QN) line-search algorithm (Nocedal and Wright 2006, chapter 7)—and a maximization in $\boldsymbol{f}$—performed with a motion compensated (MC) modified maximum-likelihood expectation-maximization (MMLEM) (De Pierro 1995) reconstruction from the gated data $\left(\boldsymbol{g}_{{l}}\right)_{l=1}^{{{n}_{\text{g}}}}$ using the current estimated motion parameter $\boldsymbol{\theta }$.

The gradient of ${{L}^{\text{TOF}}}$ in $\boldsymbol{\alpha }_{{l}}$ is similar to the PET case described in Bousse et al (2016) with a summation over the time index t:

$$\begin{eqnarray}&&{{\nabla}_{\boldsymbol{\alpha }_{{l}}}}L\left(\,\boldsymbol{f},\boldsymbol{\theta },\boldsymbol{\mu }\right)=\underset{t=1}{\overset{{{n}_{\text{t}}}}{\mathop \sum}}\,\boldsymbol{J}_{{\boldsymbol{\alpha }_{{l}}}}{{\left(\,{{\overline{\boldsymbol{g}}}_{:,t,l}}\left(\,\boldsymbol{f},\boldsymbol{\alpha }_{{l}},\boldsymbol{\mu }\right)\right)}^{\top}}\left(\frac{\boldsymbol{g}_{{:,t,l}}}{{{\overline{\boldsymbol{g}}}_{:,t,l}}\left(\,\boldsymbol{f},\boldsymbol{\alpha }_{{l}},\boldsymbol{\mu }\right)}-\mathbf{1}\right)\end{eqnarray}$$

with $\boldsymbol{A}\left(\boldsymbol{\alpha }_{{l}},\boldsymbol{\mu }\right)=\text{diag}\left\{\left({{a}_{i}}\left(\boldsymbol{W}_{{\boldsymbol{\alpha }_{{l}}}}\boldsymbol{\mu }\right)\right)_{i=1}^{{{n}_{\text{b}}}}\right\}$, $\boldsymbol{g}_{{:,t,l}}=\left({{g}_{i,t,l}}\right)_{i=1}^{{{n}_{\text{b}}}}$, ${{\overline{\boldsymbol{g}}}_{:,t,l}}\left(\,\boldsymbol{f},\boldsymbol{\alpha }_{{l}},\boldsymbol{\mu }\right)=\left({{\bar{g}}_{i,t,l}}\left(\,\boldsymbol{f},\boldsymbol{\alpha }_{{l}},\boldsymbol{\mu }\right)\right)_{i=1}^{{{n}_{\text{b}}}}$, $\boldsymbol{H}_{{t}}\in {{M}_{{{n}_{\text{b}}},{{n}_{\text{v}}}}}\left(\mathbb{R}\right)$ defined by ${{\left[\boldsymbol{H}_{{t}}\right]}_{i,j}}={{p}_{i,j,t}}$, and

$$\begin{eqnarray}&&\boldsymbol{J}_{{\boldsymbol{\alpha }_{{l}}}}\left({{\overline{\boldsymbol{g}}}_{:,t,l}}\left(\,\,\boldsymbol{f},\boldsymbol{\alpha }_{{l}},\boldsymbol{\mu }\right)\right)=\tau \boldsymbol{A}\left(\boldsymbol{\alpha }_{{l}},\boldsymbol{\mu }\right)\left(-\text{diag}\left\{\boldsymbol{H}_{{t}}\boldsymbol{W}_{{\boldsymbol{\alpha }_{{l}}}}\,\boldsymbol{f}\right\}\boldsymbol{L}\boldsymbol{J}_{{\boldsymbol{\alpha }_{{l}}}}\left(\boldsymbol{W}_{{\boldsymbol{\alpha }_{{l}}}}\boldsymbol{\mu }\right)+\;\boldsymbol{H}_{{t}}\boldsymbol{J}_{{\boldsymbol{\alpha }_{{l}}}}\left(\boldsymbol{W}_{{\boldsymbol{\alpha }_{{l}}\,}}\,\boldsymbol{f}\right)\right).\end{eqnarray} \tag{ 10 }$$

The Jacobian matrices $\boldsymbol{J}_{{\boldsymbol{\alpha }}}\left(\boldsymbol{W}_{{\boldsymbol{\alpha }}}\,\boldsymbol{f}\right)$ and $\boldsymbol{J}_{{\boldsymbol{\alpha }}}\left(\boldsymbol{W}_{{\boldsymbol{\alpha }}}\boldsymbol{\mu }\right)$ are derived in Bousse et al (2016).

The overall scheme is briefly summarized in algorithm 1, for $\Phi={{\Phi}^{\text{TOF}}}$ or $\Phi={{\Phi}^{\text{PET}}}$. Each sub-maximization algorithm (MC-MMLEM and L-LBFGS) is run until convergence. A detailed version of the scheme can be found in Bousse et al (2016).

Algorithm 1:. Summary of JRM

Input: TOF PET gated data $\left(\boldsymbol{g}_{{l}}\right)_{l=1}^{{{n}_{\text{g}}}}$, attenuation map $\boldsymbol{\mu }$ (also image and motion smoothing priors weights β and γ).

Output: Activity image $\boldsymbol{f}$, B-spline motion parameter $\boldsymbol{\theta }$

$\boldsymbol{\theta }\leftarrow \mathbf{0}$ ;

$\boldsymbol{f}\leftarrow \text{MMLEM}\left(\boldsymbol{g}_{{1}}\right)$ (reconstruction from gate 1 w/o motion compensation) ;

for $n=1,\ldots,{{n}_{\text{max}}}$ do

$\boldsymbol{\theta }\leftarrow \text{arg}\,\text{ma}{{\text{x}}_{\boldsymbol{\theta }}}\Phi\left(\boldsymbol{f},\boldsymbol{\theta },\boldsymbol{\mu }\right)$ (with L-BFGS QN line-search) ;

$\boldsymbol{f}\leftarrow \text{MC}-\text{MMLEM}\left(\boldsymbol{g},\boldsymbol{\theta }\right)$ ;

end

When the task is to reconstruct from a single gate l₀ only, the projection model (5) can be simplified by ignoring the motion on the activity and warp the μ-map only, i.e. using a Poisson model with expectation

$$\begin{eqnarray}&&\check{\boldsymbol{g}}_{i,t,l_0}(\boldsymbol{f},\boldsymbol{\alpha}_{l_0},\boldsymbol{\mu}) \triangleq \tau a_i(\boldsymbol{W}_{\boldsymbol{\alpha}_{l_0}}\boldsymbol{\mu}) [\boldsymbol{H} \boldsymbol{f}]_{(i-1){n_{{\rm t}}} + t} + s_{i,t,_{l_0}} \end{eqnarray}$$

so that the Jacobian (10) expression is simplified. The reconstructed activity $\hat{\boldsymbol{f}}$ matches the data $\boldsymbol{g}_{{{{l}_{0}}}}$ without the need of being warped. This problem was addressed in Rezaei and Nuyts (2013). However, the generalization of this model to the multi gate case requires to estimate one activity image $\hat{\boldsymbol{f}}$ per gate, as opposed to JRM which only reconstruct one image $\boldsymbol{f}$, warped to each gate.