Joint PET-MR respiratory motion models for clinical PET motion correction

Patient motion due to respiration can lead to artefacts and blurring in positron emission tomography (PET) images, in addition to quantification errors. The integration of PET with magnetic resonance (MR) imaging in PET-MR scanners provides complementary clinical information, and allows the use of high spatial resolution and high contrast MR images to monitor and correct motion-corrupted PET data. In this paper we build on previous work to form a methodology for respiratory motion correction of PET data, and show it can improve PET image quality whilst having minimal impact on clinical PET-MR protocols. We introduce a joint PET-MR motion model, using only 1 min per PET bed position of simultaneously acquired PET and MR data to provide a respiratory motion correspondence model that captures inter-cycle and intra-cycle breathing variations. In the model setup, 2D multi-slice MR provides the dynamic imaging component, and PET data, via low spatial resolution framing and principal component analysis, provides the model surrogate. We evaluate different motion models (1D and 2D linear, and 1D and 2D polynomial) by computing model-fit and model-prediction errors on dynamic MR images on a data set of 45 patients. Finally we apply the motion model methodology to 5 clinical PET-MR oncology patient datasets. Qualitative PET reconstruction improvements and artefact reduction are assessed with visual analysis, and quantitative improvements are calculated using standardised uptake value (SUVpeak and SUVmax) changes in avid lesions. We demonstrate the capability of a joint PET-MR motion model to predict respiratory motion by showing significantly improved image quality of PET data acquired before the motion model data. The method can be used to incorporate motion into the reconstruction of any length of PET acquisition, with only 1 min of extra scan time, and with no external hardware required.

MR images on a data set of 45 patients. Finally we apply the motion model methodology to 5 clinical PET-MR oncology patient datasets. Qualitative PET reconstruction improvements and artefact reduction are assessed with visual analysis, and quantitative improvements are calculated using standardised uptake value (SUV peak and SUV max ) changes in avid lesions.
We demonstrate the capability of a joint PET-MR motion model to predict respiratory motion by showing significantly improved image quality of PET data acquired before the motion model data. The method can be used to incorporate motion into the reconstruction of any length of PET acquisition, with only 1 min of extra scan time, and with no external hardware required.
Keywords: PET-MR, motion correction, lesion detection, motion models, respiratory motion (Some figures may appear in colour only in the online journal)

Introduction
PET acquisitions can be adversely affected by respiratory motion. Avid lesions in the thorax and abdomen can be blurred in images, and artefacts and tracer uptake quantification errors may be introduced. Motion correction by gating (Klein et al 1996, Boucher et al 2004, Bai and Brady 2011) requires registration of images at different motion states, but this requires long acquisition times to ensure good contrast-to-noise.
The emergence of PET-MR scanners allows the use of high contrast MR images to track respiratory motion and subsequent MR-based registration to provide deformation fields to correct simultaneously acquired PET data (Würslin et al 2013, Grimm et al 2015. To aid motion detection, MR tagging has been proposed, a technique that creates tags in MR images which can be tracked through respiration to provide the deformation information (Guerin et al 2011, Chun et al 2012. In these cases, motion tracking is done with fast dynamic MR sequences, but this means other clinical MR sequences cannot be acquired during this time. Motion models allow estimation of dense deformation fields from a simple surrogate signal that can be acquired throughout a scan (McClelland et al 2013). A motion model is built by forming a correspondence between physical deformations and surrogate data. Motion information can then be estimated at times when motion-capturing images are not being acquired, if the same surrogate is acquired continuously. This type of 'motion model' methodology allows the motion information to be gathered in a short space of time in a 'model-building' scan, allowing other clinical MR to be acquired whilst PET data are being collected. One source of surrogate data is external respiratory monitors, such as RPM (Real-time Position Management, Varian Medical Systems Inc.), pressure belt, or spirometer. Some work uses MR-derived measurements as the surrogate, by building a correspondence model between 3D MR-based deformations and either interleaved 2D navigator images (King et al 2012), or a measure derived directly from the dynamic MR images (Balfour et al 2015). In these cases, other MR sequences need to be altered to allow for a continuous acquisition of the surrogate data.
In this paper we investigate the ability of a respiratory motion model built from only 1 min of simultaneously acquired PET and MR data to capture intra-and inter-cycle variability, using the continuously acquired PET data itself as the surrogate, to predict motion during a PET scan. We then use this information to form a motion-compensated PET image. Many current clinical PET-MR protocols now use a fast 'whole body' approach, with acquisition times of only 3 min per bed position, so having minimal impact on the clinical protocol is important. We build on previous work (Manber et al 2015) in which we used a 4 min PET-MR scan to estimate respiratory motion by binning dynamic MR images that corresponded to the PET gates.
Our new methodology improves on many aspects of our previously published work. Firstly, we now use a continuous correspondence model, rather than using discrete bins. With our old method, only one MR image was chosen per slice location per respiratory bin as a representation of the moving anatomies at that point in the respiratory cycle. The continuous motion model now allows interpolation at any value of the model surrogate, and cases where a bin may contain no MR data do not occur.
The new continuous motion model allows extrapolation, estimating deformation fields at values of the surrogate signal that were not used as input to the model-building sequence, which has two advantages. Firstly, all of the PET data can be used in the image reconstruction, and secondly, the attenuation μ-map can be used in the PET reconstruction even when the surrogate during acquisition is outside of the normal breathing range.
Our new work is also more robust to registration errors. With the previous method and binning based schemes, deformation fields that mapped any motion state back to the reference state are provided by a single non-rigid registration. Now all MR images are registered back to a reference image and all deformation fields are used in the model. This means if any single registration fails, it has little effect on the overall model.
The new approach is more practical in a clinical setting. The PET and MR systems run on two separate clocks, so to temporally align the data, we previously collected extra MR sequences for clock syncing purposes. The time shift/drift between the PET and MR clocks is now accounted for as part of the optimisation scheme to build the motion model so only 1 min of extra data is needed to capture respiratory motion.
In this paper we present a novel methodology for practical PET respiratory motion correction, by acquiring only 1 min of extra PET-MR data to build a patient specific PET-MR motion model. We test four different types of models, using either 1 or 2 surrogate signals, with a linear or polynomial model, to find which provides the best estimation of motion. To evaluate our method we apply our PET-MR motion models to both MR acquired during and outside of the model-building acquisition to test model-fit and model-prediction accuracy on 45 patients. The models are then tested on PET acquisitions of 5 patients in a clinical setting by extending one 3 min bed position scan to 4 min to acquire the extra data to form the motion model, then using this in a motion-compensated reconstruction of the whole 4 min scan. We also present a patient-specific binning scheme to efficiently bin PET data prior to reconstruction based on the breathing during acquisition.

Motion models
A motion model can be defined as taking surrogate data as input and estimating motion as an output. A full review of motion models is provided by McClelland et al (2013). Model coefficients are calculated from a set of training data taking surrogate and motion data as inputs. These model coefficients allow motion estimation from new surrogate data when motion data cannot be captured. Formally, a direct correspondence model can be written as: where M is a motion estimate, φ is the model, and s is a vector containing input surrogate data. The model is made of a matrix A of scalar coefficients, the size of which depends on the number of surrogates used and the number of motion parameters (i.e. the sizes of s and M).
The motion model we present in this paper consists of only PET and MR data as input. We test 4 different correspondence models, limited to only 1 or 2 surrogate signals and either a linear or polynomial system up to an order of 2. Higher number of surrogate signals and higher model complexity lead to a larger coefficient matrix A, and this can lead to over-fitting of data and increased risk of extrapolation errors.

PET model surrogate
The surrogate signal for our motion model is a respiratory signal extracted using PCA from the raw PET data (Thielemans et al 2011); previously shown to provide similar information to an MR pencil-beam navigator to track breathing (Manber et al 2015). PET list-mode data are unlisted into short low spatial resolution sinogram frames. Sinograms are then spatially smoothed, a scale factor is applied to account for tracer kinetics, the Freeman-Tukey transformation is applied to approximately convert Poisson noise to Gaussian, then finally PCA is executed. With PCA, each sinogram in the series is approximated as where s t ( ) is the sinogram acquired at time t, t is the time at the middle of the frame, s is the mean of all sinograms, c k is principal component (PC) k and w t k ( ) is the scalar weight factor for sinogram s t ( ), PC k. For one PC, each sinogram in the time series therefore has a single weight factor, calculated as the voxel-wise multiplication of the PC with the difference between the sinogram and sinogram mean, We define our PET-derived respiratory signal as the weights corresponding to Principal Component k = 1: In this work, we define the gradient of our respiratory signal to be used as a derived surrogate, as

MR model imaging data
For the imaging input to the motion model, we use a fast 2D multi-slice gradient echo (GRE) sequence, consisting of sagittal slices at L slice locations, covering the thorax and abdomen (including lungs, liver, pancreas etc). The sagittal imaging plane was chosen, assuming through-slice motion to be minimal. Previous studies report lung lesion displacements of only 1.2 mm laterally, compared to 2.2 mm and 5.5 mm in the anterior-posterior and superior-inferior directions respectively (Seppenwoolde et al 2002).
where vector t l represents the acquisition times of slices at slice location l. At each slice location, one reference image I r l , is chosen by binning all images based on the value of the respiratory signal P t ( ) to find a single dynamic slice at the 'exhale' position. Non-rigid registration is used to find voxel-wise deformation fields d t l ( ), that map each , → ( ) are found by separate registrations. Registrations were performed with MIRT (Medical Image Registration Toolbox) (Myronenko and Song 2010) in Matlab (Mathworks, Inc.). The Residual Complexity cost function within the toolbox was chosen as it provides a degree of robustness to vessel intensity changes and spatial intensity distortions that were present in the GRE MR images. The following parameters were used: b-spline control point grid size = 5 voxels, λ (regularization weight) = 0.05, and α (trade off parameter controlling the sparseness of coefficients) = 0.01. For cases with very large deformations, lambda was reduced to 0.01.

Joint PET-MR motion model
In the case of our PET-MR motion model, we use our PET-derived signals P t ( ) and P t ( ) ′ as potential model surrogates. Using the gradient as the second surrogate allows us to capture inter-and intra-cycle variation (hysteresis). Our four different models are therefore: • 1 surrogate 2nd order (1D-poly) • 2 surrogates 2nd order (2D-poly) We now explain the model construction using the two surrogate linear model as illustration. Training data comprising of PET-derived signals P t ( ) and P t ( ) ′ , and MR-derived deformation fields d t l ( ), are used as input to the model. The d t l ( ) for each value of t l are put into column vectors containing x and y vector components at all We then form a matrix D t l ( ) for each slice location l by placing each vectorised d t l ( ) from one slice location next to each other, such that: ,..., is the MR frame number at slice location l.
We search for model coefficients that minimise the least squares fitting error to the training data: where A l are the independent model coefficient matrices for each slice location l L 1... ∈ . P t l ( ) and P t l ( ) ′ are the values of the signals P t ( ) and P t ( ) ′ at the times of each d t l ( ) in D t l ( ), found with linear interpolation. The vector 1 is a vector of 1s of the same length as P t l ( ) and P t l ( ) ′ . We also want to account for a constant unknown time shift δ between the PET and MR system clocks, so this is included as an unknown in the functional: The minimisation is solved to find each A l and δ in one optimisation with a robust multilinear regression scheme with an iteratively weighted least squares cost function (Holland and Welsch 2007). This robustness ensures a good fit to the data by assigning d t l ( ) a lower weight in the optimisation if the registration failed.
New deformation fields d l can therefore be estimated at each slice location using the model with new input surrogate values p and p′ and the model coefficient matrices Deformations for between-slice voxels are found using linear interpolation to form one 3D deformation field (the 3rd vector orthogonal to the sagittal plane is set to zero) in the PET image space per gate, using PET and MR Dicom geometry information.

Method
In this paper we present a methodology to build a patient specific PET-MR motion model by extending a 3 min scan to 4 min, with no change to the protocol for the original 3 min. Workflow steps are outlined below, and presented in figure 1 and explained in detail below: 1. Acquire 4 min of PET-MR data. 2. Extract surrogate signals from raw PET data. 3. Bin GRE MR slices based on signal to identify an 'exhale' reference image at each slice position and perform non-rigid registration of all MR slices. 4. Form motion model, providing a correspondence between deformation fields and surrogate signals. 5. Form trajectory map and plan PET data gating scheme. 6. Unlist raw PET data into several gates according to step 5. 7. Calculate 3D deformation fields for each gate with motion model. 8. Warp μ-map with motion model to form gated μ-maps. 9. Motion compensated image reconstruction (MCIR).
Details of all steps are outlined in this section.

Data acquisition
On patients undertaking clinical whole-body PET-MR scans, a 3 min PET-MR acquisition was extended to 4 min to acquire data to build the motion model, whilst acquiring extra MR during the first 3 min for model validation. The acquisition consisted of: • PET list-mode (4 min), in same FOV as one previous PET bed positions.
Scan parameters: slice thickness 10 mm, gap between slice centres 25 mm, repetition time 5.1 ms, echo time 2.5 ms, flip angle 10°, pixel bandwidth 965 Hz, matrix size 192 144 × , FOV 262 349 × mm, in-plane resolution 1.8 1.8 mm 2 × , IPAT 3. Only the last minute of the acquisition are used to form the motion models (model slices) and the rest are for validation purposes (test slices). A temporal resolution of 300 ms ensures a good temporal-spatial resolution trade off, where intra-frame motion is minimised and interframe motion is captured. The sequence was optimised to obtain good contrast in all anatomies of interest such as lung, liver and pancreas.
All data were acquired using an integrated 3T PET/MR system (Biograph mMR, Siemens Healthcare). Additional data were acquired as part of calibration and service development protocols. Patients consented to the use of their data for research purposes.
3.1.1. Calculation of gate size. Each motion model can produce a deformation field estimate for any value of the surrogate, but computational effort for motion-compensated listmode reconstruction is high due to deformations varying at each time point, therefore, we binned the PET data which have similar surrogate values (and hence deformation fields). We use a novel efficient binning scheme to group PET data based on the amplitude and range of breathing throughout each acquisition to ensure that intra-bin motion is less than half the full-widthhalf-maximum (FWHM) resolution of the scanner. A trajectory map of P and P′ is used to track the surrogate signals throughout the scan. First, a grid is drawn across the trajectory map ( figure 2(a)). Each outer square on the grid following the perimeter of the trajectory is marked ( figure 2(b)) and deformation fields at every slice location are estimated for each of these squares. For each row (constant value of P′) and each column (constant value of P) the deformation range is found for each voxel, between the 2 outermost squares (figure 2(c)). The deformation range with the highest sum over all voxels in both P and P′ directions over all slice locations is found ( figure 2(d)). In the 2 directions, voxel deformations are histogrammed (figure 2(e)). The maximum respiratory displacement is marked as the value at the 98th percentile in each of these directions. The percentile is used to disregard any local outliers in the deformation fields. Bin size in both directions is then found by dividing the maximum respiratory displacement by half the FWHM resolution of the scanner (4.5 mm). The resulting gating scheme is shown in figure 2(f), with gates marked on the trajectory.

Motion compensated image reconstruction.
Motion compensated image reconstruction (MCIR) was used to form motion corrected PET images with randoms, scatter processes and motion incorporated in the system matrix of the reconstruction (Polycarpou et al 2012, Tsoumpas et al 2013. An ordered subset expectation maximisation reconstruction algorithm was used in the STIR software, with 21 subsets, 3 iterations and 4 mm Gaussian post-filtering. Each PET acquisition was reconstructed without motion correction, and with motion correction 4 different ways using the 4 motion models. All reconstructions were carried out both with attenuation correction (AC) and without attenuation correction (NAC). NAC PET reconstructions allow examination of how well the motion correction alone is performing, without added changes due to attenuation correction. PET data processing (unlisting, reconstruction etc) was carried out with STIR (software for tomographic image reconstruction) (Thielemans et al 2012).  For each MCIR reconstruction, a stack of emission sinograms, 3D voxel-wise deformation field and background sinogram are required for each gate. Using the gate size information found with the method described above, gating was then applied to the trajectory map and raw PET data were unlisted into one stack of sinograms per gate. Deformation fields are estimated with the motion model using the surrogate values at the centre of each gate according to figure 2(f).
The attenuation μ-map is warped to the exhale position with deformation fields estimated by the motion model, using the values of the surrogate signals P and P′ during the MRAC sequence acquisition. This exhale μ-map is then warped to match each of the gates. Randoms and scatter sinograms were also calculated for each gate, and together with the μ-map at each gate, make up the background sinograms.

Evaluation of motion models on MR data
All images I in the GRE MR series (test slices and model slices) were registered to their slice location reference image I r to provide warped images I w and associated deformation field d, using the same non-rigid registration scheme as described previously. We assume these deformation fields provide a 'perfect' registration and call these the gold-standard. Deformation fields d were estimated using each type of model. Each image I can therefore also be warped to an estimated 'exhale' position to form image I ŵ .
We tested how well each model performed by comparing images I w to I ŵ with mutual information (MI) and sum of squared difference (SSD) metrics, and comparing deformation fields d to d by calculating Euclidean distance (ED). All metrics were used to compare how each model performed regarding model-fit error (testing the same set of GRE images that were used to form the model) and model-prediction error (testing using the set of GRE images not used as input to the model-acquired prior to the last 1 min model-building section).
We define a performance index M ( ) ξ of an image-based metric M as percentage improvement by dividing the improvement from model-warping by improvement from registration-warping: where n is the frame number in the set to be evaluated-either test slices or model slices, and M is either metric MI or SSD. A high value of 100% of index M ( ) ξ shows the model performs as well as registration in predicting the correct warping of all images back to their respective exhale reference images.
When using the ED metric, we compare model-based deformation fields with registrationbased deformation fields. We calculate the voxel-wise ED sum between each d and d as absolute values in mm: where V is the total number of voxels in all frames in the set. A Wilcoxon signed rank test was used to test the statistical significance of the results, with a threshold of 0.0001 (99.99% confidence), so p-values below this were considered significant.

Evaluation of motion models on PET data
A motion-compensated reconstruction of the 4 min PET acquisition was carried out using deformation fields estimated by each of the 4 motion models. This was done for 5 patient data sets, all with suspected/known lesions in the liver, pancreas or lung-patient A (4 pancreas), B (6 liver), C (10 liver), D (4 pancreas), E (6 liver/lung).
Motion-corrected images were compared to uncorrected images visually, and quantitatively with SUV measures in a region of interest (ROI) containing areas of high tracer uptake in AC and NAC PET. Measures used are SUV max , and SUV peak , defined as the maximum average activity concentration within a 12 mm diameter sphere inside the ROI (Boellaard et al  Figure 3 shows a dynamic line profile through the lung and liver (including lung vessels) for original uncorrected images I, registered images I w and model-warped images I ŵ for the 2D model in one patient. Figure 4 shows model performance for each patient in terms of MI for both model-fit and model-prediction. For model-fit error, the 2D-poly model performed best in terms of MI, SSD and ED, in 98%, 93% and 100% of patients respectively. For model-prediction error, the 2D (linear) model performed best in terms of MI, SSD and ED, in 82%, 78% and 65% of patients respectively. Figure 5 shows average performance index for all metrics over all 45 patients, in terms of both model-fit (figure 5(a)) and model-prediction ( figure 5(b)). Overall, for model-fit, the 2D-poly model performed best, and for model-prediction, the 2D (linear) model performed best, and these were statistically significant (p < 0.0001) when comparing with each other model, for each of the 3 metrics. Figure 6 summarises the mean changes in SUV peak and SUV max in lesions in the 5 patient data sets. In all patients, average SUV peak and SUV max have significantly increased due to motion correction in both the AC and NAC PET images when comparing uncorrected and motioncorrected images, but no model performs significantly better than any other.

Evaluation of motion models on PET data
Visually, the motion-corrected images from the four types of motion model were comparable, with all lesions showing an increase in sharpness. Figures 7 and 8 show examples of the improvements in image quality apparent in the motion-corrected PET reconstructions, using the 2D model. For patient A (figure 7), the blurring in the lesion marked in the uncorrected image has been reduced in the motion-corrected image and the lesion appears sharper.
For patient B (figure 8), the MRAC was acquired at a deep inhale position (possibly due to patient misunderstanding the instructions or easier due to illness), causing a mismatch between the attenuation μ-map and PET emission data. This in turn causes artefacts in the uncorrected images above the liver. In this example, the combination of misaligned μ-map, respiratory motion and avid lesions present at the top of the liver causes them to be almost invisible in the uncorrected image. Two lesions are marked on the sagittal image, and one of them is also marked on the coronal image. Both lesions are visible on the motion-corrected image. Both lesions, along with others that appear in the motion-corrected image were confirmed to be present in the patient using extra MR and CT information confirmed by an accredited radiologist.

Discussion
We have demonstrated the feasibility of building and applying a motion model with MR and PET data only; using only 1 min of data with no external devices, to build a correspondence model capable of estimating respiratory motion. For a whole body clinical PET-MR examination at our institution, PET data at each bed position are acquired for 3 min. At least 2 of these positions (thorax and abdomen) will be adversely affected by respiratory motion, and for these 2 positions, the 3 min scan can be extended to 4 min with no change needed to the protocol for the original 3 min. We explored various ways to set up the model with different numbers of surrogate signals, and different levels of complexity, and tested each model on 45 MR patient data sets, and then on 5 identified oncology patients. The methodology presented builds a robust continuous motion model with the capability of interpolation and extrapolation, by predicting motion outside of what was visible in the model data, thus allowing use of 100% of PET data.
The overall results from all metrics indicated that the 2D-poly model performed the best in terms of model-fit ( figure 5(a)). This was to be expected, as the larger the model coefficient matrix, the better the model can fit to the data. The aim of this work was to use the model to predict motion so other clinical MR sequences can be acquired in parallel with the PET scan, so the results of the model-prediction are of more interest. The general trend shows that performance drops when using the models to predict motion ( figure 5(b)). These tests also showed that the simpler 2D model (linear) performs best in terms of motion prediction. This model is able to capture inter-cycle variability (hysteresis) as well as exhale-inhale motion.
When applied to PET data, the models all showed sizeable and statistically significant improvements by using MCIR compared to a single uncorrected reconstruction, but there  were no statistically significant differences in the extent of improvement between different models. It could be assumed that the quantitative evaluation on MR data gives a much better indication of model performance as 3 different global metrics were used, as opposed to only ROI based metrics in the PET. MR images are also of a much higher spatial resolution, so small discrepancies between the models would be more readily testable. It should be noted that for the evaluation on PET data, we did not differentiate between model-fit and modelestimation as the whole 4 min PET acquisition was used in the reconstruction. However, as 3 of the 4 min of PET data were acquired outside of the model-building section, the reconstructions tested mainly model-prediction accuracy.
The difference in performance of model-fit and model-prediction in the patient-by-patient MI results of figure 5 give clues as to where the model can fail. For example, there is a big drop in MI ( ) ξ in all models for patient 5 in model-prediction compared to model-fit. When examining the dynamic MR for this patient, it is clear that bulk motion occurred during the first half of the PET scan. Currently our method cannot account for this type of bulk motion, but it could potentially be captured with intermittent MR-based checks throughout clinical scanning. Other work such as (Kolbitsch et al 2014) explores the problem of stand-alone bulk motion (without respiratory motion) in the PET-MR context.
It is clear from the PET images for patient B (figure 8) that a large source of error in the uncorrected images of this patient stemmed from the μ-map being mis-matched with the emission data. This is highlighted by the large increases in SUV peak and SUV max in the AC motion-corrected images for this patient (figures 6(a) and (b)), whereas the increase is much smaller in the NAC images. The extent of the μ-map mismatch can be seen by looking at the original signal P for this patient, where it is clear that a deep inhalation was taken at the start of the scan where the MRAC was acquired, as visible in figure 9. This shows the mis-match between an exhale reference image and the corresponding slice in the MRAC. Also shown is the alignment of the MRAC after warping to the exhale reference image. Here, the warp applied to the MRAC is outside of the model-building phase, nevertheless, the warped MRAC matches the exhale reference images reasonably well. However, more work is needed to test the extrapolation performance of the models quantitatively. This result suggests shows that motion artefacts caused by the MRAC being acquired at the incorrect respiratory position can be reduced with the proposed methodology. This could be useful for patients with lung problems that find it difficult to hold their breath at exhale, or those with communication difficulties.
There are a number of limitations to the study. Firstly, there are some confounding factors in interpreting the results from the metrics used on MR data. We assume that the nonrigid registrations worked perfectly to provide our motion ground truth, but this was not the case. Occasional registration errors may occur, but the use of the robust model means these stand-alone errors should not affect the model performance. A general, easy to use, opensource registration scheme was utilised, which we chose as a practical method, but it cannot deal with non-diffeomorphic transformations (sliding motion), as would occur between the liver and the ribs. This could be overcome by utilising a registration method that allows for non-smooth deformation fields at organ boundaries where sliding motion occurs, providing a piecewise-diffeomorphic deformation field (Risser et al 2013). Secondly, by only examining motion in the sagittal plane, we assume no lateral motion during respiration, and interpolation is used for deformations between slices. Although we did not explicitly enforce continuity of the model coefficients across slices, we did not observe any abrupt changes, which provides extra confidence in the ability of the motion model, driven by the PET-derived signal, to capture the actual patient motion. We chose a 2D multi-slice acquisition scheme to ensure good spatial resolution, but we also use a gap between slices to speed keep scan time to a minimum. A full 3D motion model including lateral motion could be found by acquiring contiguous slices (Würslin et al 2013, McClelland et al 2014 but this would extend the scan time.
Overall, results indicate that the 2D linear PET-MR motion model performs best when predicting motion as measured by MR data, and is effective in PET motion correction, suggesting that in a clinical setting the 2D linear model would be the best choice, although evaluation on a larger patient cohort is required. PET image results of patient B indicate the importance of proper motion correction through μ-map misalignment as it can lead to improved lesion detectability. For future validation work, a quantitative assessment of lesion detectability would also be of interest.

Conclusion
We have demonstrated a methodology to use only 1 min of PET and MR data to build a patientspecific respiratory motion model, capable of predicting motion at time periods outside of the model. We tested different complexities of models on MR data, demonstrated improvements in clinical PET images when the motion information was incorporated into a PET reconstruction, and showed that lesion detectability could potentially be increased. The methodology built on previous work, to now allow interpolation and extrapolation in a continuous motion model, thus allowing more efficient use of PET data.