title: Motion correction for PET using subspace-based real-time MR imaging in simultaneous PET/MR
authors: Thibault Marin Yanis Djebra Paul K. Han Yanis Chemli Isabelle Bloch Georges El Fakhri Jinsong Ouyang Yoann Petibon Chao Ma
text: 
  
  PET motion correction using real-time MR imaging]Motion correction for PET using subspace-based real-time MR imaging in simultaneous PET/MR
  
  
  
  * Equal contribution
  † Corresponding author
  1 Gordon Center for Medical Imaging, Department of Radiology, Massachusetts General Hospital, Boston MA, 02114, USA
  2 Harvard Medical School, Boston MA, 02115, USA
  3 LTCI, Télécom Paris, Institut Polytechnique de Paris, France
  
  
    Image quality of PET reconstructions is degraded by subject motion occurring
    during the acquisition. MR-based motion correction approaches have been
    studied for PET/MR scanners and have been successful at capturing regular
    motion patterns, when used in conjunction with surrogate signals
    (e.g. navigators) to detect motion. However, handling irregular respiratory
    motion and bulk motion remains challenging. In this work, we propose an
    MR-based motion correction method relying on subspace-based real-time MR
    imaging to estimate motion fields used to correct PET reconstructions. We
    take advantage of the low-rank characteristics of dynamic MR images to
    reconstruct high-resolution MR images at high frame rates from highly
    undersampled k-space data. Reconstructed dynamic MR images are used to
    determine motion phases for PET reconstruction and estimate phase-to-phase
    nonrigid motion fields able to capture complex motion patterns such as
    irregular respiratory and bulk motion. MR-derived binning and motion fields
    are used for PET reconstruction to generate motion-corrected PET images. The
    proposed method was evaluated on in vivo data with irregular motion patterns.
    MR reconstructions accurately captured motion, outperforming state-of-the-art
    dynamic MR reconstruction techniques. Evaluation of PET reconstructions
    demonstrated the benefits of the proposed method over standard methods in
    terms of motion artifact reduction. The proposed method can improve the image
    quality of motion-corrected PET reconstructions in clinical applications.
  
  
  
  
  Keywords: PET motion correction, subspace modeling, low-rank reconstruction, PET/MR.
  
  
  
  
  
  
  § INTRODUCTION
  
  
  Motion, including physiological motion (i.e., cardiac and respiratory motions)
  and involuntary bulk motion, is a major source of image quality degradation in
  Positron Emission Tomography (PET), which can result in spatial blurring
  artifacts and mismatch between emission and attenuation maps, altering
  quantification of tracer concentration and deteriorating the diagnostic value of
  PET images <cit.>. The conventional way to
  handle motion in PET is the gating method, which bins PET list-mode data to
  different cardiac and/or respiratory motion phases followed by reconstructions
  of images of each phase. However, gating results in increased noise levels due
  to the reduced number of events in each motion phase. To address the
  limitations of the gating method, many PET motion correction methods have been
  developed, which consist of two consecutive steps: motion field estimation and
  motion correction by either applying the estimated motion fields to the gated
  images or modeling it within motion-compensated PET image reconstruction
  <cit.>.
  
  PET motion correction methods can be divided into two major subcategories,
  depending on how the motion field is estimated: PET-based methods and Magnetic
  Resonance (MR)-based methods. In the PET-based motion correction methods, the
  measured emission data are first assigned to specific motion phases based on
  surrogate signals <cit.>, e.g., electrocardiogram (EKG), respiratory
  bellow, optical tracking, etc. <cit.>, or the
  PET-data themselves <cit.>, e.g., center of mass,
  time-of-flight information, frame-by-frame images, etc. Motion fields are then
  estimated by registering the reconstructed image of each phase to a reference
  phase <cit.>. However, the accuracy of the motion fields estimated
  by the PET-based methods is limited by low signal-to-noise ratio (SNR),
  especially in the case of dual gating, and the overall lack of anatomical
  structural information of PET images <cit.>.
  
  The increasing availability of hybrid PET/MR systems provides a unique
  opportunity for mitigating effects of motion in PET using MR-based motion
  correction. Because of its excellent soft-tissue contrast, high spatial
  resolution, and high SNR, MR provides more accurate estimation of motion fields
  than the PET-based methods. MR-based PET motion correction methods have been
  successfully applied to compensate respiratory and cardiac motion in various
  applications involving both static and dynamic PET imaging
  <cit.>.
  One major limitation of the MR-based motion correction methods is that the
  conventional noniterative MR imaging methods are unable to resolve cardiac or
  respiratory motion in real time due to the slow imaging speed. Binning-based MR
  imaging methods <cit.> are
  often used to address this issue, where MR k-space data are grouped into
  different motion phases based on surrogate signals (e.g., EKG), navigator
  signals, or k-space data alone, and images of each motion phase are then
  reconstructed for the estimation of motion fields. However, the binning-based
  MR imaging methods suffer from three noticeable limitations. First, they assume
  pseudo periodic motion, which does not hold well in the case of arrhythmia and
  irregular respiratory motion. Second, they rely on either surrogate signals or
  navigator signals acquired along a single direction to assign k-space data to
  specific motion phases, which cannot reliably capture involuntary bulk motion.
  Third, their performance is limited by the inherent trade-off between the number
  of motion phases (and thus the accuracy of motion field measurement) and data
  acquisition time.
  
  
  In this work, we propose a real-time MR imaging method for PET motion
  corrections in PET/MR. High resolution real-time MR imaging is achieved by a
  subspace-based imaging method, which takes advantage of a unique property of
  high-dimensional dynamic MR signals known as partial separability (PS)
  <cit.>. The PS-model takes advantages of the spatial-temporal
  correlations of dynamic MR images, significantly reduces the number of unknowns
  of the underlying spatiotemporal signal, and makes it possible to recover high
  resolution, high frame-rate dynamic MR images from highly undersampled k-space
  data <cit.>. For PET motion correction, the
  reconstructed real-time MR images are used to determine motion phases and
  estimate motion fields. PET list-mode data are binned into sinograms
  accordingly and ordered-subset expectation-maximization (OSEM) reconstruction
  <cit.> is performed integrating the estimated displacement in the
  system matrix for motion correction. We demonstrate the performance of the
  proposed method by carrying out in vivo 18F-FDG PET/MR imaging
  experiments using a 3T simultaneous PET/MR scanner.
  
  
  
  
  § METHODS
  
  
  
  
   §.§ PET/MR imaging experiment
  
  
  An 18F-FDG PET/MR scan was performed on one healthy subject
  under a study protocol approved by our local IRB. PET and MR data were
  simultaneously acquired 30 minutes after 18F-FDG injection
  (around 10 mCi) using a 3T PET/MR scanner (Siemens Biograph mMR, Siemens
  Healthcare, Erlangen, Germany).
  
  Two 5-minute MR acquisitions were performed using a spoiled gradient-recalled
  echo (GRE) sequence with stack-of-stars radial sampling trajectories in the
  coronal plane. The imaging parameters are as follows: image size = 384
  × 384 × 32, resolution = 1.9 × 1.9 ×
  5 mm^3, TR/TE = 3/1.6 ms, and flip angle = 7 degrees.
  The (k,t)-space data were acquired using a random sampling pattern shown in
  figure <ref>. A total of 35 k-space spokes were sampled in
  each frame, resulting in a frame rate of 9.5 volumes per second. For each
  frame, the first 3 spokes were respectively acquired along the k_x, k_y
  and k_z direction across the center of the k-space to estimate the temporal
  basis of the partially separable (PS) model detailed next. The remaining 32
  spokes were along a random angle in the k_x-k_y plane for every k_z.
  During the first 5-minute acquisition, the subject was instructed to move once
  to assess the effect of both respiratory and bulk motion. During the second
  5-minute acquisition, the subject was instructed to simulate an irregular
  respiratory pattern including both deep and shallow breaths. The
  vendor-provided two-point Dixon sequence was performed with breath-holding to
  obtain attenuation coefficients.
  
  
  
  
  
   §.§ Subspace-based image reconstruction
  
  
  Denote the dynamic image series as ρ(,̊ t) and its matrix representation
  ∈C^N× M such that:
  
      = [[ ρ(_̊1, t_1)           … ρ(_̊1, t_M);           ⋮           ⋱           ⋮; ρ(_̊N, t_1)           … ρ(_̊N, t_M) ]].
  
  We express ρ(,̊ t) as a PS model <cit.>:
  
      ρ(,̊ t) = ∑_l=1^L u_l()̊ v_l(t),
  
  or, equivalently C as:
  
      =  ,
  
  where ∈C^L× M concentrates in rows the temporal
  basis function v_l of the PS model up to order L and ∈C^N× L concentrates in columns the corresponding spatial
  coefficients u_l.
  
  We estimate the temporal basis functions using the training data acquired at
  every frame. Assuming p training lines, we form the so-called Casorati
  matrix _t ∈C^p N_f × M by stacking the signal of the
  p training lines at each frame, where N_f is the number of samples
  acquired in each k-space line. The temporal basis functions
  {v_l}_l=1^L can then be estimated by calculating the first L right
  eigenvectors of _t using Singular Value Decomposition (SVD).
  
  Once  is obtained (denoted by ), the image reconstruction
  problem is reduced to the determination of the spatial coefficients matrix
  . We solve this problem by fitting the PS model to the undersampled
  (k,t)-space data with additional sparsity constraints <cit.>:
  
      = min_-̣Ω()_2^2 +
        λ_1 T()_1 + λ_2 _F,
  
  where d is the measured k-space data,  is the Fourier transform
  operator in the spatial domain, i.e. Non-Uniform FFT (NUFFT) operator
  <cit.> for the stack-of-stars trajectory, Ω is the sparse
  sampling operator in the (k,t)-space, T is the finite difference operator
  along both the spatial and temporal directions, ._F is the
  Frobenius norm and the scalar variables λ_1 and λ_2 are
  regularization parameters. The first term in Eq. (<ref>) is a data
  fidelity term, the second term promotes sparsity in the reconstructed image and
  the third term favors minimal norm solutions for .
  
  We solve the optimization problem in Eq. (<ref>) using the
  Alternating Direction Methods of Multipliers (ADMM) algorithm <cit.>,
  which leads to solving the following three sub-optimization problem in an
  alternative fashion:
  
      ^(k + 1) = 𝒮_λ_1/μ(
        T(^(k) ) + η^(k)),
        
      ^(k + 1) = min_1/2d - Ω(
        )_2^2 +
        μ/2T() - ^(k + 1) +
        η^(k)_F^2
       +  λ_2 _F,
        
      η^(k + 1) = η^(k) + (T(^(k + 1) ) -
        ^(k + 1)),
  
  where  is the split variable, η is the dual variable, and μ
  is a scalar relaxation parameter. The  update
  (Eq. (<ref>)) is a soft thresholding operation and the 
  update (Eq. (<ref>)) is a convex-optimization problem, which
  is solved using the conjugate gradient algorithm.
  
  For comparison, we reconstructed MR images using the same data by a binning-based
  method, known as XD-GRASP <cit.>. The respiratory motion signal used
  for binning was processed the same way  did: the Fourier
  transform of the data at the center of the k-space at each frame was sorted into
  a 2D matrix, with data from each coil concatenated along the first dimension. A
  Principal Component Analysis (PCA) was then applied on this matrix and the
  component with the highest peak in the respiratory frequency range (0.1, 0.5 Hz)
  was selected as the binning signal. The k-space data were then regrouped into 6
  balanced bins, i.e., each bin containing the same number of spokes. The
  XD-GRASP reconstruction was performed the same way as in <cit.>.
  
  
  
   §.§ Motion estimation
  
  
  The reconstructed real-time MR images were first binned into a small number of
  phases corresponding to different body positions (respiratory and bulk motion
  phases). Binning was performed in three steps. The first step consists in
  visually determining the bulk motion phases from the MR images and discarding
  time frames corresponding to the transition between bulk motion phases. In the
  second step, a bin is assigned to each (real-time) frame by tracking the tip of
  the right lobe of the liver over time while ensuring balanced bins (i.e., all
  bins should contain a similar number of frames). Finally, a combined MR image
  is formed for each bin by averaging all real-time images in a bin. Volumetric
  image registration was then performed between all bins and a reference bin using
  the multiscale B-spline registration algorithm described in <cit.>.
  
  
  
   §.§ PET reconstruction
  
  
  The acquired list-mode PET events were first rearranged into B sinograms =
  (_1, …, _B) following the binning determined from MR images and
  discarding PET list-mode events occurring during bulk motion transitions.
  
  PET reconstruction was performed using the OSEM algorithm <cit.>
  integrating the estimated motion fields in the forward model
  <cit.>. Let  denote the PET image to reconstruct
  arranged in vector form. The system matrix, denoted by  is decomposed as
  =, where:
  
  
    * = [_1, …, _B]^⊤ is a stack of deformation
  operators estimated using the procedure described in
  Section <ref>,
  
    * = diag(, …, ) is a block-diagonal geometrical
  projection matrix constructed by repeating the static projection matrix 
  implemented using Siddon's algorithm <cit.>,
  
    * = diag(_1, …, _B) is a diagonal matrix with
  time-varying attenuation coefficients,
  
    * = diag(, …, ) is a diagonal matrix with detector
  sensitivity coefficients  repeated for all bins.
  
  
  With these notations, the motion-corrected OSEM update for a given subset l is
  given by:
  
      ^(n + 1) = ^(n)/_l^⊤1_l^⊤_s/_l ^(n)  + _l,
  
  where _l is the system matrix for the l-th subset and _l is the
  combined additive correction sinogram for subset l including randoms and
  scatter. Correction sinograms were constructed as follows. Random coincidences
  were estimated using the delayed window method. Scatter was estimated using the
  single scatter simulation algorithm <cit.> from an initial
  reconstruction performed without motion correction. Scatter was estimated
  separately for each bulk motion phase. Attenuation coefficients were obtained
  from a vendor-provided Dixon sequence during breath-holding. The attenuation
  map was deformed to each bin and forward projected to calculate sinogram-domain
  attenuation coefficients.
  
  The OSEM used 12 subsets and 5 iterations. This motion-corrected reconstruction
  is denoted by MC in the rest of the paper. For comparison, two other
  reconstruction methods were considered: a traditional OSEM without motion
  correction (NMC) and a gated reconstruction where only list-mode events
  occurring in a given motion phase are reconstructed without motion correction
  (Gated). Both reference methods used 4 iterations to account for the difference
  in convergence speed, aiming to match the noise level in MC and NMC
  reconstructions.
  
  
  
   §.§ Quantitative analysis
  
  
  In order to compare PET reconstructions, two evaluation measures were used: the
  contrast-to-noise ratio (CNR) and target-to-background ratio (TBR). The
  contrast-to-noise ratio is defined as:
  
      CNR(x, ℛ) = x̅_ℛ -
                                           x̅_ℛ_0/σ_0,
  
  where σ_0 is the standard deviation in the background region
  ℛ_0 (located in the lung) and x̅_ℛ is the average
  activity of image x in region ℛ. With the same notations, the
  target-to-background ratio is given by:
  
      TBR(x, ℛ) = x̅_ℛ/x̅_ℛ_2,
  
  where ℛ_2 is a region located in the liver. Both metrics were
  evaluated in a small region located in the kidney (shown in
  figure <ref>(a)). In the absence of ground truth, the
  contrast-to-noise and target-to-background ratios were used as indicators of
  image quality.
  
  
  
  
  
  
  § RESULTS
  
  
  
  
   §.§ Correction of bulk motion
  
  
  In this experiment the subject was instructed to move after around 2.5 minutes
  in the 5-minute acquisition. Images reconstructed by XD-GRASP and the proposed
  method at the end inhalation and end-exhalation phases are shown in
  figure <ref>. The images obtained by XD-GRASP method show
  noticeable blurring artifacts largely because the bulk motion was not detected
  from the navigator signal. More specifically, figure <ref>(a)
  shows the navigator signal obtained from the training line along the k_z
  direction in each frame as in <cit.>. Since the bulk motion of the
  subject was along the x direction (left to right), the navigator signal only
  recorded abnormal changes during the bulk motion, i.e., the red region
  in figure <ref>(a), but did not contain sufficient information
  to indicate what type of motion occurred. Therefore, six motion bins were chosen
  in XD-GRASP while the k-space data acquired in the red region of
  figure <ref>(a) were discarded, resulting in blurring artifacts.
  
  
  
  The images obtained by the proposed method shown in
  figure <ref> successfully capture both respiratory motion (as
  indicated by the red and green dashed lines) and bulk motion (as indicated by
  the yellow arrows). Figure <ref>(b) to (d) show the temporal
  basis functions of the PS model estimated from the three training lines, where,
  intuitively, the first component (figure <ref>(b)) shows
  respiratory patterns and the second component (figure <ref>(c))
  indicates bulk motion. To further demonstrate the real-time capabilities of the
  proposed method, figure <ref> shows images at multiple time
  frames along with a 1D profile through the liver along time. The images from
  before and after bulk motion demonstrate the ability to capture both respiratory
  and bulk motion. The yellow overlay emphasizes the body displacement between
  Stage 1 and Stage 2. The profile plot shows the respiratory motion, captured
  for both bulk motion phases. The transition portion between the two bulk motion
  phases corresponds to the frames that were discarded in the PET reconstruction.
  Videos showing the reconstructed MR images are available in the Supplementary
  Material M1.
  
  
  
  
  
  To account for the two body positions and for respiratory motion in PET
  reconstruction, real-time MR images obtained by the proposed method were grouped
  into 12 bins (6 bins for each body position) for motion field estimation.
  Frames in the transition between the two bulk motion phases were excluded (a
  total of 15 seconds were discarded). Motion was estimated between all bins and
  the bin corresponding to the end-exhalation, which was used as reference bin.
  Estimated motion fields are shown in figure <ref>. The top
  left image shows motion caused by respiration, mostly visible as a vertical
  displacement near the liver (shown with more details in the inset image). The
  left column images show the bulk motion, which is mostly lateral. Finally, the
  bottom right figure shows a combination of respiratory motion between end
  inhalation and end exhalation and bulk motion.
  
  
  
  
  Reconstructed PET images are shown in figure <ref>.
  Figures <ref>(a) and <ref>(b) show coronal
  and axial slices using different reconstruction methods. Motion-blur is clearly
  visible on the NMC reconstructions, primarily in the lateral direction,
  corresponding to bulk motion but also in the vertical direction due to
  respiratory motion. The gated reconstruction, which uses one sixth of the PET
  counts at a single body position, shows sharper features but is severely
  corrupted by noise. The proposed method compensates both respiratory and bulk
  motion, significantly reducing motion-blur, while exhibiting a low noise level.
  Figure <ref>(c) shows line profiles through the kidney.
  Without motion correction (NMC), the activity peak is lowered by motion blur.
  Instead, two distinct peaks are visible, which correspond to the two bulk motion
  phases. Gated and motion-corrected reconstructions both preserve the peak
  activity, but gated reconstructions exhibit a high level of noise due to the
  reduced amount of data used for reconstruction.
  
  
  
  Evaluation measures are reported in table <ref>. The
  table shows that motion corrected reconstruction leads to the highest CNR: a 83%
  improvement was observed over reconstruction without motion correction and 198%
  over gated reconstruction. The noise level in NMC and MC is similar (within
  15%) but the contrast is substantially improved by motion correction, while the
  noise level in the gated reconstruction is three times higher leading to the low
  CNR. For the TBR, the gated reconstruction achieves the highest ratio, because
  gated reconstructions favor high contrast (at the expense of high noise). The
  proposed motion compensation method approaches the gated TBR (25% decrease) and
  outperforms NMC (20% increase).
  
  
  
  
  
  
   §.§ Correction of irregular respiratory motion
  
  
  The second experiment was designed to evaluate the performance of the proposed
  method in the case of irregular respiratory motion. The subject was instructed
  to alternate between slow deep and fast shallow breaths throughout the 5 minutes
  PET/MR acquisition.
  
  MR images obtained by the proposed method are shown in
  figure <ref>. The top row shows images at different frames:
  two at the end of inhalation and two at the end of exhalation taken from
  different breathing patterns (deep/shallow), respectively. The full extent of
  the respiratory motion is captured and the images are artifacts-free.
  Figures <ref>(e) and <ref>(f) show 1D
  profiles of the image through the liver changing over time. Both the images and
  the plot in figure <ref> clearly show the breathing patterns,
  alternating between deep slow breaths and fast shallow ones. Based on the
  reconstructed real-time MR images, 12 bins were determined through analysis of
  the liver displacement in the MR images and were consequently used for motion
  field estimation and motion corrected PET reconstruction. Sequences of MR
  reconstructions are shown in Supplementary Material M2.
  
  
  
  Corresponding PET reconstructions are shown in figure <ref>.
  Images reconstructed without motion correction (NMC) exhibit blurring artifacts.
  This is particularly visible on the left kidney (see the green line on the gated
  coronal image) where the bright spot visible on other images is elongated in the
  vertical direction, due to the large amplitude of the respiratory motion. The
  gated reconstruction uses one sixth of the total number of counts and therefore
  is degraded by noise, despite resulting in a sharper image. The proposed motion
  correction method results in the best image quality, in terms of noise and
  resolution. Corresponding line profiles are plotted in
  figure <ref>(c). The NMC peak is elongated along the y-axis,
  due to the large extent of the mostly vertical respiratory motion. The Gated
  line profile is sharper near its peak but has a large noise level. The proposed
  MC method results in a good compromise between sharpness and low noise.
  Contrast-to-noise and target-to-background ratios (defined in
  Eq. (<ref>) and Eq. (<ref>) respectively) are reported in
  table <ref> (regions of interest are shown in
  figure <ref>(b)). Metrics show the superior performance of
  the proposed motion correction method. The improvement in CNR is around 163%
  over NMC and over 200% over gated reconstruction. The TBR for the proposed
  method is within 15% of the gated TBR and around 95% larger than NMC.
  
  
  
  
  
  
  
  
  § DISCUSSIONS
  
  
  We have demonstrated the performance of the proposed MR-based motion correction
  for PET in two challenging cases: bulk motion and irregular respiratory motion.
  The proposed subspace-based MR imaging method allows for reconstruction of
  high-resolution 3D volumes at a rate of 9.5 volume/s, which enables accurate
  motion field estimation even in the case of irregular motions. Another
  important benefit of the proposed approach is the ability to perform informed
  binning for PET motion correction, rather than relying on navigators or external
  markers which offer limited information on the subject motion. With full
  real-time volumetric MR images, detecting motion becomes straightforward, and
  the process of determining an appropriate number of bins is greatly simplified.
  
  The key assumption of the subspace-based imaging method is the low-rank property
  of dynamic MR signals. We performed a simulation study to investigate this
  property in the case of regular and irregular respiratory motion. Two phantoms
  (shown in the Supplemental Material M3 and M4) were generated
  using the XCAT software <cit.> to simulate regular and irregular
  respiratory motion. Respiratory and cardiac cycles were divided into
  respectively 30 and 40 phases and 3D volumes were computed for each respiratory
  and cardiac phase combination (i.e. 1200 volumes). A 4D (3D space + time)
  phantom was then built by selecting and concatenating frame by frame the 3D
  volumes based on simulated EKG and respiratory signals. Both the breathing
  frequency and diaphragm expansion were varied while keeping a constant heart
  rate in the simulation of the irregular respiratory motion. Each phantom
  contained 6 respiratory cycles. The contrast was designed to simulate a
  Balanced Steady-State Free Precession (bSSFP) signal for several compartments
  such as fat, muscles, etc. using T1 and T2 values from the literature
  <cit.>. SVD was then performed to investigate the effect of an
  irregular respiratory pattern on the rank, and its corresponding approximation
  error with low-rank truncation (see Supplemental Material figure M5).
  The decay of the calculated singular values from both phantoms was very similar,
  indicating that the breathing pattern does not substantially affect the rank of
  the data.
  
  The proposed method utilizes an MR acquisition which fully overlaps with the PET
  acquisition and provides real-time MR images for motion correction. The
  proposed method can still have benefits for other commonly used acquisition
  protocols. It is common in practice to reserve a first part of the PET
  acquisition to perform MR motion field measurements and use the remaining PET
  acquisition time to perform additional MR measurements (e.g. using T1 or T2
  contrast sequences) that can be used for other diagnostic tasks
  <cit.>. The proposed method can advantageously replace the motion
  field measurement sequence, possibly reducing the acquisition time while
  preserving image quality. A gating signal (e.g. navigator or external marker)
  can then be used in subsequent MR sequences to select an appropriate bin for
  each PET frame. Another approach is to integrate contrast sequences into the
  motion field estimation sequence described in this paper. This is under
  investigation and will be reported in separate publications.
  
  
  The study reported in this paper has several limitations. First, the computation
  time for the low-rank reconstruction with sparsely sampled non-Cartesian k-space
  data could be a concern. The current MATLAB (The MathWorks, Inc., Natick,
  Massachusetts, United States) implementation performs reconstruction of one
  slice and one coil in around one hour. We anticipate that using a lower level
  programming language and parallel computing devices (e.g. GPU) will help achieve
  reasonable runtimes <cit.>.
  Second, the proposed method does not have sufficient temporal resolution to
  resolve the motion in the transition phase between the two bulk motion phases of
  experiment 1 (figure <ref>(e)). The time-varying
  profile plot shows that the image quality in the transition is severely
  degraded. The corresponding list-mode data were excluded from the PET
  reconstruction. Since the duration of the bulk motion was short, only about 5%
  of the list-mode data were discarded and thus should not be a significant
  limitation. Third, this study focuses on demonstrating the feasibility of using
  subspaced-based real-time MR for PET motion correction. We showed the
  performance of our method in two cases (bulk motion and irregular respiratory
  motion) from in vivo PET/MR experiments on a healthy subject. More subjects are
  needed to fully evaluate the performance of the proposed method in clinical
  settings.
  
  
  
  
  § CONCLUSION
  
  
  We proposed an MR-based method for PET motion correction using a subspace-based
  real-time MR imaging for motion field estimation. We demonstrate the
  feasibility of the proposed method using 18F-FDG-PET/MR
  studies on a healthy subject. Our results show that the proposed method can
  capture and correct for normal and irregular respiratory motions as well as
  bulk body motion. The proposed method can be beneficial to a range of clinical
  applications where irregular motion patterns are expected.
  
  
  
  
  § ACKNOWLEDGMENTS
  
  
  This work was supported in part by the National Institutes of Health under award
  numbers: T32EB013180, R01CA165221, R01HL118261, R21MH121812, R01HL137230 and
  P41EB022544.
  
  
  jphysicsB
  
  
  
  