4D-PET reconstruction using a spline-residue model with spatial and temporal roughness penalties

A TAC can quite generally be modelled as the convolution of the AIF with a residue function,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f\left(t \right)=~K\int_{0}^{t}{{{C}_{I}}\left(s \right)R\left(t-s \right)}ds,\nonumber \end{align} \tag{ 1 }$$

where ${{C}_{I}}\left(t \right)$ is the AIF, $K$ is a proportionality constant interpreted as overall flow and $R(t)$ is the residue function, which describes the fraction of tracer remaining in the region at time t after entering it, and so provides information about the kinetics of radiotracer transport and metabolism processes (Gunn et al 2001).

The spline-residue model represents the residue function as a weighted sum of ${{N}_{S}}$ cubic B-spline functions, ${{\zeta }_{k}}\left(t \right)$, with a Dirac delta function added to account for the finite blood volume in the region:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle R\left(t \right)=~{{\mu }_{0}}\delta \left(t \right)+\sum\limits_{k=1}^{{{N}_{S}}}{{{\mu }_{k}}{{\zeta }_{k}}\left(t \right)},\nonumber \end{align} \tag{ 2 }$$

where ${{\mu }_{k}}$ is the coefficient of the kth B-spline function ${{\zeta }_{k}}\left(t \right)$, the coefficient ${{\mu }_{0}}$ equals ${{V}_{b}}/K$, and ${{V}_{b}}$ is the fractional blood volume. Thus from (1) the TAC is given by a weighted sum of spline-residue basis functions $\newcommand{\e}{{\rm e}} {{\eta }_{l}}\left(t \right)$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f\left(\boldsymbol{\theta },t \right)=\sum\limits_{l=0}^{{{N}_{S}}}{{{\theta }_{l}}{{\eta }_{l}}(t)},\nonumber \end{align} \tag{ 3 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\eta }_{l}}\left(t \right)=\left\{\begin{array}{@{}lllllllllllllllllllllllll@{}} {{C}_{I}}\left(t \right), & l=0 \nonumber \\ {{C}_{I}}\left(t \right)\otimes {{\zeta }_{l}}\left(t \right), & l>0 \nonumber \end{array} \right.,\nonumber \end{align} \tag{ 4 }$$

$\otimes$ denotes a convolution, and the coefficients ${{\theta }_{l}}$ are given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\theta }_{l}}=\left\{\begin{array}{@{}lllllllllllllllllllllllllllll@{}} {{V}_{b}}, & l=0 \nonumber \\ K{{\mu }_{l-1}}, & l>1 \nonumber \end{array} \right..\nonumber \end{align} \tag{ 5 }$$

Spline-residue, cubic B-spline and spectral model basis functions are compared in figure 1. The early spectral and spline-residue model basis function are very similar, while the later basis functions are rather different.

Zoom In Zoom Out Reset image size

Linear kinetic models represent the number of positron annihilation events in a given voxel j at time-frame m, ${{x}_{jm}}$, as a linear combination of pre-defined basis functions:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{x}_{jm}}\left({{\boldsymbol{\theta }}_{j}} \right)=\sum\limits_{k=1}^{{{N}_{B}}}{{{B}_{km}}{{\theta }_{jk}}},\nonumber \end{align} \tag{ 6 }$$

where ${{N}_{B}}$ is the total number of basis functions, ${{\theta }_{jk}}$ is the weighting factor of the kth basis function ${{B}_{k}}\left(t \right)$ in voxel j, and ${{B}_{km}}$ is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{B}_{km}}=\int_{{{t}_{ms}}}^{{{t}_{mf}}}{{{B}_{k}}\left(t \right)\exp \left(- \lambda t \right)}~dt,\nonumber \end{align} \tag{ 7 }$$

where ${{t}_{ms}}$ and ${{t}_{mf}}$ are respectively the start and end times of time-frame m, and λ is the radiotracer decay constant. The basis functions are pre-defined and therefore only the weighting factors ${{\boldsymbol{\theta }}_{j}}$ need to be calculated when fitting the models.

The expected number of photon counts on detector pair i in time-frame m, $\left\langle {{y}_{im}} \right\rangle$, can be estimated as a function of the model parameters using

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\langle {{y}_{im}}\left(\boldsymbol{\theta } \right) \right\rangle =\sum\limits_{j=1}^{{{N}_{V}}}{{{P}_{ij}}~{{x}_{jm}}\left({{\boldsymbol{\theta }}_{j}} \right)+~{{\varepsilon }_{im}}},\nonumber \end{align} \tag{ 8 }$$

where ${{N}_{V}}$ is the total number of voxels, ${{\varepsilon }_{im}}$ represents the erroneous counts measured by detector pair i in time-frame m (random coincidences and scattered photons), ${{P}_{ij}}$ are the elements of the ${{N}_{D}}\times {{N}_{V}}$ system matrix P, and ${{N}_{D}}$ is the number of detector pairs. The element ${{P}_{ij}}$ represents the probability of a pair of photons originating in voxel j being detected by detector pair i. The system matrix used here is independent of time, though time-dependent effects such as detector dead-time can be and sometimes are included in the system matrix calculation (Qi et al 1998).

Modelling the measured counts ${{y}_{im}}$ as independent Poisson-distributed variables, the log-likelihood function of the measured scanner data $L\left(\boldsymbol{y}|\boldsymbol{\theta } \right)$ (with a constant term omitted) is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle L\left(\boldsymbol y|\boldsymbol \theta \right)=\sum\limits_{m=1}^{{{N}_{T}}}{\sum\limits_{i=1}^{{{N}_{D}}}{({{y}_{im}}\ln \left(\left\langle {{y}_{im}}\left(\boldsymbol{\theta } \right) \right\rangle \right)-\left\langle {{y}_{im}}\left(\boldsymbol{\theta } \right) \right\rangle)}},~\nonumber \end{align} \tag{ 9 }$$

where ${{N}_{T}}$ is the number of time-frames. The most likely parameter values are obtained by iteratively maximizing $L\left(\boldsymbol{y}|\boldsymbol{\theta } \right)$ with respect to $\boldsymbol{\theta }$. Many basis functions can be included in linear kinetic models, potentially leading to overfitting of the data, and so we modify the objective function to include a temporal regularization term:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Phi \left(\boldsymbol{\theta } \right)=L\left(\boldsymbol{\theta }|\boldsymbol{y} \right)-\gamma \Gamma \left(\boldsymbol{\theta } \right),\nonumber \end{align} \tag{ 10 }$$

where $\gamma$ is a parameter controlling the trade-off between temporal smoothness and accurate TAC description, and $\Gamma \left(\boldsymbol{\theta } \right)$ is a temporal roughness penalty defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Gamma \left(\boldsymbol{\theta } \right)=\sum\limits_{j=1}^{{{N}_{V}}}{\Lambda \left({{\boldsymbol{\theta }}_{j}} \right)},\nonumber \end{align} \tag{ 11 }$$

with $\Lambda \left({{\boldsymbol{\theta }}_{j}} \right)$ being the penalty function for voxel j.

To efficiently maximize $\Phi \left(\boldsymbol{\theta } \right)$ we propose a methodology based on optimization transfer (Lange et al 2000), in which determination of $\boldsymbol{\theta }$ at iteration n is transferred to a surrogate function $Q\left(\boldsymbol{\theta }|{{\boldsymbol{\theta }}^{n}} \right)$ which minorizes the original log-likelihood function:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle Q\left(\boldsymbol{\theta }|{{\boldsymbol{\theta }}^{n}} \right)\leqslant \Phi (\boldsymbol{\theta }),\nonumber \end{align} \tag{ 12 }$$

with equality if and only if $\boldsymbol{\theta }={{\boldsymbol{\theta }}^{n}}$. By choosing ${{\boldsymbol{\theta }}^{n+1}}$ as the $\boldsymbol{\theta }$ value maximizing $Q\left(\boldsymbol{\theta }|{{\boldsymbol{\theta }}^{n}} \right)$, Lange et al (2000) showed that $\Phi \left(\boldsymbol{y}|{{\boldsymbol{\theta }}^{n+1}} \right)\geqslant \Phi \left(\boldsymbol{y}|{{\boldsymbol{\theta }}^{n}} \right)$.

We obtain the surrogate objective function by subtracting $\gamma \Gamma \left({{\boldsymbol{\theta }}_{j}} \right)$ from the surrogate function proposed by Wang and Qi (2010) for the minorization of $L\left(\boldsymbol{y}|{{\boldsymbol{\theta }}^{n}} \right)$ (with a constant term omitted). Doing so gives:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle Q\left(\boldsymbol{\theta }|{{\boldsymbol{\theta }}^{n}} \right)=\sum\limits_{j}{\left(\sum\limits_{i}{{{P}_{ij}}} \right)\left(\left[ \sum\limits_{m}{\left(\hat{x}_{jm}^{n+1}\ln \left({{x}_{jm}}\left({{\boldsymbol{\theta }}_{j}} \right) \right)-{{x}_{jm}}\left({{\boldsymbol{\theta }}_{j}} \right) \right)} \right]-{{\gamma }^{\prime }}\Lambda \left({{\boldsymbol{\theta }}_{j}} \right) \right),}\nonumber \end{align} \tag{ 13 }$$

where ${{\gamma }^{\prime }}=\frac{\gamma }{\left(\sum\nolimits_{i}{{{P}_{ij}}} \right)}$, n is the current iteration number and $\hat{x}_{jm}^{n+1}$ is the image obtained when updating the current image estimates using the MLEM algorithm (Shepp and Vardi 1982):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat{x}_{jm}^{n+1}=\frac{{{x}_{jm}}\left(\boldsymbol{\theta }_{j}^{n} \right)}{\sum\nolimits_{i}{{{P}_{ij}}}}\sum\nolimits_{i}{{{P}_{ij}}\frac{~{{y}_{im}}}{\left\langle {{y}_{im}}\left({{\boldsymbol{\theta }}^{n}} \right) \right\rangle }}.\nonumber \end{align} \tag{ 14 }$$

Because $Q\left(\boldsymbol{\theta }|{{\boldsymbol{\theta }}^{n}} \right)$ is separable in voxels the parameter values $\boldsymbol{\theta }_{j}^{n+1}$ can be obtained using

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \boldsymbol{\theta }_{j}^{n+1}=\underset{{{\boldsymbol{\theta }}_{j}}}{\mathop{\max }}\,\left[ \sum\limits_{m}{\hat{x}_{jm}^{n+1}\ln \left({{x}_{jm}}\left({{\boldsymbol{\theta }}_{j}} \right) \right)-{{x}_{jm}}\left({{\boldsymbol{\theta }}_{j}} \right)-{{\gamma }^{\prime }}\Lambda \left({{\boldsymbol{\theta }}_{j}} \right)} \right].\nonumber \end{align} \tag{ 15 }$$

Instead of maximizing (15), we adapt a method proposed by Matthews et al (2010) for the un-regularized case, which uses a weighted least squares approach to calculate $\boldsymbol{\theta }_{j}^{n+1}$. For the temporally regularized case, $\boldsymbol{\theta }_{j}^{n+1}$ is obtained by minimizing the penalized weighted least square error

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \boldsymbol{\theta }_{j}^{n+1}=\underset{{{\boldsymbol{\theta }}_{j}}}{\mathop{\min }}\,\left[ \underset{m}{\mathop \sum }\,{{w}_{jm}}{{\left(\hat{x}_{jm}^{n+1}-{{x}_{jm}}\left({{\boldsymbol{\theta }}_{j}} \right) \right)}^{2}}+{{\gamma }^{\prime }}\Lambda \left({{\boldsymbol{\theta }}_{j}} \right) \right],\nonumber \end{align} \tag{ 16 }$$

with weighting factors

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{w}_{jm}}=\frac{1}{{{x}_{jm}}\left(\boldsymbol{\theta }_{j}^{n} \right)}.\nonumber \end{align} \tag{ 17 }$$

Wang and Qi (2012a) have noted that monotonic convergence to the maximum-likelihood solution is not guaranteed when the model fitting step is modified to a weighted least squares problem, however convergence has been observed in practice (Matthews et al 2010).

Temporal roughness penalties have been used in previous 4D-PET reconstruction studies, but there is no consensus on what form $\Lambda \left({{\boldsymbol{\theta }}_{j}} \right)$ should take. Two penalty functions are explored in this work, the first being

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Lambda \left({{\boldsymbol{\theta }}_{j}} \right)={{\left| {{\boldsymbol{\theta }}_{j}} \right|}^{2}},\nonumber \end{align} \tag{ 18 }$$

which corresponds to L₂ regularization, and the second

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Lambda \left({{\boldsymbol{\theta }}_{j}} \right)=\int_{0}^{{{T}_{scan}}}{{{\left(\frac{{{\partial }^{2}}f\left({{\boldsymbol{\theta }}_{j}},t \right)}{\partial {{t}^{2}}} \right)}^{2}}}dt,\nonumber \end{align} \tag{ 19 }$$

where ${{T}_{scan}}$ is the dPET scan duration. The second penalty function is often used to fit splines to noisy data and has been applied to spline-based 4D-PET reconstruction (Nichols et al 2002, Li et al 2007), though not in an optimization transfer framework.

Both penalty functions can be expressed as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Lambda \left({{\boldsymbol{\theta }}_{j}} \right)=\boldsymbol{\theta }_{j}^{T}\boldsymbol{\Omega }{{\boldsymbol{\theta }}_{j}},\nonumber \end{align} \tag{ 20 }$$

where the superscript T indicates a matrix transpose. For penalty function (18) $\boldsymbol{\Omega }$ is the identity matrix I, while for (19) the elements of $\boldsymbol{\Omega }$, ${{\Omega }_{ab}}$, are

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\Omega }_{ab}}=\int_{0}^{{{T}_{scan}}}{{{{\ddot{B}}}_{a}}(t)~{{{\ddot{B}}}_{b}}(t)}dt,\nonumber \end{align} \tag{ 21 }$$

where ${{\ddot{B}}_{k}}(t)$ is the second-order time-derivative of basis function k. Using this notation the cost function in equation (16) can be re-expressed as a Tikhonov regularization problem

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \boldsymbol{\theta }_{j}^{n+1}=\underset{{{\boldsymbol{\theta }}_{j}}}{\mathop{\min }}\,\left[ \sum\limits_{m}{{{w}_{jm}}{{\left(\hat{x}_{jm}^{n+1}-{{x}_{jm}}\left({{\boldsymbol{\theta }}_{j}} \right) \right)}^{2}}+{{\gamma }^{\prime }}\boldsymbol{\theta }_{j}^{T}\boldsymbol{\Omega }{{\boldsymbol{\theta }}_{j}}} \right],\nonumber \end{align} \tag{ 22 }$$

and solved using the equation (Tikhonov et al 1995)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \boldsymbol{\theta }_{j}^{n+1}={{\left({{\boldsymbol{B}}^{T}}\boldsymbol{WB}+~{{\gamma }^{\prime }}\boldsymbol{\Omega } \right)}^{-1}}{{\boldsymbol{B}}^{T}}\boldsymbol{W}\hat{\boldsymbol x}_{j}^{n+1},\nonumber \end{align} \tag{ 23 }$$

where $\boldsymbol{B}$ has the ${{B}_{km}}$ elements defined in equation (7), $\boldsymbol{W}$ is the diagonal matrix $diag\left({{w}_{j1}},{{w}_{j2}},\ldots ,{{w}_{j{{N}_{T}}}} \right)$ and $\hat{\boldsymbol x}_{j}^{n+1}$ is an ${{N}_{T}}\times 1$ vector with mth element $\hat{x}_{jm}^{n+1}$.

To determine the best ${{\gamma }^{\prime }}$ to use in (22) for voxel j, $\gamma _{j}^{\prime }$, a range of ${{\gamma }^{\prime }}$ values can be defined in advance, and the one producing the fit with the lowest generalized cross validation (GCV) score (Wahba 1990) taken as optimal for this voxel:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \gamma _{j}^{\prime }=\underset{{{\gamma }^{\prime }}}{\mathop{\min }}\,\left[ GCV\left({{\gamma }^{\prime }} \right) \right]=\underset{{{\gamma }^{\prime }}}{\mathop{\min }}\,\left[ \frac{\sum\nolimits_{m}{{{w}_{mj}}{{\left(\hat{x}_{jm}^{n+1}-{{x}_{jm}}\left(\theta _{j}^{n+1} \right) \right)}^{2}}}}{({\rm Tr}{{\left(\boldsymbol{I}-\boldsymbol{B}({{\boldsymbol{B}}^{T}}\boldsymbol{WB}+~{{\gamma }^{\prime }}\boldsymbol{\Omega } \right)}^{-1}}{{\boldsymbol{B}}^{T}}\boldsymbol{W}){{)}^{2}}} \right],\nonumber \end{align} \tag{ 24 }$$

where ${\rm Tr}$(...) denotes the matrix trace, and the denominator of (24) corresponds to the effective degrees of freedom. The model fitting step is much faster than the image update step, and so selection of the ${{\gamma }^{\prime }}$ value for each voxel TAC in this semi-automatic manner does not greatly slow down the reconstruction.

Spatial regularization can be built into the reconstruction by replacing the image update step in equation (14) with a corresponding step from an iterative 2D- or 3D-maximum a posteriori (MAP) algorithm. These image updates are designed to maximize objective functions of the form $L\left(\boldsymbol{x}|\boldsymbol{y} \right)-\beta U(\boldsymbol{x})$ with respect to the image $\boldsymbol{x}$, where $\beta$ is a tunable parameter controlling the trade-off between resolution and noise, and $U\left(\boldsymbol{x} \right)$ is a concave function designed to penalize rough images,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle U\left(\boldsymbol{x} \right)=\frac{1}{4}\sum\limits_{j}{\sum\limits_{k\in {{\mathcal{N}}_{j}}}{{{z}_{jk}}\psi \left({{x}_{j}}-{{x}_{k}} \right)}},\nonumber \end{align} \tag{ 25 }$$

where ${{\mathcal{N}}_{j}}$ is the set of nearest neighbours of voxel j and ${{z}_{jk}}$ is a weighting factor equal to the normalized inverse distance between voxels j and k (Wang and Qi 2012b). Here we use the Lange function (Lange 1990)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \psi \left(\xi \right)=\left(\frac{\left| \xi \right|}{\delta }-\ln \left(1+\frac{\left| \xi \right|}{\delta } \right)~ \right)\delta ,\nonumber \end{align} \tag{ 26 }$$

which contains a further smoothing parameter $\delta$, and achieves good noise suppression in fairly uniform regions, while preserving edges better than the more widely used quadratic function, $\psi \left(\xi \right)={{\xi }^{2}}$ (Lange 1990).

The proposed temporally regularized 4D-PET reconstruction algorithm for linear kinetic models, subsequently referred to as nested-MAP reconstruction, can be summarized as follows.

• 1.  
  Start with an initial dPET image sequence estimate, in this work a sequence of uniform images with the radioactivity concentration in each voxel set to 100 Bq/cc.
• 2.  
  Update each image with one iteration of the MAP algorithm. Here MAP updates were performed via algorithm 1 of Wang and Qi (2012b) using their pixel-based rather than patch-based approach. Values of $\beta$ and $\delta$ are pre-selected.
• 3.  
  Fit a temporal model to each voxel TAC via (23), using either a fixed value of $\gamma$ or a range of $\gamma$ values and subsequently selecting the best value via GCV.
• 4.  
  Return to step 2 using the image voxel values predicted by the fitted model, ${{x}_{jm}}\left(\boldsymbol{\theta }_{j}^{n+1} \right)$, as the seed for next MAP update, and continue for either a fixed number of iterations or until the images have converged.

The 4D-XCAT2 digital phantom package (Segars et al 2010) was used to simulate a single slice of an NSCLC patient injected with the FMISO hypoxia tracer. Tracer activity concentrations in different regions and time-frames were chosen to match smoothed TACs taken from a clinical FMISO-dPET image of a patient with stage IV NSCLC. Specifically, lung and bone TACs were obtained from spherical ROIs of 3 cm diameter placed in healthy lung and spine regions in the real patient image, and an AIF TAC was taken from a cylindrical ROI of diameter 10 mm located in the centre of the descending aorta on five consecutive PET axial slices. Tumour TACs were obtained from irregularly-shaped ROIs considered to contain hypoxic and normoxic tumour tissue. The ROIs were drawn by an experienced radiologist and checked by a second radiologist.

For smoothing, each TAC except the AIF was fitted with cubic splines, adaptively placing the knots according to the algorithm proposed by Ralli et al (2017), and with irreversible two- and three-tissue compartment models having 3 and 5 rate-constants respectively. These compartment models are schematically drawn in figure 2 and subsequently referred to as 2C3K and 3C5K. The AIF TAC was fitted with the phenomenological three-exponential model of Feng et al (1993) alone. Weighted least squares was used for all fitting, with the weighting factors

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{w}_{m}}=~\frac{\Delta {{T}_{m}}{{e}^{-\lambda {{T}_{m}}}}}{{{a}_{m}}},\nonumber \end{align} \tag{ 27 }$$

where ${{a}_{m}}$ is the decay-corrected average activity concentration within the region during frame m, $\Delta {{T}_{m}}$ is the frame duration, and T_m is the mid-point of the mth frame (Chen et al 1991). Compartment model fitting was carried out using the Levenberg–Marquart algorithm, available in the MATLAB optimization toolbox (Mathworks).

Zoom In Zoom Out Reset image size

Model fit quality was assessed for each TAC using leave-one-out cross-validation, calculating the weighted residual sum of squares (RSS) error of the fit ith the weighting factors defined in (27). The Wald-Wolfowitz runs test was used to check whether any significant structure remained in the residuals of each model fit at the 5% significance level. Of the models that passed the runs test for a given patient TAC, the best model fit was taken to be the one with the lowest weighted RSS value.

Two patient phantoms were created with identical spatial geometries and voxel dimensions of 3.1  ×  3.1  ×  2.0 mm³. Lung, bone, blood, normoxic and hypoxic tumour regions were filled with noise-free ground-truth activity concentrations that varied with time according to model fits to the TACs obtained from the corresponding regions in the patient image. For the first 'realistic' phantom the best model fits to the different TACs were used, while for the second 'simplified' phantom, 2C3K model fits were used instead. The fitted curves were binned into a (1  ×  30 s, 6  ×  5 s, 6  ×  20 s, 7  ×  60 s, 10  ×  120 s, 3  ×  300 s) time-frame sequence followed by two additional 600 s frames at 2 and 4 h post-injection. This frame sequence matches the dynamic imaging protocol of the clinical dPET scan from which the phantom TACs were derived: following this protocol, the patients were injected with the FMISO tracer 30 s into scanning.

To illustrate the phantom geometry, an image of the final time-frame of the realistic phantom is shown in figure 3(a). We have used the phantoms to study the performance of linear model-based versus 2C3K-based 4D-PET reconstruction when all the underlying TACs take realistic shapes (realistic phantom), and when they are all described by fits of the 2C3K model (simplified phantom).

Zoom In Zoom Out Reset image size

dPET sinograms representative of those produced by an mMR PET-MR scanner (Siemens Healthcare, Erlangen, Germany) were generated for both phantoms using the PET-SORTEO Monte-Carlo simulation package (Reilhac et al 2004), which has been validated for the Siemens mMR scanner (Reilhac et al 2016). Fifty noise realizations of dynamic-PET sinogram data were generated for each phantom, including effects of scattered photons, random co-incidences and attenuation. No patient motion was simulated, the focus of the current study being to evaluate the effectiveness of noise suppression using 4D-PET reconstruction. On average, the total number of counts in each noise realization was approximately 3.5 million for both single-slice phantoms studied.

Attenuation and normalization correction sinograms were obtained respectively from an attenuation map of the patient phantom, and from simulated detector counts generated using PET-SORTEO for a 20 min scan of a cylindrical phantom containing a uniform activity concentration. The attenuation and normalization corrections were then modelled as part of the system matrix. Numbers of scattered photons and random coincidences were estimated using the single-scatter simulation algorithm (Watson 2000) and a delayed co-incidence window respectively.

From each simulated realization of FMISO patient phantom PET scanner data, dPET image sequences were reconstructed in 3.1  ×  3.1  ×  2.0 mm³ voxels using both conventional MAP and 4D nested-MAP algorithms, running each algorithm for 30 iterations. Nested-MAP 4D reconstructions were performed using the non-linear 2C3K model, and the linear adaptive-knot cubic B-spline, spectral and spline-residue models. For each linear model, reconstructions were carried out using the temporal regularisation penalties of equations (18) and (19), and with no temporal regularisation. For 2C3K model-based reconstructions no temporal regularization was used because this model is much more constrained than the linear ones. For nested-MAP reconstruction based on the 2C3K model, the non-linear model fitting step was performed using the Levenberg-Marquardt algorithm instead of (23). Spectral, spline-residue and 2C3K model-based 4D reconstructions require image-derived or blood-sampled AIFs, which were obtained here by fitting the three-exponential model of Feng et al (1993) to TACs obtained from conventionally (not 4D) reconstructed MAP images, for ROIs placed in the left ventricle.

Spectral model-based reconstructions were carried out using 100 basis functions with exponential decay constants spaced logarithmically between 1.1  ×  10⁻⁴ s⁻¹ (the decay constant of ¹⁸F) and 0.01 s⁻¹. For spline-based reconstructions, voxel-specific knot locations were selected using the adaptive-knot placement algorithm proposed by Ralli et al (2017), which for cubic splines places knots along equal segments of the integral of the 4th root of the 4th derivative of a TAC according to theorem XII.5 of De Boor (1978). For each voxel 11 free knots were positioned by applying the algorithm to the TAC obtained for that voxel from conventional MAP-reconstructed images.

Basis functions for the spline-residue model were obtained by placing 4 knots at the beginning and end of each TAC, as well as the point where the TAC starts to rise, to handle discontinuities, and positioning an additional 6 uniformly spaced knots between the initial rise and end points of the TAC. The B-splines associated with these knots were convolved with the AIF to calculate the spline-residue basis functions. Preliminary work fitting the spline-residue model to synthetic noisy FMISO TACs led to the choice of 6 additional knots, this number generating fits that best matched the ground-truth. Uniformly spaced knots performed well, perhaps because the residue function for a given TAC varies considerably less than the TAC itself.

The spatial regularization parameters $\beta$ and $\delta$, defined in (26), were both set to 0.1, a choice that produced the best contrast-to-noise ratio in images of digital phantom similar to the NEMA image quality phantom (National Electrical Manufacturers Association 2013) reconstructed using the MAP algorithm from simulated PET detector data generated with PET-SORTEO.

For each linear model and temporal roughness penalty used in the 4D reconstruction algorithms, the range of ${{\gamma }^{\prime }}$ values {${{\gamma }^{\prime }}=$ 0.001, 0.002, ..., 0.01} all produced good fits to TACs from conventional MAP-reconstructed patient phantom images. At each iteration of the temporally regularized nested-MAP reconstructions, therefore, the ${{\gamma }^{\prime }}$ value used for each voxel was individually selected from those ten as the one that minimized the GCV score of the model fit to that voxel's TAC, as described in section 2.2.

2.5.1. Image quality metrics

To characterize the accuracy of the reconstructed images, the average absolute bias of imaged activity concentrations over the scan time-course was calculated for every voxel j:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\left[ {\rm Image}\,{\rm Bias} \right]}_{j}}=\frac{1}{{{T}_{scan}}}\sum\limits_{m}{\Delta {{T}_{m}}|{{{\bar{a}}}_{jm}}-a_{jm}^{true}|},\nonumber \end{align} \tag{ 28 }$$

where ${{\bar{a}}_{jm}}$ is the mean activity concentration in voxel j at time-frame m in all 50 repeat image sequences, and $a_{jm}^{true}$ is the true activity concentration.

The noise in each voxel j at every time-frame m was calculated using the weighted standard deviation ${{\sigma }_{w,jm}}$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\sigma }_{w,jm}}=~{{\left(\frac{\sigma _{jm,~measured}^{2}\Delta {{T}_{m}}}{a_{jm}^{true}} \right)}^{\frac{1}{2}}},\nonumber \end{align} \tag{ 29 }$$

where $\sigma _{jm,measured}^{2}$ is the variance amongst the 50 repeat ${{a}_{jm}}$ values, and the weighting factors $a_{jm}^{true}/\Delta {{T}_{m}}$ are based on the dynamic-PET noise model of Chen et al (1991) and nominally account for intrinsic variations in noise between frames. Then the average noise for a given voxel j was calculated across all ${{N}_{T}}$ time-frames

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\left[ {\rm Image}\,{\rm Noise} \right]}_{j}}={{\sigma }_{wj}}=\frac{1}{{{N}_{T}}}\sum\limits_{m=1}^{{{N}_{T}}}{{{\sigma }_{w,jm}}}.\nonumber \end{align} \tag{ 30 }$$

Overall bias and noise were characterised as mean absolute bias and mean image noise ($\left\langle {{\sigma }_{w}} \right\rangle$) averaged over all voxels within the patient (the whole patient region). The same measures were also averaged over the tumour region alone, usually the primary focus of oncological dPET studies. Normalised mean absolute bias and $\left\langle {{\sigma }_{w}} \right\rangle$ values were expressed as percentages of the mean ground-truth activity averaged over all time-frames and all voxels in the whole patient or tumour regions. To facilitate algorithm inter-comparison, normalised mean image noise and absolute bias values were computed at every iteration of all 4D reconstructions. Corresponding values were also calculated for just the first 120 s of the scans, to assess algorithm performance at early time-points.

To assess the convergence of the reconstruction process, the mean square error (MSE) was calculated for each image voxel j at time-frame m and iteration n of the nested-MAP reconstructions

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm MSE}_{jm}^{n}=\frac{1}{{{N}_{r}}}\sum\limits_{k=1}^{{{N}_{r}}}{{{\left(a_{jm}^{true}-a_{jm,k}^{n} \right)}^{2}}},\nonumber \end{align} \tag{ 31 }$$

where ${{N}_{r}}$ is the number of noise realizations. At each iteration, a weighted sum of ${\rm MSE}_{jm}^{n}$ over all time-frames was calculated for every individual voxel, and these values were summed over all image voxels to create a single total MSE measure, TMSE:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm TMS}{{{\rm E}}^{n}}=\sum\limits_{j=1}^{{{N}_{V}}}{\sum\limits_{m=1}^{{{N}_{T}}}{\frac{\Delta {{T}_{m}}{\rm MSE}_{jm}^{n}}{a_{jm}^{true}}}},\nonumber \end{align} \tag{ 32 }$$

where ${{N}_{V}}$ is the number of image voxels. To check that the images did not change substantially after the first 30 reconstruction iterations explored throughout most of this study, each nested-MAP reconstruction was run for a further 10 iterations. Fractional changes in TMSE from one iteration n to the next n  +  1 were calculated as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle [{\rm fractional}\,{\rm TMSE}\,{\rm change}\,{\rm at}\,{\rm iteration}\,n]=\frac{{\rm TMS}{{{\rm E}}^{n}}-{\rm TMS}{{{\rm E}}^{n+1}}}{{\rm TMS}{{{\rm E}}^{n}}},\nonumber \end{align} \tag{ 33 }$$

and plotted as a function of iteration number for n  =  1 to n  =  39, positive TMSE changes corresponding to reductions in the total error, and negative changes to increases. Additionally, the fractional change between iterations 30 and 40 was calculated as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle [{\rm fractional}\,{\rm TMSE}\,{\rm change}\,{\rm at}\,{\rm iteration}\,{\rm 30/40}]=\frac{{\rm TMS}{{{\rm E}}^{30}}-{\rm TMS}{{{\rm E}}^{40}}}{{\rm TMS}{{{\rm E}}^{30}}}.\nonumber \end{align} \tag{ 34 }$$

2.5.2. Parametric map quality metrics

FMISO uptake kinetics are often determined using the 2C3K model (Wang et al 2009, McGowan et al 2017). Following this approach, we have fitted the 2C3K model to voxel TACs obtained from all the image sequences reconstructed using the different algorithms. The voxels studied were those lying within the phantom sub-volume shown in figure 3(b), which contains the tumour-like region and surrounding lung and is therefore of the greatest interest. Then we determined the accuracy and precision of uptake kinetics as characterised by the 2C3K model fits to the reconstructed images versus fits of the same model to the ground-truth phantom TACs, studying the 2C3K model rate-constants shown in figure 2(a) together with the flux constant

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{k}_{flux}}=\frac{{{K}_{1}}{{k}_{3}}}{{{k}_{2}}+{{k}_{3}}}.\nonumber \end{align} \tag{ 35 }$$

The bias and noise in fitted values of each parameter q  =  K₁, k₂, k₃, k_flux were calculated for each voxel j lying within the image sub-volume:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\left[ {\rm Parameter}\,{\rm Bias} \right]}_{j}}={{\bar{q}}_{j}}-q_{j}^{true},\nonumber \end{align} \tag{ 36 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\left[ {\rm Parameter}\,{\rm Noise} \right]}_{j}}=~{{\sigma }_{{{q}_{j}}}},\nonumber \end{align} \tag{ 37 }$$

where $q_{j}^{true}$ is the ground-truth value of kinetic parameter q in voxel j, ${{\bar{q}}_{j}}$ is the mean of the q_j values obtained for voxel j from each of the 50 reconstructed image sequences, and ${{\sigma }_{{{q}_{j}}}}$ is the standard deviation of these 50 q_j values. Then the [Parameter Bias]_j and [Parameter noise]_j measures were averaged over all the voxels within the analysed phantom sub-volume to characterize the overall performance of each reconstruction algorithm. These metrics were also averaged over a region containing hypoxic tumour alone, to assess the performance of the algorithms specifically within the tumour region.

Figure 4 shows fits of the 2C3K, 3C5K and cubic spline models to the real TAC data obtained from normoxic and hypoxic tumour and healthy lung and spine regions of the imaged NSCLC patient. Runs test results and leave-one-out cross-validation weighted RSS scores are listed in table 1.

Zoom In Zoom Out Reset image size

Fits of the adaptive-knot spline and 3C5K models had the lowest cross-validation scores for the healthy tissue and tumour regions respectively, and were therefore used to represent the ground-truth TACs for these regions in the realistic phantom. All the model fits used in the realistic phantom passed the runs test. The worst leave-one-out cross-validation scores were obtained for the 2C3K model fits, which only passed the runs test for the hypoxic tumour region.

Iteration-by-iteration plots of image bias versus noise, averaged over all time-frames and patient voxels, are shown in figure 5 for realistic and simplified phantom image sequences reconstructed using temporally regularized nested-MAP 4D algorithms based on the cubic spline, spectral and spline-residue linear models. Each plot compares results obtained for one phantom and one reconstruction algorithm using either no temporal roughness penalty or the penalty functions of equations (20) or (21).

Zoom In Zoom Out Reset image size

For spectral and spline-residue-based 4D reconstructions, the ${{\left| \boldsymbol{\theta } \right|}^{2}}$ temporal roughness penalty of equation (18) produced substantially less noisy images than the other penalty options, at similar levels of bias, and was therefore considered the best penalty function for these algorithms. For cubic spline-based 4D reconstructions, however, the integrated square derivative penalty function of equation (19) was viewed as the best penalty function, since it produced the least biased images at noise-levels only slightly higher than obtained using the ${{\left| \boldsymbol{\theta } \right|}^{2}}$ penalty.

In figure 6, noise is plotted against bias for images of the realistic phantom reconstructed using each linear model-based 4D algorithm and its associated optimal temporal roughness penalty, and using with the 2C3K-based 4D algorithm. Separate plots are shown for noise and bias measures averaged over the whole patient or the tumour regions, and averaged over the whole scan time or just the first 120 s. Image quality curves for reconstructions based on the spline-residue model were substantially better than those obtained for reconstructions based on cubic splines or the 2C3K model. And at early time-points and within the tumour region, the image quality curves for the spline-residue-based algorithm were also better than those for reconstructions based on the spectral model. However, the spectral model had a slight edge when the image quality metrics were averaged across the whole phantom and scan duration.

Zoom In Zoom Out Reset image size

Corresponding data are shown in figure 7 for reconstructions of the simplified phantom. For this phantom the spline-residue model produced slightly better results than the spectral model when considered across the whole phantom and scan duration; but 4D reconstructions based on the 2C3K model achieved the lowest bias values, unsurprisingly since the phantom kinetics are 2C3K-based.

Zoom In Zoom Out Reset image size

Figure 8 shows fractional changes in TMSE as a function of iteration number for the 2C3K- and linear model-based reconstructions of the realistic phantom, using the optimal temporal roughness penalties for the linear models. After 30 iterations the fractional change in TMSE per iteration was very small for reconstructions based on the spectral, spline-residue and 2C3K models. Total fractional changes in TMSE between iterations 30 and 40, calculated with (34), were 0.013 for the spectral model, 0.017 for the spline-residue model, 0.007 for the 2C3K model and  −0.030 for cubic splines. Thus continuing reconstruction beyond 30 iterations led to small improvements at best, and in the case of the spline-based reconstruction a 3% worsening in TMSE, most likely due to noise amplification at the later iterations.

Zoom In Zoom Out Reset image size

Figure 9 shows voxel-by-voxel spatial plots of absolute bias and noise in k_flux parametric maps of the realistic phantom sub-region shown in figure 3(b), obtained from image sequences reconstructed using the conventional (non-4D) MAP algorithm and nested-MAP 4D algorithms based on the 2C3K model and linear models used with their optimal temporal roughness penalties. The voxel-by-voxel noise and bias plots calculated for the simplified phantom were similar. The results in figure 9 are summarized in table 2, which lists the signed bias and noise in the k_flux parametric maps averaged over the entire phantom sub-region. In both figure 9 and table 2 the bias and noise are expressed as percentages of the ground-truth k_flux value, averaged over the entire sub-region. Of all the algorithms compared, spline-residue model-based reconstruction achieved the lowest average bias and second lowest noise levels in the k_flux parametric maps, the spectral model-based algorithm achieving lower noise but greater bias.

Table 2. Mean signed parameter bias and standard deviation (S.D.) in k_flux parametric maps of the realistic phantom, averaged over the sub-region shown in figure 3(b). The maps were obtained from image sequences reconstructed using the conventional (non-4D) MAP algorithm, and the 2C3K and linear model-based nested-MAP reconstructions (with their optimal temporal roughness penalties). Bias and S.D are expressed as percentages of the average ground-truth k_flux value across the phantom sub-region. The lowest bias and S.D values are highlighted in bold and underlined.

Algorithm	Kinetic model fitted during reconstruction	Bias (%)	S.D. (%)
MAP	None	7.9	15.8
Nested-MAP	2C3K	−4.2	8.4
Nested-MAP	Adaptive-knot cubic B-splines ($\int{{{\left[ [{{f}^{\left(2 \right)}}\left(t \right) \right]}^{2}}}dt$ penalty)	2.6	12.1
Nested-MAP	Spectral model (${{\left| \theta \right|}^{2}}$ penalty)	−  3.2	5.7
Nested-MAP	Spline-residue model (${{\left| \theta \right|}^{2}}$ penalty)	−0.1	7.7

Zoom In Zoom Out Reset image size

The average normalised bias and noise (standard deviation) in 2C3K model parameters, for fits to the TACs of every voxel of the realistic phantom sub-region of figure 3(b) in images reconstructed using the different algorithms, is shown in figure 10 for the individual 2C3K rate constants and the composite k_flux parameter. For each parameter, these values are normalised as fractions of the values achieved by conventional (non-4D) MAP reconstruction for the same parameter, thus showing the extent to which each nested-MAP reconstruction improves on the conventional MAP reconstruction for each kinetic parameter. Equivalent plots for the simplified phantom are shown in figure 11. Bias and noise values averaged over the smaller hypoxic tumour sub-volume alone are shown in figure 12 for the realistic and simplified phantoms.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the realistic phantom it can be seen from figure 10 that for five of the eight combinations of bias/noise and kinetic parameters analysed, bias or noise averaged across the whole phantom was reduced more than 50% by using 4D reconstruction based on the spline-residue model rather than conventional (non-4D) MAP-reconstruction. Furthermore, for 5/8 combinations spline-residue 4D reconstruction produced better results than any of the other reconstruction algorithms studied, while spectral-based 4D reconstruction produced the best results for 2/8. Compared to conventional reconstruction, spline-residue-based 4D reconstruction did not notably increase bias or noise for any of the eight combinations, whereas the bias and noise in fitted K₁ parameter values rose substantially above the conventional MAP levels when 4D reconstruction was carried out using the other temporal models. 4D reconstruction based on the 2C3K model generated the most biased kinetic parameters, and higher levels of noise than spline-residue-based reconstruction. From figure 12 it can be seen that for the hypoxic tumour region alone, the spline-residue- and spectral-based 4D reconstructions each achieved the best results for 3/8 of the bias/noise and kinetic parameter combinations analysed, and cubic spline-based 4D reconstruction the best results for 2/8.

For the simplified phantom, whose kinetics entirely follow the 2C3K model, 4D reconstruction based on the 2C3K model unsurprisingly achieved the lowest levels of bias for fitted 2C3K model kinetics parameters (see figures 11 and 12). This algorithm also achieved useful reductions in average noise levels for 2C3K parameter values throughout the simplified phantom compared to conventional MAP reconstruction, although not within the hypoxic tumour region. 4D reconstruction based on the spline-residue and spectral models achieved the greatest reductions in noise, but outside of the hypoxic tumour region the bias in fitted K₁ and k₂ values was larger for these algorithms than for conventional MAP reconstruction, particularly so for the spline-residue model.