Evaluation of a direct motion estimation/correction method in respiratory‐gated PET/MRI with motion‐adjusted attenuation

Purpose Respiratory motion compensation in PET/CT and PET/MRI is essential as motion is a source of image degradation (motion blur, attenuation artifacts). In previous work, we developed a direct method for joint image reconstruction/motion estimation (JRM) for attenuation‐corrected (AC) respiratory‐gated PET, which uses a single attenuation‐map (μ‐map). This approach was successfully implemented for respiratory‐gated PET/CT, but since it relied on an accurate μ‐map for motion estimation, the question of its applicability in PET/MRI is open. The purpose of this work is to investigate the feasibility of JRM in PET/MRI and to assess the robustness of the motion estimation when a degraded μ‐map is used. Methods We performed a series of JRM reconstructions from simulated PET data using a range of simulated Dixon MRI sequence derived μ‐maps with wrong attenuation values in the lungs, from −100% (no attenuation) to +100% (double attenuation), as well as truncated arms. We compared the estimated motions with the one obtained from JRM in ideal conditions (no noise, true μ‐map as an input). We also applied JRM on 4 patient datasets of the chest, 3 of them containing hot lesions. Patient list‐mode data were gated using a principal component analysis method. We compared SUVmax values of the JRM reconstructed activity images and non motion‐corrected images. We also assessed the estimated motion fields by comparing the deformed JRM‐reconstructed activity with individually non‐AC reconstructed gates. Results Experiments on simulated data showed that JRM‐motion estimation is robust to μ‐map degradation in the sense that it produces motion fields similar to the ones obtained when using the true μ‐map, regardless of the attenuation errors in the lungs (< 0.5% mean absolute difference with the reference motion field). When using a μ‐map with truncated arms, JRM estimates a motion field that stretches the μ‐map in order to match the projection data. Results on patient datasets showed that using JRM improves the SUVmax values of hot lesions significantly and suppresses motion blur. When the estimated motion fields are applied to the reconstructed activity, the deformed images are geometrically similar to the non‐AC individually reconstructed gates. Conclusion Motion estimation by JRM is robust to variation of the attenuation values in the lungs. JRM successfully compensates for motion when applied to PET/MRI clinical datasets. It provides a potential alternative to existing methods where the motion fields are pre‐estimated from separate MRI measurements.


INTRODUCTION
Patient respiratory motion is a source of quantitation errors in positron emission tomography (PET), due to the degradation of the image resolution 1 and potential misalignment with the attenuation map (l-map) used for attenuation correction (AC), which is derived from computed tomography 2 (CT) or magnetic resonance imaging (MRI). 3 Both phenomena can compromise accurate detection and quantitation of lesions in the reconstructed PET. 4 To reduce the effect of motion, motion-free images can be reconstructed from a single gate corresponding to one respiratory phase, 5 but the reduction in the number of counts increases noise and bias. The reconstructed gates can be coregistered and averaged to a single image to reduce noise, but the final output suffers from bias induced by noise from each gate. 6,7 Ideally, motion correction (MC) is achieved using an estimated respiratory motion field directly incorporated into the reconstruction in order to maximize the number of usable counts. Most MC approaches in PET can be classified into two categories: indirect and direct.
Indirect approaches consists of pre-estimating a motion field that is then incorporated in the PET system matrix for MC iterative reconstruction. For example, the motion can be derived by registering non-AC-gated PET images, 8,9 but registration suffers from noise and low contrast from the individually non-AC-reconstructed gates. Alternatively, the motion can be pre-estimated from a separate dynamic MRI sequence 4,10,11 or gated CT, 12 but they require additional MRI or CT measurements.
Direct approaches consist of jointly estimating the activity distribution and the motion field directly from the PET data, without reconstructing "intermediary" PET images. A single reference activity image is estimated together with the motion by performing penalized maximum-likelihood (PML) of the PET data. Such approaches have been proposed in the literature [13][14][15][16] but they ignored attenuation. Recently, we developed a probabilistic model where the activity distribution is deformed alongside the l-map, 17,18 which produces a motion-free AC-reconstructed activity image. Moreover, this approach incorporates l-map misalignments within the motion so that the final output activity image is free of attenuation artifacts. This approach, namely joint activity reconstruction/motion estimation (JRM), was successfully applied on PET/CT patient data. 17 In principle, JRM requires an accurate forward model, and in particular an estimate of the attenuation that can be deformed to fit each respiratory gate. This raises the question of its applicability in PET/MR. The standard protocol for the thorax is to perform AC with a l-map derived from a segmented Dixon MRI sequence, [19][20][21] with generic attenuation values allocated to each of the segmented classes (generally air and soft adipose tissue). When the input l-values do not correspond to the patient values, the JRM optimization algorithm will attempt to match PET projection data that are outside of the range of the forward model, which can affect motion estimation. To our knowledge, a direct motion compensation approach for AC reconstruction with a single lmap has never been used in PET/MRI.
In this paper, we evaluate the impact of a deteriorated MRI-derived l-map on motion estimation by JRM with a non-time-of-flight (TOF) PET/MRI system. We first analyzed the existence of a solution to the JRM problem with a locally deteriorated l-map. We performed XCAT simulations to generate a gold-standard, and applied JRM using simulated MRI-derived l-maps with a range of deterioration (wrong lung and bone l-values and truncated arms). We then applied JRM on four sets of respiratory-gated PET/MRI clinical data, and investigated the effects of JRM motion compensation on hot lesions, as well as comparing the JRM reconstructed images with non-AC individually reconstructed gates.

2.A.1. Joint image reconstruction/motion estimation with attenuation correction in gated PET
Positron emission tomography gating is achieved by regrouping the raw list-mode data into n g sinograms fg ' g¼fg 1 ; ...; g n g g, each sinogram vector g ' 2 R n b þ comprising n b entries (detector bins pairs) and corresponding to a single gate on which the patient is assumed to be static (i.e., no motion). At each gate ' the collected data g ' results from the emission of the radiopharmaceutical tracer, modeled by an activity concentration volume f ' 2 R n v þ comprising n v voxels, which represents the tracer distribution at gate '. The probabilities of detecting of the annihilations are altered by the presence of an attenuation medium (patient tissues, bones, the bed, etc.), also distributed in a volume l ' 2 R n v þ . Each entry [f ' ] j and [l ' ] j , respectively, corresponds to the activity concentration and the attenuation at voxel j and gate '. At this stage, the gated activity distribution volumes f ' can be reconstructed from the corresponding data g ' , accounting for the corresponding attenuation l ' , with the help of iterative algorithms such as maximum-likelihood expectation-maximization (MLEM) 22,23 and penalized-MLEM. 24 In Bousse et al. 17 we adopted a model where the activity f ' and attenuation l ' at each gate ' are produced by the deformation of two volumes f and l: (1) f(r) and l(r) being interpolated versions of f and l respectively, and r j the center of voxel j. This simultaneous warping activity/attenuation model was used for joint image reconstruction/motion estimation (JRM) with attenuation correction by maximization of the penalized log-likelihood Φ with respect to the activity image f and the motion fields Uðf ; fu ' g; fg ' g; lÞ; (2) Uðf ; fu ' g; fg ' g; lÞ , X n g '¼1 Kðg ' j g ' ðf ; u ' ; lÞÞ À Uðf Þ À Vðfu ' gÞ; where Kðgj gÞ , P n b i¼1 g i log g i À g i is the Poisson log-likelihood (the similarity term), which depends on the gated data g ' and on the expected PET data g ' ðf ; u ' ; lÞ defined as Medical Physics, 44 (6), June 2017 with a(l) ≜ diag{e ÀRl }, and U and V are smoothness penalty terms on the activity image and the motion fields, respectively. In (4), P 2 R n b Ân v represents the nonattenuated, non-TOF PET system matrix, that is, defined by the system geometry and detector sensitivity, R 2 R n b Ân v is an operator that computes the linear attenuation along each line of response, s ' is the acquisition time corresponding to gate ', and b ' 2 R n b is a background vector which is the sum of the expected random and scatter counts.
In Bousse et al., 17 solving (2) was achieved by alternating the optimization in f and {φ ' }, with the help of a MC-MLEM algorithm (using the current motion estimate) for f and a quasi-Newton optimization 25 for {φ ' }. A single imagef is reconstructed using the entire gated data {g ' }, with MC using the estimated deformation fields fû ' g. This reconstructed image does not correspond to a reference gate. The reconstructed gated activity images at each gate ' aref ' ,Ŵ 'f , whereŴ ' is the warping operator derived fromû ' and are considered as the final output.
In this approach, the attenuation map l is assumed to be known. A solution of (2) is achieved with an activity imagef and motion fieldsû ' such that the expected projections g ' ðf ;û ' ; lÞ match the gated data g ' . In particular, the warped volumesf ' ¼Ŵ 'f andl ' ¼Ŵ ' l are aligned with the data.

2.A.2. The effect of a "Wrong" l-map input on motion estimation with JRM
In PET/MRI, the usual procedure is to utilize an MRI Dixon sequence segmented into three classes 19 (air, soft adipose tissue) or four classes 20 (air, soft and adipose tissue and bone). Generic (nonpatient dependent) values are attributed to each class to form a "piece-wise constant" l-map, used for attenuation correction. However, it has been observed that the lung attenuation may vary significantly across the respiratory cycle 26 (AE20%). In addition, MRI images may suffer from truncation due to the smaller field of view. 21 In the forward model (4), the attenuation map l is a fixed input image. The deformation fields {φ ' } can change its shape but not its values. This can be potentially problematic when the values of the input attenuation l are incorrect.
We now analyze the existence of a solution to JRM in non-TOF PET, for a single gate case (the multiple gate case can be similarly addressed as in Bousse et al., 17 ) when the input attenuation values are wrong. More precisely, we will demonstrate that with a locally deteriorated input l-map it is possible to find an approximate solution to the JRM problem (2) with the correct motion.
Assume the true activity and attenuation are respectively f H and l H , and the input l-map for attenuation correction is l. For a single-gate noiseless reconstruction problem, in the absence of penalty terms, the JRM task (2) reduces to finding an activity image f and a warping operator W [defined by some diffeomorphism φ as in (1)] that fit the observed data

Case 1
When the true attenuation l H is a deformed version of the input attenuation l, that is,Wl ¼ l H for some invertible warping operatorW, a solution to (5) is achieved with With this solution, the reconstructed activity f is "aligned" with the input attenuation l, and the estimated motion W realigns f and l to the observed data. 17

Case 2
On the other hand, if l H is not a deformed version of the input attenuation l (which happens for example when the values of the input attenuation l are incorrect), then l H À Wl 6 ¼ 0 for all W. Equation (6) indicates that a(l H À WlÞPf H should be in the range of P. In general, this will not be possible. However, there are some examples where this is feasible, such as when aðl H À Wl) is a multiple of the identity matrix. This condition is satisfied if and only if l H À Wl forward projects to a uniform sinogram, which is generally untrue.

Case 3
We now assume that input attenuation l satisfies for some warping operatorW and an error g with ‖g‖ 1 small, which corresponds to a situation where the input attenuation is a deformed version of the truth, with an additive error (local deterioration). It is possible to find an approximate solution to (5), along the lines of previous work from Thielemans et al., 27 which investigated the influence of using the wrong attenuation map in image reconstruction. We seek for a solution f,W to A first order Taylor expansion of (7) gives Now assume that R=P (which assumes the PET scanner is non-TOF) and g is supported on a small area located far from the edges of f H , then it is reasonable to assume 27 that Pf H is fairly constant on the lines of response intersecting the support of g and (8) becomes where q is the mean value of Pf H on supp(Pg). Taking W ¼W gives Thus, is an approximate solution to (9)-and therefore to (5)-with the correct motion. In addition, the reconstructed image f ¼ f H þ qg is identical to the one derived from Thielemans et al. 27 (without motion). This reasoning can be extended to the multigate case by combining the misalignmentW with the motions W 2 ; ...; W n g at the other gates. 17 In practice, attenuation errors g in MRI-derived l-maps do not have a small support (e.g., lungs), so that the existence of a solution of the form (10) is not guaranteed. In the next section we set out to evaluate this potential problem using MRI-derived l-maps with incorrect lung values.

2.B. Simulation study
The aim of this simulation study was to experimentally assess the effect of a wrong input l-map on the motion estimation by JRM. We investigated two types of defects: wrong l-values in the lungs and truncation of the arms (in addition to the absence of bones).

2.B.1. Simulated data
We simulated a sequence of five activity and attenuation images (1839183952 volumes, 3.125 mm 3 cubic voxels), f H ' and l H ' (Figs. 1(a) and 1(c), respectively) using the XCAT phantom software. Each volume corresponds to a respiratory gate (' = 1 and ' = 5 correspond to inspiration and expiration, respectively). Only respiratory motion was simulated (no cardiac motion), with a 3 cm amplitude diaphragm motion, and a 1.2 cm amplitude anterior/posterior motion. The volumes contain a hot lesion in the lung.
Noiseless respiratory-gated data were obtained by projecting each activity image f H ' with the corresponding attenuation coefficients, ...; 5; and noisy data were generated as a Poisson process, ...; 5: We used a uniform background for b ' (randoms and scatter). The PET system matrix P is a standard 2-dimensional "sliceby-slice" projector modeling a 5 mm FWHM point spread function for resolution and we assumed P = R. The total number of counts for {g ' } was 3910 9 , including background effects (accounting for 30% of the total counts).
We generated a collection of l-map mimicking Dixon MRIderived l-maps, consisting of an XCAT attenuation image (at inspiration, that is, corresponding to gate ' = 1) with only two classes (soft tissues and lungs, no bones), with different possible values for the lungs: from À100% of the original value (no attenuation) to +100% (double of the original value). Each of these l-maps is denotedl x , x 2 [À100,100]. Figure 2(b) showŝ l 0 (x = 0), which has the correct lung l-values [i.e., same value as in the ground truth l H 1 ,F i g .2(a)]. In addition, an MRIderived l-map, denotedl trunc , was generated froml 0 by truncating the arms (Fig. 2(c)).

2.B.2. Reconstructions
Joint activity reconstruction/motion estimation was applied by maximizing the penalized-likelihood Φ (3) with respect to the activity image f and the deformation fields {φ ' }, with the attenuation maps described in the above section. The joint reconstructions (activity images and motion fields) are gÞ , arg max f ;fu ' g Uðf ; fu ' g; fg ' g;l trunc Þ: Note that solving (11) will produce errors inf x for large absolute values of x, but the aim of this work is to assess the robustness of the estimation of the motionû x ' . Due to the nonuniqueness of the solution in image registration (two motion fields can produce the same deformed image, especially in the presence of uniform regions), the estimated motion fields fû x ' g and fû trunc ' g were compared to the motion fields obtained from JRM in ideal conditions, using the true attenuation l H 1 at gate '=1 and the non-noisy sinograms fg H ' g: This approach shows how degraded l-maps may affect motion estimation by JRM. The motion fields {φ ' } were modeled using the same cubic B-spline linear combinations as in previous work, 17 with one control point every three voxels along each axis. We did not use a penalty term on the activity image (i.e., we used U = 0), but we included a "small" penalty term on the motion fields {φ ' }, consisting of a quadratic smoothing function with a small weight (c = 10 À2 ), similar to the one used in Bousse et al. 17 Each JRM reconstruction was run until convergence.
In the rest of this work,Ŵ x ' ,Ŵ trunc ' andŴ ideal ' denote the warping operators associated toû x ' ,û trunc ' andû ideal ' , respectively, as defined in (1). The corresponding reconstructed activity images and deformed l-maps at each gate ' are denoted l H 1 : Figure 1(b) shows the ideal reconstructed activity imageŝ f ideal ' at each gate '. As observed in Bousse et al., 17 they look similar to the true activity images f H ' (Fig. 1(a)). The deformed l-mapsl ideal ' are shown in Fig. 1(d). They also look similar to the ground truth l-maps l H ' (Fig. 1(c)). A non-MC reconstruction, reconstructed from the non-noisy data fg H ' g, is shown in Fig. 3.
Results were assessed by comparing the estimated motion fieldsû x ' andû ideal ' , both in the hot lesion and the entire volume, calculated as where r j 2 R 3 is the center of the j-th voxel and R ' is either the lesion's support at gate ' or the entire volume.

2.C. Patient data
We investigated four patient datasets, acquired for research purpose (following a clinical scan). Patients consented to the use of their data for research purposes. The PET data were acquired with a Siemens Biograph PET/MRI scanner 28 at University College London Hospital (UCLH), London, UK. Each acquisition consists of a 4 min list-mode PET-scan, single bed position, and an end-expiration Dixon MRI sequence for AC. The l-maps (Fig. 4) were derived from a segmented Dixon MRI sequence (air, soft, and adipose tissue), and the system included a maximum-likelihood for activity and attenuation (MLAA) algorithm 29 17 for a detailed description of the algorithm). Similarly to simulations, we did not use a penalty term on the activity image (i.e., we used U = 0), but only a "small" quadratic penalty term 17 on the motion fields (c = 0.2). For JRM, the list-mode PET data were regrouped into n g ¼ 5 respiratory gates following the strategy adopted by Thielemans et al.: 30 list-mode data were unlisted into 600 sinograms of 0.4 s duration. The dimension of each sinogram was reduced to n 0 b ( n b (n b is the original number of bins), in order to reduce computational time. We then performed principal component analysis on the resulting n 0 b Â 600 matrix in order to extract a frame by frame temporal respiratory signal, obtained by projection on the first principal component. This signal was then used to unlist the original list-mode file into n g ¼ 5 respiratory gates based on its amplitude. Contrary to simulations, gates ' = 1 and ' = 5 correspond to end-expiration and end-inspiration, respectively.
Projection and backprojection for MLEM-noMC and JRM, as well as background (scatter and random) estimation were performed using the open-source software for tomographic image reconstruction STIR 31 .A5 m m FWHM resolution model was used (for both P and R). For JRM, 5 background sinograms b ' (scatter and random), corresponding to the five gates, were derived by rescaling the total (ungated) background sinogram according to the gate durations s ' , '=1,...,5. The reconstructions are 289 9 289 9 127 volumes with 2.031 9 2.031 9 2.045 mm 3 voxels.
Although the activity values are different due to the incorrect l-values, the organ boundaries off x ' at each gate look similar to the ones of the ideal reconstructed activity imageŝ f ideal ' (Fig. 1(b)) regardless of x, suggesting that JRM motion estimation is not affected by the lung l-value defects. Similar results are observed with the deformed l-mapsl x ' (Figs. 6(a)-6(c)). Figure 7 shows the normalized absolute difference betweenû x ' andû ideal ' , as defined in (14), with x 2 {À100,À50,À20,0,20,50,100}. Results show, for the current choice of B-spline parametrization and penalty, that when the lung attenuation error does not exceed AE50% the difference is negligible, and does not exceed 0.45% with AE100% error. This shows that the impact of the deterioration the l-values on motion estimation by JRM is limited for reasonable errors on the l-map. Note that the minimum value is not reached by x = 0, possibly due to the noise in the data and/or the fact that the MRI-derived l-mapl 0 differs from the truth l H 1 (because of the absence of bones, the uniform soft tissues and lungs). Figure 8 shows the average reconstructed activity in the lungs versus the percentage of lung attenuation error x. One can observe that the lung activity increases linearly with the attenuation error, in a similar fashion to the approximate solution (10) we derived in Section 2.A.2. This shows that despite the nonlocal character of the lungs attenuation error, the reconstructed activity appears to be similar to the predictions.
To finish, we displayed the displacement fields vectors fr j Àû ' ðr j Þg in Fig. 9 (at gate ' = 5), for the ideal case (true l-map, no noise), MRI-derived l-mapsl x (with x = 0, À100 and 100), and with the truncated l-map (discussed in the next section). We observed that the displacement fields using MRI-derived l-maps (Figs. 9(b)-9(d)) are fairly similar to the one obtained from the ideal case ( Fig. 9(a)).

3.A.2. Truncated arms
Reconstructed activity volumesf trunc ' using the truncated MRI-derived l-map (Fig. 2(c)) by solving (12) are shown in Fig. 5(d). They appear similar to the ideal reconstructed activity imagesf ideal ' (Fig. 1(b)), and more importantly, the warped l-mapsf trunc ' ¼Ŵ trunc 'l trunc (Fig. 6(d)) appear similar tol 0 ' (Fig. 6(b)). This means that JRM estimated deformation fieldsû trunc ' that are not only matching the respiratory motion but also "stretching" the truncated arms in order to match the projected data {g ' }. This phenomenon can be observed on the displacement field ( Fig. 9(e)).l trunc and l trunc 1 ¼Ŵ trunc 1l trunc (gate '=1) are displayed in Fig. 10 for comparison. Similar results were observed in previous work, 17,18 where the activity image was reconstructed in the "input l-map space", then warped alongside the activity to match the observed PET data.

3.B. Patient data
For each patient, we display the images corresponding tô f 1 ¼Ŵ 1f , that is, the motion compensated reconstructed image at gate ' = 1 (end expiration), except Patient-C for which all gates are displayed.

3.B.1. Patient-A
Patient-A was acquired with a 68 Ga DOTA-TATE 132 MBq injection, 71 min prior to the scan. The patient was scanned for neuroendocrine tumor deposit in the liver. Figure 11 shows the transaxial and coronal slices of the reconstructed PET volumes with MLEM-noMC ( Fig. 11(a)) and FIG. 7. Mean absolute difference between fû x ' g and fû ideal ' g calculated following (14), for x varying between À100% and +100%.  (Fig. 11(b)). Although the lesion is visible on both reconstructions, it has a higher contrast on the JRM image and is less blurry. The SUV max analysis (Table I,r o w1 ) shows that JRM reconstruction produced an image with a 25.86% SUV max increase on the lesion in comparison to MLEM-noMC.

3.B.2. Patient-B and B
0 Patient-B and B 0 consists of a baseline and follow-up scan of the same patient, suffering from multiple endocrine neoplasia type 1 (pancreatic neuroendocrine tumors). The first scan was performed with a 68 Ga DOTA-TATE 150 MBq injection 144 min prior to the scan, and the follow-up scan was performed one year later with an 170 MBq injection, 61 min prior to the scan. The noise increase between the two scans may be attributed to the decay. Three lesions were investigated (denoted 1, 2, and 3). Figures 12 and 13 show respectively the transaxial and coronal slices of the reconstructed PET images from the baseline scan, using MLEM-noMC (Figs. 12(a) and 13(a)) and JRM (Figs. 12(b) and 13(b)). Each transaxial slice contains one hot lesion. The hot lesions are hardly visible on the MLEM-noMC images but they can be seen on the JRM images. Quantitative analysis (Table I, r o w s2 -4) shows that the lesions SUV max increased significantly when using JRM (from 16.37% to 60.71%).
This case study is particularly relevant, as on the noMC images only one tumor was detected by the clinicians, for which the default plan is a surgical intervention. However, the presence of three tumors (visible on the JRM images) changes the risk benefit balance of surgery, as the entire pancreas would need to be removed.  The follow-up scan reconstructions are shown in Figs. 14 and 15. The three lesions are now clearly visible on the MLEM-noMC images (Figs. 14(a) and 15(a)), but the contrast is greatly improved on JRM images (Figs. 14(b) and 15(b)), as compared to MLEM-noMC. SUV max analysis confirms this observation (Table I,

3.B.3. Patient-C
Patient-C was acquired with an 18 F-FDG 152 MBq injection, 129 min prior to the scan. The patient suffered from squamous cell carcinoma of the left axilla. No hot lesions were detected but the data suffer from high amplitude respiratory motion. The MLEM-noMC images ( Fig. 16(a)) are degraded due to the motion. In contrast, the JRM-reconstructed images (Fig. 16(b)) appear sharper and have higher contrast. More particularly, the shape of the JRM-reconstructed heart appears well defined, and appears similar to the non-AC gate 1 reconstruction (Fig. 17(a), first row), suggesting that the motion was successfully compensated.
Although the change in the images was not judged to be clinically relevant for this particular case study, 18 F-FDG is often used for assessment of myocardial viability and inflammation (such as myocarditis in sarcoidosis), and can benefit from the image quality improvement provided by JRM.  Figure 17(a) shows the non-AC-reconstructed PET images at each respiratory gate ' = 1,...,5 and the amplitude of the motion can be observed on the heart as well as on the right lung lower boundary (a green line was plotted for visualization purposes). In Bousse et al. 17 , we compared the warped l-mapsl ' ¼Ŵ ' l, derived from the JRM-estimated motion fieldsû ' , to the gated CT images (cine-CT 32 ). In absence of gated l-maps, we displayed the warped reconstructed imagesf ' ¼Ŵ 'f ( Fig. 17(b)), wherê f andŴ ' are respectively the single MC reconstructed image (using all the gates) and the estimated motion at gate ', and the warped MRI-derived l-mapsl ' ¼Ŵ ' l ( Fig. 17(c)), alongside the non-AC-gated PET images, at each gate '=1,...,5. Although both images are not directly comparable due to absence of AC in Fig. 17(a), we can observe that the deformed heart obtained withŴ ' appears similar to the one in the gated non-AC reconstruction. Similarly, we can observe that the lower boundary of the right lung in the non-AC-gated PET images seems to match withf ' andl ' , with perhaps the exception of gate 5 which was the shortest gate (hence the noisiest).
It is worth noting that for this patient data the arms were almost absent from the l-map (despite the MLAA correction), so the estimated motion fields were not able to stretch them down, as we observed with simulations in Section 3.A.2.

DISCUSSION
The quality of MC reconstruction is determined by the accuracy of the motion estimation. Although dynamic MRIderived motion models are under investigation, 4,10,11 their accuracy depends on the quality of the dynamic MRI, and they require modifications of the acquisition protocol (mainly addition of a dynamic MRI acquisition). No standard procedure has therefore emerged. On the contrary, JRM can be applied in principle to any PET/MRI protocols, although the motion estimation can be subject to noise (e.g., low counts gates). In addition, since JRM relies on PET data only,i ti s important to understand the effect of using an incorrect input l-map, which can be the case with an MRI-derived l-map. In order to fully assess the potential of JRM in PET/MRI, it should be compared with indirect MC approaches, for example, using motion fields pre-estimated from an MRI sequence. 4,10,11 This investigation showed that motion estimation with JRM appears not to be affected by the utilization of a MRIderived l-map. Simulation results presented in Section 3.A show that the estimated motion fields were similar to the ones obtained with the true l-map, used to generate the data. In the presence of large errors in the lung l-values, the estimated motion remained unchanged, as predicted by our analysis. In case of partial truncation of the arms, the estimated motion tried to match the PET projection data by stretching the arms. This result is not surprising since it was demonstrated in previous work 17,18 that JRM reconstructs an activity imagef aligned to the input l-map, and estimates a deformation field which warps both the reconstructed activity imagê f and the l-map in order to match the PET data. However, this result must be interpreted cautiously, as JRM cannot "recreate" missing features in the l-map, but only deform already existing features (c.f. the results on Patient-C, Section 3.B.3). For heavy arms truncation, an MLAA-type algorithm 29 should be applied to reconstruct the arms in the lmap.
Results on real data confirmed our observations from simulations. Hot lesions on JRM-reconstructed images are significantly more visible than on non-MC images (c.f. Patient- ): the SUV max values are higher when using JRM and do not suffer from blur due to motion. Similarly, the reconstruction of the heart (Patient-C) is also greatly improved when using JRM as compared with non-MC reconstruction. JRM-reconstructed images at each gate, obtained by applying the estimated motion operatorŴ ' to the MC reconstructed activityf returns images with features matching the ones of the gated non-AC reconstruction. These results show that the motion estimation from JRM is accurate despite the utilization of a segmented Dixon MRI sequencederived l-map.
We did not investigate the effects of a poorly aligned lmap (e.g., deep breath-in or deep breath-out) on patient data, as we showed in Bousse et al. 18 it requires a large number of iterations in non-time-of-flight PET. At the time of implementation, the STIR Matlab interface we used did not incorporate ordered subset EM (OSEM), and therefore such experiments were not feasible.

CONCLUSION
In this work, we implemented a direct AC-PET joint image reconstruction/motion estimation (JRM) for PET/MRI using a segmented Dixon MRI sequence-derived l-map for attenuation correction. As previous work showed that JRM was successful for PET/CT, 17 we showed here that it can also be used in PET/MRI, and that the motion estimation was still accurate even using a heavily degraded l-map. Future work includes: (a) assessing the potential of JRM to deal with a misaligned l-map in clinical PET/MRI, as we have previously done in PET/CT, 17 and, (b) comparison with indirect motion-estimation approaches.