Pulmonary imaging using respiratory motion compensated simultaneous PET/MR

Purpose: Pulmonary positron emission tomography (PET) imaging is confounded by blurring artifacts caused by respiratory motion. These artifacts degrade both image quality and quantitative accuracy. In this paper, the authors present a complete data acquisition and processing framework for respiratory motion compensated image reconstruction (MCIR) using simultaneous whole body PET / magnetic resonance (MR) and validate it through simulation and clinical patient studies. Methods: The authors have developed an MCIR framework based on maximum a posteriori or MAP estimation. For fast acquisition of high quality 4D MR images, the authors developed a novel Golden-angle RAdial Navigated Gradient Echo (GRANGE) pulse sequence and used it in conjunction with sparsity-enforcing k - t FOCUSS reconstruction. The authors use a 1D slice-projection navigator signal encapsulated within this pulse sequence along with a histogram-based gate assignment technique to retrospectively sort the MR and PET data into individual gates. The authors compute deformation fields for each gate via nonrigid registration. The deformation fields are incorporated into the PET data model as well as utilized for generating dynamic attenuation maps. The framework was validated using simulation studies on the 4D XCAT phantom and three clinical patient studies that were performed on the Biograph mMR, a simultaneous whole body PET / MR scanner. Results: The authors compared MCIR (MC) results with ungated (UG) and one-gate (OG) reconstruction results. The XCAT study revealed contrast-to-noise ratio (CNR) improvements for MC relative to UG in the range of 21%–107% for 14 mm diameter lung lesions and 39%–120% for 10 mm diameter lung lesions. A strategy for regularization parameter selection was proposed, validated using XCAT simulations, and applied to the clinical studies. The authors’ results show that the MC image yields 19%–190% increase in the CNR of high-intensity features of interest a ff ected by respiratory motion relative to UG and a 6%–51% increase relative to OG. Conclusions: Standalone MR is not the traditional choice for lung scans due to the low proton density, high magnetic susceptibility, and low T ∗ 2 relaxation time in the lungs. By developing and validating this PET / MR pulmonary imaging framework, the authors show that simultaneous PET / MR, unique in its capability of combining structural information from MR with functional information from PET, shows promise in pulmonary imaging. C 2015 American Association of Physicists in Medicine. dx.doi.org 10.1118


INTRODUCTION
Respiratory motion poses a significant challenge to pulmonary positron emission tomography (PET) imaging by causing blurring artifacts, which degrade both image quality and quantita-tion. Motion-induced blurring, which may lead to underestimation of tumor activity and overestimation of tumor volume, 1 is a major deterrent for therapeutic monitoring, which relies on the consistency of quantitative information across scans. Additionally, motion artifacts may have a tremendous impact on lesion detectability and often result in small lesions remaining undetected, 2,3 thereby jeopardizing diagnosis and staging. As PET plays an important clinical role in the staging and treatment evaluation of thoracic lesions, 4,5 such inaccuracies in quantitation hinder cancer management. A variety of motion tracking systems are available, which generate external chest motion trajectories (most commonly through visual markers monitored with cameras or through pressure sensors mounted on bellows) that can be used as surrogates for internal motion measures. 6 The respiratory signals generated by these systems can be used to perform respiratory gating, a technique in which emission events detected during short time windows corresponding only to a specific respiratory phase or gate (within which motion effects can be assumed to be negligible) are used for reconstruction. 7,8 This technique, however, sacrifices signal-to-noise ratio (SNR) as it only uses a small fraction of the emission events.
A more sophisticated solution to the motion problem is by means of determining volumetric deformation fields and using them to compensate for motion. Motion compensation via pre-reconstruction affine repositioning of the lines of response, which has been successfully applied to neuroimaging studies involving rigid motion of the head, 9,10 is unsuitable for pulmonary applications as respiratory motion causes non-rigid deformations of the thorax and upper abdomen. 11,12 Instead, accurate compensation of respiratory motion requires volumetric deformation fields computed through non-rigid registration typically of 4D anatomical images, gated in synchrony with the PET. The final PET image corresponding to a chosen reference gate is either computed via a reconstruct-transform-average (RTA) approach using individual gated PET images [13][14][15] or by incorporating the deformation fields directly into a motion compensated image reconstruction (MCIR) framework that outputs a single PET image at the reference gate. [16][17][18][19][20] Some works have taken these approaches one step further by performing joint estimation of images and deformation fields. [21][22][23][24] Studies comparing the RTA approach with MCIR generally concur that the latter approach yields higher quantitative accuracy. [25][26][27][28][29] While in principle it is possible to perform MCIR using motion information derived directly from 4D PET images, either without attenuation correction 30 or with attenuation correction using a respiratory-averaged computed tomography (CT) image, 31,32 the accuracy of these methods is limited by the poor anatomical detail in the PET images. 33 Hence, despite the added radiation dose imposed by 4D CT, the majority of MCIR studies reported are based on motion information derived from 4D CT images. 18,19,25,34 Due to the sequential nature of the PET and CT scans in PET/CT systems, accurate synchronization of PET and CT gates can be challenging and is further complicated by the inherent variability in the breathing patterns of patients over time. [35][36][37] Simultaneous whole body PET/MR in an emerging technology that unifies PET (a physiological imaging technique) with magnetic resonance (MR, primarily used for imaging the anatomy), thereby enabling synchronized capture of both structure and function. PET/MR successfully addresses several of the limitations of PET/CT in MCIR. Unlike CT, MR yields superb soft-tissue contrast and high spatial resolution without any radiation burden. Additionally, the simultaneity of the PET and MR components eliminates the problems due to mismatch between the PET and the CT in terms of position and respiratory phase frequently encountered with PET/CT. Despite the initial success of PET/MR in MCIR, [38][39][40][41][42] high quality pulmonary imaging with this new technology is still an open problem since the low proton density, large susceptibility, and very short T * 2 in the lungs make pulmonary MR imaging challenging. While tagged MR studies have been successfully used in phantoms, rabbits, and primates, [39][40][41] tags fade too quickly to capture the longer breathing cycle in humans. A variety of methods have been proposed to generate 4D MR image sequences in the presence of breathing motion. A slicestacking technique based on a 2D navigator was described in Ref. 43. 2D and 3D MR protocols were compared in Ref. 44. A navigator-based MR sequence with radial phase encoding was presented in Ref. 45 and was used in PET/MR simulation studies for nonrigid bulk motion correction. PET/MR phantom imaging results based on a real-time MR sequence with 1D navigators were presented in Ref. 46. These studies highlight the merits of MCIR in simulation and phantom studies.
In this paper, we present an MCIR framework for pulmonary PET/MR imaging and validate it using clinical datasets. This is an extension of our prior reports, 47,48 which, to our knowledge, present the first clinical results using PET/MR pulmonary imaging with MCIR. Two subsequent studies have shown preliminary clinical results demonstrating the utility of PET/MR in capturing and correcting for respiratory motion in the lungs. 49,50 4D MR images were assembled from 2D slices using a 1D pencil beam navigator and used to generate thoracic PET images using the RTA method in Ref. 49. In contrast with this work, which uses Cartesian sampling, we have designed a radial MR sequence, which is more wellsuited for 4D imaging applications. A self-gated 3D radial stack-of-stars MR sequence was presented in Ref. 50 and was used for RTA-based motion correction of PET images.
One key difference between this method and our approach is that we use a navigator, which, compared to self-gating, is less prone to corrupting factors such as noise, scanner gradient delay, field changes etc.
The main contributions of this work are as follows: We have designed a navigator-encapsulated golden-angle radial MR pulse sequence and combined it with sparsity-enforcing iterative reconstruction to generate high-quality 4D MR images of the lungs under free-breathing conditions. We have developed an optimized histogram-based data binning scheme for retrospective gating that leads to roughly equal amounts of PET and MR data per gate. We have also developed an MCIR framework for PET based on maximum a posteriori (MAP) estimation with gradient-based optimization. The advantages of regularization in MCIR have been demonstrated in Ref. 51 using the one-step-late expectation maximization or OSL-EM approach. Unlike OSL-EM, which has been shown to converge only for certain forms of the regularizer, 52 our MAP technique with gradient ascent guarantees convergence of the inverse problem for a wide range of priors.
Section 2 provides details on all the data processing steps, the key components of the MCIR framework, the simulation setup, the details of the clinical experiments, and the figures of merit used for evaluation. In Sec. 3, we present results from a simulation study and three clinical patient studies. A discussion of our results is presented in Sec. 4.

2.A. Overall framework
The data acquisition and processing framework shown in Fig. 1 was developed for a Biograph mMR whole body simultaneous PET/MR scanner (Siemens Medical Systems, Erlangen, Germany). The PET detector assembly in this scanner, consisting of arrays of MR-compatible avalanche photodiodes (APDs) and detector blocks with 4 × 4 × 20 mm 3 lutetium oxyorthosilicate (LSO) crystals, is housed in the bore of a 3 T superconducting magnet. The PET field of view (FOV) is 59.4 cm in the transaxial direction and 25.8 cm in the axial direction. For MCIR, PET data are acquired in list mode. MR k-space data are collected simultaneously with the PET. Both the PET and the MR raw data are binned to gates retrospectively. The gated MR images are used to compute deformation fields and gated attenuation maps, which serve as inputs to the PET reconstruction module. The MAP reconstruction module then outputs the final PET image.

2.B. PET data model and MAP reconstruction
The PET data model for MCIR can be obtained by concatenating the models for the individual gates and incorporating deformation information. The unknown activity image is assumed to correspond to a reference gate, a discrete phase of motion selected as the reference. Deformation fields map this reference gate to all other gates. The mean of the data vector corresponding to the kth gate,ȳ k ∈ R m , can be modeled as 53 y k = P sens P blur P k attn P geom T k x +r k +s k .
Here, x ∈ R n is the reference gate image, P sens ∈ R m×n is a diagonal matrix representing the detector sensitivity, P blur ∈ R m×n accounts for detector blurring due to scattering, random coincidences, and crystal penetration, P k attn ∈ R m×n is a diagonal matrix consisting of attenuation coefficients for the kth gate, P geom ∈ R m×n is the geometric projector, the (i, j)th element of which represents the probability that a photon pair produced in voxel j reaches the front faces of the detector pair i, T k is a trilinear interpolation matrix generated from voxelwise deformation fields, which transforms x to the kth gate, andr k ands k represent the mean of the random and scattered events for the kth gate. This model assumes motion effects to be negligible within each gate.
MAP estimation involves maximization of the posterior probability density function, p(x| y) = p( y |x)p(x)/p( y). 53 For PET imaging, p( y |x) is the Poisson likelihood with meanȳ (mean data vector concatenated over all gates). The prior probability is usually assumed to follow a Gibbs distribution of the form p(x) = (1/Z)e −βU(x) , where U(x) is a Gibbs energy function and β is the regularization parameter. Taking logarithms of both sides and dropping the term which is constant with respect to the unknown image, the MAP objective function can be represented as F. 1. Data acquisition and processing framework for MCIR using simultaneous PET/MR.
where L( y |x) is the Poisson log likelihood and U(x) is the uniform quadratic penalty, which represents a uniform Gaussian distribution on the nearest-neighbor intensity differences within the local neighborhood, N j , of a voxel j. In this paper, we compare three different image generation approaches: • OG: Conventional sinogram-based MAP reconstruction using counts for one gate-the reference gate (k = 1) with • UG: Conventional sinogram-based MAP reconstruction using all (ungated) counts with • MC: MAP MCIR using gated sinograms with In each case, we minimize the objective function in (2) using the preconditioned conjugate gradient algorithm with a bent Armijo line search to impose non-negativity. 54

2.C. MR image acquisition and reconstruction
MR imaging of the lungs is challenging due to the following facts: (1) The lungs contain only about 800 g of tissue and blood, distributed over a volume of 4-6 l, leading to very low proton density and hence low MR signal intensity, and (2) the high magnetic susceptibility at the alveolar surfaces leads to a very short T * 2 and hence causes rapid dephasing of the MR signal. 55 Additionally, thoracic motion necessitates the use of fast MR pulse sequences to prevent motion artifacts, further contributing to the challenge. For real-time motion capture, we have developed a novel MR pulse sequence entitled Golden-angle RAdial Navigated Gradient Echo (GRANGE). This sequence uses radial sampling, which enables oversampling of the central part of k-space (where most of the useful image information is typically concentrated) with gradual azimuthal undersampling of the extremities. Each radial line in k-space captures both low and high frequency information. A reduction in the number of k-lines therefore does not directly translate to a loss of resolution or ghosting artifacts unlike some Cartesian trajectories. Instead undersampling for radial sequences, as we will demonstrate, typically translates to radial streaking artifacts. Three key enabling characteristics of our approach described below, including two features of the GRANGE sequence and the use of sparsity-enforcing iterative reconstruction, make our approach ideal for pulmonary imaging with free breathing and retrospective data binning and processing.

2.C.1. Slice-projection navigator
For each repetition time (TR), our sequence first acquires a fast 1D slice-projection navigation signal to track the diaphragm position. 56 It then acquires radial k-space lines in a slice-interleaved fashion as illustrated in Fig. 2. The same sequence of k-lines is acquired from all the prescribed coronal slices in each TR. Each radial line is assigned to a respiratory bin based on its nearest navigator signal. If required, the spatial resolution can be further improved by interpolation. It is also possible to use an intrinsic signal from the DC point of the k-space lines for respiratory gating purposes. 57 However, a navigator is used in this work because it is less prone to corrupting factors such as noise, scanner gradient delay, field changes etc.

2.C.2. Golden-angle phase encoding
From one TR to the next, the radial k-lines are prescribed according to the golden angle formula: θ k = mod(2.399 96k, 2π). 58 As shown in Fig. 3, with phase encoding based on a fixed angular increment, a free breathing setting may lead to nonuniform and partial coverage of k-space. With golden-angle phase F. 2. GRANGE: A slice-interleaved golden-angle radial gradient echo pulse sequence with slice-projection navigator encapsulation. encoding, the coverage of k-space is more uniform even for an arbitrary number of k-lines, thereby imparting binning flexibility for gate assignment. We note that, while continuous golden-angle radial lines yield a more uniform distribution, the coverage might not be optimal when combining multiple respiratory cycles together. 59,60 However, this approach is still superior to sequential acquisition.

2.C.3. Iterative reconstruction
Binning of the k-lines into a larger number of gates causes progressive undersampling of k-space data per gate. Reconstruction using conventional filtered backprojection (FBP) leads to increasing levels of streaking artifacts as the number of gates is increased. Gridding reconstruction, another recon-struction algorithm widely used in radial MR reconstruction, suffers from a similar problem due to insufficient sampling. We adopted k-t space FOCal Underdetermined System Solver (k-t FOCUSS), a compressive sensing based reconstruction technique. 61,62 This iterative method commences by finding a low-resolution estimate and then prunes this solution to a sparse signal representation. This method takes advantage of the spatiotemporal sparsity of the images to reconstruct cine MR images from undersampled data.

2.D. Gate assignment
The binning of PET data into gates is usually either timebased or amplitude-based. 7 Time-based methods define gate indices by dividing the respiratory signal into equal-sized temporal bins. In contrast, amplitude-based methods define equal-sized bins along the amplitude axis. Time-based gating techniques fail when the breathing patterns are not perfectly regular. Several studies have shown that amplitude-based methods are more robust against inter-and intracycle variability of respiratory signals than time-based approaches. 7,63 However, traditional amplitude-based gating leads to varying numbers of event counts per gate. We have implemented a histogram-guided, amplitude-based gating strategy that incorporates information from MR navigator projections and assigns gate indices while ensuring that each gate has roughly the same amount of data. In comparison with traditional amplitude-based gating, where the gate assignment relies on equal-width amplitude bins, we use varying-width amplitude bins determined from histogram equalization. Figure 5 illustrates our gate assignment strategy. The navigator projections yield a set of 1D MR images along the superior-inferior direction, perpendicular to the diaphragm. Figure 5(A) shows a typical MR navigator image time series from a human. We extract a 1D correlation trace from this image as shown in Fig. 5(B) by computing Pearson's correlation coefficients with respect to a reference time point. This 1D trace is treated as our respiratory signal. The points on this trace are then sorted according to the correlation magnitude and divided into a chosen number of bins, such that each bin contains the same number of discrete data points. Figure 5(C) shows the corresponding equalized histogram, where the x-axis represents the fraction of data points in each bin and the y-axis indicates the corresponding unequal bin widths for the correlation values. The upper and lower limits for each bin [indicated using asterisks on the y-axis of Fig. 5(C) and using horizontal dotted lines in Fig. 5(B)] are used to quantize the correlation values [y-axis of Fig. 5(B)] to discrete gate indices. The dotted vertical lines (corresponding to the time points where the correlation trace intersects the horizontal lines) segment the time axis into separate gates. The result is shown in Fig. 5(D), which shows the division of the data acquisition time scale into a set of discrete gates. Since the PET and MR scans are synchronous, time stamps derived from the raw data can be used to align the time scales for the two modalities. The gate indices in Fig.  5(D) can, therefore, be directly used to bin the PET events into individual gates.

2.E. Deformation field computation and attenuation correction
Non-rigid registration based on diffeomorphic demons is used to compute the deformation fields mapping the reference gate MR image to the other gates. 64 In addition to being directly used by the MAP reconstruction module for PET image warping, the deformation fields are required for accurate attenuation correction of the PET images. Compton interactions with tissue scatter photons and diminish their energy leading to an attenuation effect which is accounted for in the data model in (1) as P k attn , a diagonal matrix containing survival probabilities representing the attenuation map. The attenuation correction factors for each gate are dependent on the phase of respiratory motion associated with that gate. 65,66 It was been shown that the use of either an average or a static attenuation map leads to prominent attenuation artifacts. 32,67, 68 We therefore use dynamic attenuation maps. The Biograph mMR generates separate attenuation maps for the scanner hardware dependent component and the patient or object dependent component. For our experiments, the latter is generated by the scanner using an MR Dixon scan acquired in breath-hold mode at the end of exhalation followed by segmentation into four tissue classes (air, lung, fat, and nonfat soft tissue) and assignment of constant attenuation values to each class. This static attenuation map was first registered to the MR reference gate image. The deformation fields were then applied to this image to generate attenuation maps for the remaining gates.

2.F. Simulation and patient studies
The PET/MR MCIR framework was evaluated using a simulation study and three patient studies. For each study, we compared the performance of the OG, UG, and MC cases. The study designs and figures of merits used for quantitative assessment are described below.

2.F.1. Simulation study
To validate the MAP MCIR framework for PET, we performed a simulation experiment using the 4D XCAT phantom, a spatiotemporal human torso phantom derived from CT datasets with cardiac and respiratory motion models obtained using 4D tagged MRI data and 4D high-resolution respiratory-gated CT data, respectively. 69 12 spherical pulmonary lesions were added, with six lesions each of diameter 14 mm in the right lung and liver and six of diameter 10 mm in the left lung and upper left abdomen. Standardized uptake values (SUVs) were assigned to emulate typical [ 18 F]FDG distributions: lungs 0.5, muscle 1, bones 2, liver 2.5, myocardium 4, and lesions 8 (all in units of g/ml). The resultant phantom is shown in Fig. 6. We simulated a 5 s breathing cycle divided into eight gates. Poisson deviates of the gated images were projected using the system model in (1) to generate 100 noisy sinograms mimicking realistic clinical count rates. Since the goal of the simulation study is to test the MAP reconstruction framework, we assumed an accurate system model. Intragate motion was therefore ignored in this setup.
For smoothing regularizers such as the uniform quadratic penalty, regularization-imposed bias leads to loss of spatial resolution. In setting the MAP regularization parameter, we seek to ensure fair comparison of the OG, UG, and MC cases (both in terms of image quality and quantitation) by ensuring similar regularization-imposed bias and hence similar image resolution in corresponding voxels. This is achieved by setting consistent relative weights between the data-fit and regularization terms for the three cases. Any differences in the observed F. 6. Coronal slice from the XCAT phantom showing the assigned SUVs in g/ml and simulated lesions labeled 1-12 in the thorax and upper abdomen.
bias (e.g., bias due to motion artifacts) would be attributable to inaccuracies in the data model alone. For MC, we set the regularization parameter β MC to 10 3 , 3×10 3 , 10 4 , 3×10 4 , 10 5 , 3×10 5 , and 10 6 . For fair comparison of OG and UG with MC, we select the regularization parameters β OG and β UG to ensure similar relative weighting between the data-fit term and the prior relative to MC. The log likelihood for MC in (5) is based on counts from all the gates, while the prior corresponds to only the reference gate. In comparison, for OG, the counts in the log likelihood in (3) are reduced roughly by a factor of p, the number of gates, while the prior remains unchanged since the reconstructed image x OG has roughly the same number of counts as x MC . The regularization parameter for this case was set to β OG = β MC /p. For UG, the log likelihood is based on the same number of detected events as MC and is therefore comparable. But the counts in the unknown image x UG are greater by a factor of p. The regularization parameter for this case was also set to β UG = β MC /p.

2.F.2. Patient studies
To validate the overall data acquisition and processing framework, we conducted PET/MR clinical experiments on three patients. The subjects scanned were a 24 yr old male (P 1 ), a 65 yr old female (P 2 ), and a 62 yr old female (P 3 ). PET list-mode data and GRANGE MR k-space data were simultaneously acquired on the Biograph mMR scanner. All procedures were performed following protocol and institutional standards with prior approval from the institutional review boards at Massachusetts General Hospital, Boston, MA, and SDN, Naples, Italy, the two data acquisition sites. The MR acquisition parameters in the studies were as follows: TR 3.3 ms/slice, bandwidth 1 kHz/pixel, readout oversampling factor of 2, and 256 samples per radial line. For P 1 , 24 coronal slices with a slice thickness 5 mm were used to cover the lung volume, and 4096 radial lines were acquired per slice. For P 2 and P 3 , 24 and 26 coronal slices, respectively, with a slice thickness of 8 mm were used to cover the lung volume with 4000 radial lines per slice. The excitation angle was set to 30 • based on the Ernst angle to provide higher SNR. A short-duration (∼3 ms) slice-projection navigator was inserted for every TR. The patients were injected with between 5 and 12.5 mCi [ 18 F]FDG radiotracer. For P 1 , P 2 , and P 3 , PET listmode data were acquired for 6, 7, and 7 min, respectively, simultaneously with the MR acquisition. The PET and MR data were processed following the steps described in Fig. 1. Gates were assigned by means of histogram equalization on correlation coefficients derived from the navigator projections. The PET list-mode events and MR k-lines were grouped into six amplitude-based gates per respiratory cycle. As our gating strategy does not distinguish between "wax" and "wane" phases of the respiratory cycle, these six gates therefore effectively represent 12 phases of motion per cycle. The binned kspace data were used to reconstruct gated MR images, from which the deformation fields were computed. The geometric component of the system matrix, P geom , for the Biograph mMR was computed by incorporating the block structure of the scanner and the gaps between the detector blocks. The point spread function of the scanner (in P blur ) was modeled using Monte Carlo simulations that accounted for blurring phenomena such as photon pair nonlinearity, intercrystal scatter, and crystal penetration as well as the block structure. 70 The gated list-mode data were binned into span 11 sinograms of size 344 × 252 × 837. The gated sinograms and the deformation fields were used to perform MAP MCIR. The reconstructed image dimensions were 256 × 256 × 129 and the voxel size was 2.025 × 2.025 × 2.025 mm 3 . A set of 30 MAP iterations was used in each case.

2.F.3. Evaluation metrics
For the simulation study, in which the ground truth is known, we rely on the bias-variance trade-off for validating the MAP reconstruction framework and for testing the proposed strategy for regularization parameter tuning. Accordingly, we plot absolute bias vs standard deviation curves for the overall volume. 71 As a figure of merit for both the simulation and clinical studies, we compute the contrast-to-noise ratio (CNR) 72,73 of high-intensity features in the lungs as Here,s and σ represent the mean and standard deviation, respectively, of intensities in the region of interest (ROI) and the background muscle tissue (a region with relatively uniform tracer uptake).
We also compute the mutual information (MI) between the features of interest in the UG, OG, and MC PET images with the corresponding feature in the MR image for the reference gate. This metric is a measure of the similarity between the intensity distributions of the MR and PET images. Both noiseinduced variations (as anticipated in the OG case) and motioninduced blurring (as anticipated in the UG case) tend to reduce this metric.
Overestimation of the volume of high-intensity features such as lesions is a well-known consequence of motioninduced blurring. We therefore compute the volume in ml for each feature of interest as another figure of merit for the clinical images. The volumetric mask is determined by using 50% of the peak intensity of the feature as a threshold.

3.A. Simulation results
The absolute bias and standard deviation images based on 100 noise realizations for each of the three cases with β MC = 3 × 10 4 , 10 5 , and 3 × 10 5 are shown in Fig. 7. The spatial F. 7. (A) Absolute bias and (B) standard deviation images for three different values of the MC regularization parameter, β MC = 3 × 10 4 , 10 5 , and 3 × 10 5 . For the OG case, we set β OG = β MC /p and for the UG case, we set β UG = β MC /p to ensure similar relative weighting of the data-fit and regularization terms and hence similar levels of regularization-induced bias. For each value of β MC , the OG and MC cases have similar bias and the UG and MC cases have similar standard deviation. The OG standard deviation is much higher due to fewer emission events. The UG bias is high due to the inaccuracy in the data model, which does not account for motion effects. averages of the bias and standard deviation for the three cases are tabulated in Table I. For each value of β MC , OG and MC show comparable bias (0%-10% difference) while OG exhibits increased standard deviation (658%-726% difference). The difference in standard deviation is as expected since the OG result is based on 1/8th the total counts. The UG and MC results exhibit similar variance behavior for each value of β MC (2%-5% difference). But the UG result shows a much stronger bias effect (16%-71% difference) manifesting as smearing artifacts caused by respiratory motion. This increase in bias is the result of the systematic error arising from not modeling the motion effects. The plot of the overall standard deviation against the overall absolute bias (computed as spatial averages) for all values of β MC is shown in Fig. 8(A).
We treated the 12 lesions as the ROIs and muscle tissue as the background. The CNR comparison for β MC = 3 × 10 4 is shown in Fig. 8(B). For the pulmonary lesions 1 through 4 and 7 through 10, the CNR improvement for MC relative to OG was in the range of 10%-15% without any clear pattern of dependence on location (and hence degree of deformation). In comparison, relative to UG, MC led to a consistent increase in CNR gain from the apex toward the base of the lungs. For the lesion pairs (1,7), (2,10), and (3,11), which undergo similar degrees of displacement and have similar background activity levels, the CNR gains for MC relative to UG were (21%,39%), (69%,100%), and (107%,120%), respectively. Thus the smaller lesions (7, 10, and 11) showed a larger margin of CNR gain than the larger lesions (1, 2, and 3). For the abdominal lesions (and also for lesion 4, located near the dome of the liver), the CNR values were influenced by the spillin from the liver and abdominal tissue, which were assigned higher activity levels than the lungs.

3.B. Experimental results
For MAP reconstruction of the patient data, the MC regularization parameter was set to β MC = 6. Since the data were binned into six gates, the regularization parameters for both the OG and UG cases were set to β OG = 1 and β UG = 1, respectively. Coronal slices from the MR and PET images of patient P 1 are shown in Fig. 9. Figure 9(A) shows the GRANGE MR image corresponding to the reference gate while Fig. 9(B) shows an overlay of the deformation field mapping the reference gate (end-exhalation) to the end-inhalation gate on the MR image. Figs. 9(C)-9(E), show the OG, UG, and MC PET images, respectively. As expected, the OG image is noisier than the UG and MC images, because it uses only about 1/6th the total photon count. Magnification of a blood vessel located near the diaphragm, marked using a green box in Fig. 9(F), indicates motion-induced blurring in the UG image, an artifact that is significantly reduced in the MC image. Profiles along the blue lines in Fig. 9(B) marked a, b, and c are shown in Figs. 10(A)-10(C), respectively. The profile in Fig. 10(A) corresponds to a region with less movement. We observe good matching between the UG and MC profiles, while the OG profile shows sharper noise-induced variations. Figure  10(B) shows a profile across the blood vessel that exhibits clear motion effects. We observe roughly a 34% improvement in peak recovered activity in the MC image relative to UG. Figure 10(C) shows a profile across the dome of the liver, another key region affected by respiratory motion. Compared to the MC image profile, the UG profile here shows consistent underestimation of activity in all voxels, while the OG result shows more pronounced noise effects. Coronal slices from the MR and PET images of patients P 2 and P 3 are shown in Figs. 11 and 12, respectively. High-intensity features near the diaphragm (including a lung lesion in the case of P 2 ) exhibiting distinctive motion blur patterns are magnified in Figs. 11(F) and 12(F).
We computed the volume, MI, and CNR for the features in Figs. 9(F), 11(F), and 12(F) to compare the OG, UG, and MC cases. These figures of merit are shown in Table II. Relative to OG, UG leads to an overestimation of feature volume by 11%-38%. In comparison, the overestimation for MC is 2%-19%. The overestimation of feature volume in the MC case is mainly attributable to inaccuracies in deformation fields, intragate motion, and interpolation errors. The OG and UG cases exhibit consistently lower MI with the reference gate MR image feature compared to the MC image. Although the OG image is morphologically the most similar to the reference gate MR image, it exhibits more fluctuations in intensity values due to noise effects, which very likely lead to a lower MI than the MC case. In the UG case, the lower MI can be attributed by the motion artifacts which impact the intensity distribution. For CNR computation, a uniform activity area of muscle tissue near the shoulder was treated as the background. Our results indicate that MC leads to a CNR improvement of 19%-196% relative to OG and 6%-51% relative to UG.  Fig. 9(B) showing the intensity distribution across a moving blood vessel. MC leads to a 34% improvement in peak recovered activity in this high-intensity feature relative to UG. (C) Profile along blue line c in Fig. 9(B) across the dome of the liver. Profiles (B) and (C) are from regions close to the diaphragm undergoing noticeable movement caused by respiration.

DISCUSSION
We have presented a complete data acquisition and processing framework for pulmonary PET imaging. Standalone MR is not the traditionally preferred choice for pulmonary applications. However, given the unique capabilities of simultaneous PET/MR in correcting for voluntary or involuntary patient motion, MR-based motion correction is advantageous for pulmonary PET imaging, which is greatly affected by breathing motion. We presented GRANGE, a dedicated golden-angle radial gradient echo MR pulse sequence that can be used in conjunction with k-t FOCUSS reconstruction to generate high-quality 4D MR images under free-breathing condition. While 3D encoding is also implementable, 74 the GRANGE sequence currently uses 2D encoding. Although we have not used parallel imaging in this work, the GRANGE technique can also be combined with parallel imaging.
Our PET reconstruction approach is based on MAP estimation with gradient ascent type optimization. Unlike unregularized reconstruction methods based on expectation maximization, our approach guarantees convergence. Additionally, it is feasible to derive closed form estimates of the resolution and variance of the MAP estimate for further statistical analysis. [75][76][77][78] We have also proposed and validated a strategy for regularization parameter selection based on consistent relative weighting of the data-fit and regularizer terms for the OG, UG, and MC cases. In the absence of consistent weighting, corresponding voxels could have arbitrary resolution and variance properties, making both quantitative and qualitative comparison of the three cases less meaningful.
We presented a retrospective gating strategy that is guided by a slice-projection navigator encapsulated in the GRANGE sequence. This method ensures roughly equal distribution of PET and MR data across the gates. Patient breathing may  vary substantially in course of an imaging scan. Our gate assignment technique does not rely on perfect periodicity of the respiratory signal and is therefore well-adapted to handle intercycle variability in breathing. While this gating method is capable of distinguishing between the "wax" and "wane" phases of each breathing cycle, we avoided it in this work for simplicity. Separating these phases in the gate assignment for future studies would ensure more robust performance in the presence of intracycle variability. In our studies, we used the end-exhalation gate as the reference gate. Previous studies have shown that, in a typical breathing cycle, patients spend the longest amount of time in an end-expiration quiescent phase. 36,79 Additionally, the scanner generated attenuation maps typically correspond to this phase. While our framework allows any gate to be used as the reference, the choice of end-exhalation as the reference reduces the overall susceptibility of the reconstruction result to errors in deformation field estimation and interpolation. During patient data acquisition, a spinal coil and a flexible surface coil were used on the mMR scanner. In our current reconstruction framework, the attenuation due to the spinal coil and the scanner bed is taken into account, while the attenuation of the surface coil is ignored. The errors caused by this coil are expected to be small as the coil is specially designed for this application. However, to reduce these errors and further improve accuracy, the attenuation of the flexible coil can be taken into account in the reconstruction by using markers to measure its location and shape. 80 In the clinical experiment, the PET and the MR data were simultaneously acquired for 5.5 min. For future studies, it may be feasible to exploit the robustness of k-t FOCUSS reconstruction to further reduce the MR acquisition time to <3 min. The saved time can, for instance, be used to run different diagnostic MR sequences while the PET data are collected in the background. As long as the slice-projection navigator is embedded in the other MR sequences, it would be feasible to use the simultaneously acquired PET events for MCIR. Compared to higher dimensional navigator signals, 81 the 1D slice-projection navigator is less accurate but faster. The latter characteristic makes it easily embeddable in MR sequences and hence a more practical choice for this application.