A 3D MR-acquisition scheme for nonrigid bulk motion correction in simultaneous PET-MR

Purpose: Positron emission tomography (PET) is a highly sensitive medical imaging technique commonly used to detect and assess tumor lesions. Magnetic resonance imaging (MRI) provides high resolution anatomical images with different contrasts and a range of additional information important for cancer diagnosis. Recently, simultaneous PET-MR systems have been released with the promise to provide complementary information from both modalities in a single examination. Due to long scan times, subject nonrigid bulk motion, i.e., changes of the patient’s position on the scanner table leading to nonrigid changes of the patient’s anatomy, during data acquisition can negatively impair image quality and tracer uptake quantiﬁcation. A 3D MR-acquisition scheme is proposed to detect and correct for nonrigid bulk motion in simultaneously acquired PET-MR data. Methods: A respiratory navigated three dimensional (3D) MR-acquisition with Radial Phase Encoding (RPE) is used to obtain T1- and T2-weighted data with an isotropic resolution of 1.5 mm. Healthy volunteers are asked to move the abdomen two to three times during data acquisition result-ing in overall 19 movements at arbitrary time points. The acquisition scheme is used to retrospectively reconstruct dynamic 3D MR images with different temporal resolutions. Nonrigid bulk motion is detected and corrected in this image data. A simultaneous PET acquisition is simulated and the effect of motion correction is assessed on image quality and standardized uptake values (SUV) for lesions with different diameters. Results: Six respiratory gated 3D data sets with T1- and T2-weighted contrast have been obtained in healthy volunteers. All bulk motion shifts have successfully been detected and motion ﬁelds describing the transformation between the different motion states could be obtained with an accuracy of 1.71 ± 0.29 mm. The PET simulation showed errors of up to 67% in measured SUV due to bulk motion which could be reduced to less than 10% with the proposed motion compensation approach. Conclusions: A MR acquisition scheme which yields both high resolution 3D anatomical data and highly accurate nonrigid motion information without an increase in scan time is presented. The proposed method leads to a strong improvement in both MR and PET image quality and ensures an accurate assessment of tracer uptake. © 2014 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License.


INTRODUCTION
Positron emission tomography (PET) is a commonly used tool to detect and quantify FDG uptake in tumors. 1,2 Recently, hybrid PET-magnetic resonance (MR) scanners have been introduced to combine the excellent sensitivity of PET with high-resolution multicontrast MR images. Sequential PET-MR systems consist of PET and MR components which are spatially separated to avoid any mutual interference allowing for a coregistered PET-MR acquisition. 3 For simultane-ous PET-MR systems, the PET and MR gantry are integrated such that PET and MR cover the same field of view. After synchronizing the clocks which control data acquisition of PET and MR, these systems yield simultaneously obtained image information. 4,5 The advantage of simultaneous PET-MR systems is the potential to speed up the acquisition of PET-MR data and providing additional information about motion from MRI.
In contrast to computed tomography (CT), MR imaging (MRI) does not just yield anatomical but also semifunctional information such as diffusion-weighted images. Such multimodality data optimally complements each other for a highly accurate detection and assessment of metastases. [6][7][8] PET data are usually acquired in multistations, where several 3D data sets are obtained at different bed positions. This acquisition mode is also applied in simultaneous PET-MR allowing for the acquisition of different 3D MR datasets at each bed position. These datasets can be images with different contrast as well as image information used for generation of attenuation correction (AC) maps. Advances in data acquisition and image reconstruction allow for a spatial resolution of 1 mm in MRI and 3-5 mm in PET within a clinically feasible scan time. 9,10 However, patient motion during 3D data acquisition in both imaging modalities can be a major challenge for PET-MR. 11,12 Apart from respiratory motion, which can be compensated for with respiratory gating techniques, bulk patient motion, i.e., a change of the patient's position on the scanner table leading to nonrigid changes of the patient's anatomy, during data acquisition is challenging to detect. The quality and therefore diagnostic information of high-resolution images can be severely impaired by wholebody shifts of only a few millimeters during data acquisition. In particular, the shift of the imaged object can result in severe artifacts in MR images and can cause blurring of lesions and a reduction in measured standardized uptake values (SUV) in PET. Furthermore, it also leads to a mismatch between AC map and acquired PET emission data which can affect SUV evaluations as well. 13 In contrast to PET-CT, where data are obtained sequentially, PET-MR allows for simultaneous acquisition of PET and MR data with different contrasts. [14][15][16][17][18][19] Therefore, MRI can be used for detection and compensation of motion. Several techniques have been proposed which use this idea for motion correction of the heart or abdomen. Nevertheless, the majority of these methods relies on the motion to be periodic and is therefore not applicable to bulk motion correction. Other approaches require dedicated MRI sequences which are optimized for motion estimation but do not yield high resolution anatomical data required for medical diagnosis and thus strongly restrict information available from MR.
To overcome these limitations we present the use of a 3D high-resolution Radial Phase Encoding (RPE) MR acquisition scheme, which allows for the reconstruction of images with different temporal resolutions from the same acquired data. 20,21 The aim of this study is to use this acquisition technique as an image-based motion scout to automatically detect and estimate nonrigid bulk motion from the acquired MR data. The obtained motion information is then used to correct for motion in both PET and MR images. In contrast to previous applications of the RPE sampling scheme, the presented technique is not limited to periodic motion such as respiratory motion or to affine motion compensation frequently used to correct for the movement of the heart during breathing.
The proposed approach was successfully assessed in six datasets obtained in three healthy subjects for T1-weighted (T1w) gradient (GRE) and T2-weighted (T2w) turbo spin (TSE) echo images. In addition, due to a not available PET-MR system, simulations were carried out to assess the effect of bulk motion and motion correction on PET image quality and lesion assessment.
FIG. 1. Overview of the proposed bulk motion correction approach. Bulk motion shifts during data acquisition are detected and corrected in the obtained MR data which strongly improves image quality. Additionally, the information on bulk motion shifts can also be used to improve the image quality of a simulated simultaneous PET acquisition. For more details, please refer to Sec. 2. Steps which are required for a bulk motion compensated simultaneous PET-MR acquisition are indicated as filled gray boxes. All other steps are only necessary for PET simulation and data evaluation. AC map: attenuation correction map; Em map: emission map; MCIR: motion corrected image reconstruction; and OSEM: ordered subsets expectation maximization. Bulk motion of the body during MR data acquisition can cause severe artifacts in the obtained MR image. The illustrations used here visualize a nonrigid bulk motion shift of the abdomen after 50% of the total scan time which leads to deformation of the spine and inner organs. (Motion detection) The acquired RPE data can be split up into multiple time frames (D i ) each covering a short time interval of the total scan time. Each D i is reconstructed using only part of the acquired data but covers the same field of view. This can lead to undersampling artifacts which are visible as streaking in the obtained images. Despite these artifacts a reliable affine image registration can be carried out. The obtained affine motion fields aMF i-j describe the affine transformation between time frame i and j. To detect bulk motion shifts, the average over the length of all motion vectors (i.e., the average displacement) in each motion field aMF i-j is calculated as R i (j). In the presented example, D 1 and D 2 are in the same motion state; therefore, the length of all motion vectors in aMR 1-2 is 0 and R 1 (2) is 0. D 1 and D 3 describe two different motion states; therefore, the length of the motion vectors in aMR 1-3 is larger than 0 and R 1 (3) is greater than 0. Plotting all functions R i (j) shows that these functions cross when nonrigid bulk motion shifts occur which allows for the detection of the time point TBM (R 2 and R 3 are not shown because they are identical to R 1 and R 4 , respectively). (Motion estimation) Based on the motion detection, the acquired data are split into multiple bulk motion states (in this case two motion states) and images B k describing each bulk motion state are reconstructed. A 3D nonrigid image registration algorithm is used to determine the motion between the individual bulk motion states. The obtained nonrigid motion fields MF k-l can then be applied to a motion compensated MR or PET reconstruction.

METHODS AND MATERIALS
An overview of the methodology is given in Fig. 1. Figure 2 depicts the proposed motion detection and estimation approach in more detail. A brief outline of the nonrigid bulk motion compensation method is given below and a thorough discussion follows in Secs. 2.A-2.E. MR data are obtained with a RPE scheme which allows for the reconstruction of 3D data sets with different temporal resolutions from the same acquired data. 20 To compensate for bulk motion, images with three different temporal resolutions are reconstructed. In a first step the recorded data are split into time frames D i with a high temporal resolution by combining data from a short interval of the total scan time. These so-called dynamics describe the entire FOV at different time points during data acquisition. The temporal resolution of D i can be adapted retrospectively and is in the range of several seconds compared to a total acquisition time of several minutes. Images are reconstructed from this dynamic data to detect the time point(s) T BM when bulk motion occurred. Since the image quality of the dynamic images is low, they are only used for motion detection and not for motion estimation. Second, images are reconstructed before and after the time point(s) T BM describing each bulk motion state BM. More data are used for this image reconstruction and therefore the image quality is high enough to estimate nonrigid motion between each state BM using a local affine registration approach. 22 Finally, the estimated motion fields are used to correct for bulk motion shifts in both MR and PET data.

2.A. MR data acquisition and reconstruction
Data acquisition are carried out using a bit-reversed RPE scheme (Fig. 3). This sampling scheme is similar to the scheme presented by Boubertakh et al. 20 for cardiac MRI. The RPE sequence is a modification of a standard 3D Cartesian MR acquisition approach. For a 3D Cartesian MR sampling scheme, frequency encoding is carried out along parallel lines. The individual frequency encoding lines are arranged on a Cartesian grid in the 2D phase encoding (PE) plane. RPE uses Cartesian frequency encoding but the PE points are located along radial lines. This sampling scheme has been shown to maintain good image quality even for high undersampling factors (i.e., few radial PE lines) if images are reconstructed with an iterative non-Cartesian SENSE reconstruction approach. 20 The radial PE lines are recorded in a bit-reversed order. This ensures that the angle between two successively recorded PE lines is large, i.e., PE lines are sampled at 0 • , 90 • , 45 • , 135 • , 22.5 • , 112.5 • , etc. The bit-reversed order leads to a homogeneous covering of the 2D PE plane over time and allows for the reconstruction of images from temporal subsets (dynamics) of the acquired data.
After the MR scan is finished, a high-resolution 3D image can be reconstructed using all acquired radial PE lines. In addition, dynamic images with different temporal resolution can be obtained from the same data.

2.B. Bulk motion detection
In order to detect the time points T BM when bulk motion has occurred during data acquisition, recorded data are split into N dynamics D i . Each D i is reconstructed using a small number of subsequently acquired radial PE lines and therefore describes a temporal snapshot much shorter than the total scan duration (Fig. 2).
It is important to note that each radial PE point represents a fully sampled frequency encoding line. Each dynamic covers the entire FOV with less data compared to the recorded data which leads to undersampling artifacts. Nevertheless, due to the undersampling properties of RPE combined with a non-Cartesian SENSE reconstruction scheme the image quality is sufficient for bulk motion detection.
For motion detection, each of the dynamics D i is registered to all other dynamics D j using an affine registration algorithm. 22 This registration yields motion fields aMF i-j describing affine transformations between D i and D j . For each motion field aMF i-j , the mean displacement R i (j) is calculated as the average length of all motion vectors of aMF i-j . If no motion has occurred between D i and D j , then the length of all motion vectors is zero and the mean displacement value is also zero. If on the other hand a bulk motion shift has occurred, the mean displacement value will be larger than zero. This approach is repeated for all dynamics which leads to N displacement functions R i . It is important to note that this technique does not require a prior definition of a threshold but all required information is obtained from the image registration. Figure 2 shows two displacement functions R 1 and R 4 for dynamic D 1 and D 4 which describe two different bulk-motion states. Each of them shows a sharp increase or decrease when bulk motion occurs which allows for the detection of the time points T BM as the crossings of the displacement functions.
We compared this approach to using normalized image correlation coefficients obtained between the individual dynamics for bulk motion detection.

2.C. Bulk motion estimation
Based on T BM , a second set of images B k is reconstructed describing each of the bulk motion states (Fig. 2). Depending on the motion, data acquired close to T BM can be excluded from image reconstruction to ensure the obtained images are not corrupted by bulk motion.
The images B k are registered using a 3D hierarchical nonrigid registration. 22 The obtained motion fields M k describe the transformation of each bulk motion state BM k to a reference motion state BM Ref . Here we selected the bulk motion state with the highest number of radial PE lines as BM Ref .

2.D. Bulk motion compensation of MR images
The final 3D motion compensated MR images are obtained by using a motion compensated non-Cartesian iterative SENSE reconstruction 23, 24 using the motion fields M k obtained in the previous step.

2.E. Bulk motion compensation of PET images
Similar to the compensation of the MR data, T BM and M k are used in a motion compensated PET reconstruction to compensate for bulk motion shifts. A motion corrected image reconstruction (MCIR) (Ref. 25) approach using the motion fields M k is applied to correct for subject motion. Nonrigid bulk motion shifts are compensated for by directly applying M k to the activity distribution in each voxel during a conventional iterative PET image reconstruction and transform the voxels to the same motion state. In contrast to approaches which carry out motion correction in the image domain (reconstruct-transform-average), MCIR has been shown to yield more accurate uptake values. 26 As discussed previously, AC data required by the PET reconstruction need to be adapted for each bulk motion state in addition to the correction of the PET emission data.

2.F. Volunteer experiments
A bit-reversed RPE sampling scheme was implemented on a 3T MRI scanner (Philips Healthcare, Best, The Netherlands) for both gradient echo (GRE) and turbo spin echo (TSE) sequences. In three volunteers, T1w GRE (flip angle: 7 • , repetition(TR)/echo time(TE): 6.0/3.0 ms, total scan time: 4 min) and T2w TSE (flip/inversion angle: 90 • /120 • , repetition(TR)/echo time(TE): 365/7.1 ms, turbo factor: 16, total scan time: 7 min) data were acquired using a 32 channel cardiac phased array coil: 288 mm 3 FOV, acquisition matrix: 192 × 192 × 128 (undersampling factor of 1.5 along the angular direction), and 1.5 mm 3 isotropic resolution. In order to reduce scan times, a partial Fourier factor of 0.75 was used along the radial PE direction for the T2w TSE acquisition. Respiratory navigation (navigator window: 6 mm) was employed in both sequences to minimize respiratory motion artifacts using a pencil beam MR navigator placed on the right hemi-diaphragm. During data acquisition, volunteers were advised to carry out several bulk motion shifts of the lower abdomen. Written informed consent was obtained from all participants and the study was carried out following an approved protocol from our local hospital.
All images were reconstructed offline using MATLAB (The MathWorks, Inc., Natich, MA) and the required coil sensitivity maps were calculated from a separate reference scan. In order to compensate for the partial Fourier acquisition, the T2w TSE data were reconstructed using a homodyne approach. 21 In addition to the free-breathing high-resolution scans, a dual-echo DIXON acquisition was carried out between the T1w GRE and the T2w TSE sequence to obtain emission and attenuation maps for the PET simulation: 27 repetition(TR)/echo times(TE): 3.3/1.2/2.1 ms, 288 × 240 × 288 mm 3 FOV, 1.5 × 3 × 3 mm 3 resolution. This scan was obtained in one breathhold of approximately 20 s.
The acquired MR data were split into 55 dynamic time frames D i (each with 20 radial PE lines) using a sliding window approach (window shift: two radial PE lines) each covering the entire FOV. In order to avoid artifacts from corrupted data during the actual bulk motion shift, one radial PE line before and after T BM were excluded from the image reconstruction. This means approximately 4 s (T1w GRE) and 7 s (T2w SE) of data acquisition were rejected around each T BM .

2.G. Data analysis
To quantify the bulk motion shifts eight anatomical landmark points, LM were manually selected in the 3D MR images describing each bulk motion state BM using distinct anatomical features. The difference of the landmark position between different BM provides a quantitative measure of the displacement in different regions of the FOV due to bulk motion shifts. The maximum displacement of these landmark points was determined for each scan and for each volunteer: where i and j describe different bulk motion states, . is the Euclidean norm and lm = 1. . . 8 are the indices of the different landmarks.
The accuracy of the obtained motion fields was assessed using the target registration error (TRE) between different pairs of landmarks. 28 For this study, the 3D TRE was determined by transforming the LMs to the reference motion state BM Ref and comparing them to the landmarks LM Ref manually selected in the reference motion state: where ms refers to the different N bulk motion states described by the motion fields M ms .

2.H. Accuracy of nonrigid motion detection and estimation
The accuracy achieved with the proposed motion compensation approach depends on the amplitude and duration of the bulk motion states. In order to assess the shortest motion state which can still be accurately detected and estimated with the proposed method, nonrigid bulk motion shifts were simulated using image data and motion fields acquired in a volunteer.
Small motion leads to little difference in the individual dynamic images D i and is more challenging to detect. Therefore, the volunteer scan with the smallest overall landmark displacement was selected for this study.
An MR data acquisition with three different bulk motion states was simulated with different durations of the second bulk motion state ranging from 25% to 3% of the total scan time. The duration of the third bulk motion state was adapted accordingly.
The error of the motion detection approach (E Det ) was determined as the difference between the detected time points T D BM and the ground truth time point T BM used in the simulation for different motion state durations: with p = [3%, 5%, 6%, 9%, 13%, 16%, 19%, 22%, 25%] describing the duration of the second bulk motion state relative to the total scan time. This assessment was carried out for two different temporal resolutions of the dynamic images.
To assess the motion estimation approach the same simulated MR data were used. Nonrigid motion fields M D were determined for different durations of the second bulk motion state. The error of these motion fields E Est was determined by calculating the mean target registration error (TRE) over eight landmarks LM similar to Eq. (2): The main difference to the TRE from above is that the landmarks LM* were not manually selected but calculated by applying the motion fields used for creating the simulated data to the reference landmarks LM Ref . This ensures that any errors are due to inaccuracies of the motion fields and not due to the selection of the landmarks.
In addition, the image reconstructed from the data of the second bulk motion state was also transformed to the reference motion state using M D . The cross-correlation value between the transformed and the reference motion state was also calculated to evaluate if it could be used as a surrogate of the accuracy of the motion estimation.

2.I. PET simulation
Realistic PET data acquisition was simulated using semiautomatically segmented MR images with manually inserted virtual lesions. 29 For this, the 3D DIXON MR images were separated into bone, liver, lung, kidneys, soft tissue, and background signal. The 3D DIXON MR scans yielded two 3D data sets acquired at two different echo times which allowed for an automatic segmentation of water, fat, background, and lung image information. Kidneys and liver were separated from other soft tissue manually. Bone information was segmented using an approach based on signal inversion as suggested by Helle et al. followed by a manually removal of regions which were falsely classified as bone due to their low signal such as air in the bowels. 30 Typically measured standardized uptake values (SUV) for FDG were assigned to the different tissue types to create 3D emission maps: air 0, lung tissue 0.5, soft tissue and fat 1, bone 0.3, liver 2.5, and kidneys 3.0. These maps were then transformed to the different bulk motion states using the motion fields obtained with the RPE acquisition scheme.
Soft tissue tumor lesions with different diameters varying between 7.5 and 19.5 mm were placed manually in the emission maps of each bulk motion states at accurately determined anatomical landmarks. The uptake values assigned to the lesions were 11 for lesions in the kidney and 8 for lesions close to the spine. AC maps were created by assigning the following attenuation values: 31 background 0 cm −1 , bone 0.15 cm −1 , fat 0.09 cm −1 , lung tissue 0.03 cm −1 , and soft tissue (also including liver and kidneys) 0.1 cm −1 . The same AC information was used for the PET simulation and PET image reconstruction to ensure that any errors of the measured uptake values are due to the effects of bulk-motion.
The image resolution of the AC and emission maps was reduced from the high MR resolution to a realistic PET resolution using 3D Gaussian filtering with a kernel width of 6 mm.
AC and emission maps were used in an analytic 4D PET simulation approach based on the open-source library STIR (Software for Tomographic Image Reconstruction, release 2.4) (Ref. 32) to obtain PET sinograms for each bulk motion state. This very fast approach generates realistic 4D PET-MR datasets from measured MR acquisitions. It provides simulations of simultaneously acquired data which are useful for the validation of motion correction strategies and has been previously compared against Geant4 based Monte Carlo simulations. 29 A Philips Gemini TF scanner was simulated with a cylindrical bore (90 cm diameter, 18 cm length) and 28 detector blocks with 44 × 23 crystals each. The crystals measured 22 × 4 × 4 mm 3 .
Compton scattering was simulated analytically in 3D such that it accounts for 35% of the obtained events. 33,34 The emission maps were scaled to a total number of 75 × 10 6 counts corresponding to a 5 min 3D FDG-PET respiratory-gated acquisition. To generate multiple realizations of the acquisition process, Poisson noise was added taking 20 different seeds from a random number generator for each realization.
Image reconstruction was performed with STIR using the iterative three-dimensional ordered subsets expectation maximization (OSEM) algorithm with 23 subsets and two full iterations and 4 mm isotropic 3D Gaussian post-filtering. Subject motion during PET data acquisition impairs the final image quality in two different ways. It leads to a blurring of lesions and to a misalignment between attenuation maps and reconstructed emission data. To study the impact of these two effects on the final image quality, the following PET images were reconstructed: r Motion free (MF) reference image: All PET data were acquired using the emission map of the reference bulk motion state. were obtained for each of the bulk motion states. A MCIR approach using the motion fields M k obtained from the MR data were applied to correct for subject motion. For MCIR, the motion information is used to transform the recorded PET data to a common reference motion state during PET image reconstruction which yields a final motion corrected PET image.
r Motion correction of sinograms and attenuation correction maps (MCSA): Identical to MCS but the attenuation correction maps were also transformed to the different bulk motion states before they were applied in the MCIR approach.

2.J. PET data analysis
The relative error (RE) between the motion free reference image (MF) and noMC, MCS, and MCSA was calculated as where C is the SUV of each image voxel. In order to assess the effect of motion correction on tumor visibility and uptake of the radiotracer, the maximum (SUV max ) uptake value was calculated within a spherical region-of-interest (ROI) centered at each lesion. 35 In addition, the contrast recovery coefficient (CRC) relative to the motion free case was calculated for all lesions: where SUV H and SUV B describe the average uptake value for lesions (H) and background (B), respectively. All parameters were assessed in MF, noMC, MCS, and MCSA images. All measurements, definition of spherical ROIs and analysis were carried out using MATLAB (The MathWorks, Inc., Natich, MA).

RESULTS
The MR data acquisition and PET simulations were successfully carried out for all volunteers. Every bulk motion shift in any of the six datasets could be accurately detected and corrected.
The individual steps of the bulk motion detection and estimation are depicted in Fig. 4 for one volunteer. Figure 4 The white outline describes the position of the kidneys at the beginning of the MR scan and is used to visualize the displacement between these three time frames describing the beginning, middle, and end of data acquisition. (c) An affine registration is used to determine the average displacement functions R of all dynamics. R of dynamic 10, 30, and 52 are shown here. A low value of R 10 between dynamic 1 and 20 indicates that these dynamics are in the same motion state as dynamic 10. The following sharp increase of R 10 suggests that bulk motion has occurred. The time points when bulk motion shifts took place are detected as crossings of the displacement functions (black arrows). Based on this information, the acquired data are split into multiple bulk motion states (three motion states in this case). (d) One slice of the 3D image volume for each of the bulk motion states B 1 , B 2 , and B 3 . Depending on the duration of each bulk motion state, more data are used for the image reconstruction of B 1 , B 2 , and B 3 compared to the dynamics D leading to a superior image quality. (e) Difference images between motion state B 1 , B 2 , and B 3 , respectively. Bulk motion leads to displacements between the different motion states which are clearly visible in the difference images. A local affine registration algorithm is used to obtain motion fields describing the transformation between the different bulk motion states (i.e., transformation between B 1 and B 2 and B 1 and B 3 , respectively). (f) Applying the obtained motion fields leads to the corrected data cB 1 , cB 2 , and cB 3 which are now all describing the same motion state. (g) The difference images of the corrected data do not show any visible displacement anymore. shows in vivo images with high temporal resolution used for motion detection. Each image is reconstructed using only 20 radial phase encoding lines. Despite this high undersampling factor, the image quality of the dynamic images was sufficient to obtain affine motion estimations which allowed for a reliable detection of the time points when bulk motion occurred [ Fig. 4(c)]. The images used for motion estimation are reconstructed from a variable number of radial phase encoding lines depending on the duration of each bulk motion state [ Fig. 4(d)]. For the volunteer shown here, these images were obtained using at least 34 radial phase encoding lines for each image leading to a better image quality compared to the dynamic images shown in Fig. 4(b). They allowed for accurate nonrigid motion estimation [Figs. 4(e)-4(g)]. Figure 5 compares obtained displacement (R) and normalized correlation coefficient (Corr) functions for two volunteers. Both R and Corr show a similar behavior with a high correlation coefficient corresponding to a low displacement. Using thresholds of 2 mm for R and 0.9 for Corr, leads to a reliable detection of the three bulk motion states BM 1-3 in volunteer 1. BM 1-3 can also be determined with the same threshold for R but not for Corr in volunteer 2. With the proposed approach of determining the crossings of R and Corr [Figs. 5(c) and 5(d)], BM 1-3 can be detected in both volunteers. The proposed approach is subject independent and does not require the manual selection of a threshold value. For the volunteer experiments, T BM was determined from the crossing of R rather than Corr because R provides quantitative information which does not just show the occurrence of bulk motion but also its amplitude in millimeters.
Motion corrupted and motion corrected MR images are depicted in Fig. 6. Bulk motion shifts during data acquisition led to strong artifacts and severely impaired MR image quality. The proposed approach could accurately detect and correct for bulk motion shifts which strongly improved the images.
The maximum displacement of the selected landmarks ranged from 10 mm to more than 30 mm. The target registration error was 1.71 ± 0.29 mm over all volunteers and scans. Figure 7 shows the displacement and the TRE for each volunteer and each scan. There was no visible dependency of the TRE on either the scan type (T1-or T2-weighted) or on the motion amplitude (i.e., the maximum displacement of the landmark points).
The accuracy of the proposed motion compensation approach for different durations of a bulk motion state is depicted in Fig. 8. For both temporal resolutions (55 dynamic and 115 dynamic time frames), the accuracy decreases with decreasing duration of the bulk motion state. For the 55 dynamic cases, the motion detection fails for motion states shorter than 9% of the total scan time (i.e., the time it takes to acquire 128 RPE lines). Splitting the data into 115 dynamic time frames increases the temporal resolution of the motion detection approach and bulk motion states with less than 5% of the acquired data can be detected. For the T1w scan, this means a bulk motion state with a duration of approximately 11 s can still be correctly identified.
The error of the motion estimation E Est is below the spatial resolution of 1.5 mm even for motion states with durations as short as 6% of the total scan time [ Fig. 8(b)]. Furthermore, the cross-correlation value (XCorr) shows excellent inverse FIG. 6. Result of bulk motion correction for MRI. T1-weighted (T1w) and T2-weighted (T2w) 3D MR images of the abdominal area were obtained in three volunteers. Volunteers were advised to carry out two to three bulk motion shifts during data acquisition which caused severe image artifacts (corrupt). The proposed motion detection and correction scheme successfully compensated for the bulk motion shifts and led to a strong improvement in image quality (corrected).
correlation with the accuracy of the obtained motion fields. Bulk motion during data acquisition also led to severe artifacts in PET images (Fig. 9). The misalignment of the sinograms caused a blurring of tumor lesions and of anatomical features such as the outline of the kidneys. Correcting for motion before combining the individual PET motion states strongly improved the final PET image. Nevertheless, the relative error could still locally reach values of up to 80% if the AC maps were not adapted for each bulk motion state. That adaption was especially important in areas where the AC values were changing locally such as along the spine which can be clearly seen in Fig. 9(a) for MCS and the difference image MCS-MF. This effect is also visible in tumor profiles where the misalignment between AC maps and image data led to a wrong increase of measured SUV values for a lesion close to the spine [ Fig. 9(b)]. For a lesion surrounded by soft tissue, such as a lesion in one of the kidneys, the adaption of AC maps to the different motion states only led to a small improvement in SUV [ Fig. 9(c)]. In summary, the best image quality was achieved if both the AC maps and the PET image data were motion corrected. Figure 10 shows SUV max and CRC measured for different lesion sizes. Bulk motion shifts of the subject led to an underestimation of uptake values of up to 67% depending on the maximum displacement and the lesion size. The difference between MCS and MCSA is clearly visible for lesions close to the spine, where misalignment of AC maps and image data in the MCS approach can distort the measured SUV values. That distortion caused errors of more than 25%.

DISCUSSION
We have successfully demonstrated a technique which provides automatic image-based nonrigid motion detection and correction for T1w and T2w 3D MR images with high isotropic resolution. In addition, the obtained motion information can be used to motion correct simultaneously acquired PET and MR images. Our approach does not require any additional acquisition of motion information but utilizes the obtained image data and can be applied to images with different contrasts. High-resolution 3D RPE data were acquired using a bitreversed order of the PE lines. This sampling scheme ensures that successively obtained radial PE lines are spread far apart and that recorded data are optimal distributed at the end of the scan. It also restricts the number of radial PE lines to a power of two. This limitation could be overcome by using a Golden-angle RPE sampling scheme. 36 The dynamics used to detect the bulk motion shifts were reconstructed using an iterative non-Cartesian SENSE algorithm. Although this approach yielded a sufficient image quality, a higher temporal resolution (i.e., more dynamic frames reconstructed using fewer radial PE lines) could be achieved by incorporating further constraints such as total variation into the image reconstruction algorithm. This approach could allow for the detection of even very fast and short bulk motion shifts.
The accuracy of the motion detection is higher if the data are split into 115 dynamic time frames rather than into 55 dynamics. Nevertheless, for the volunteer experiments we used the latter approach because it requires less computational time as each dynamic time frames requires approximately 2 min of reconstruction time using MATLAB (The MathWorks, Inc., Natich, MA). In the in vivo study, the motion state with the shortest duration and thus least amount of data contained 17.15% of the total scan. For this case, splitting the data into 55 dynamics provided sufficiently accurate motion detection.
The assessment of the accuracy of the motion estimation showed that accurate nonrigid bulk motion fields can be obtained even for motion states as short as 10% of the total scan time. Data from shorter motion states could be rejected to ensure that inaccuracies in the nonrigid motion estimation do not negatively impair the final image reconstruction. A more sophisticated approach would be to utilize the cross-correlation information between the reference and motion corrected images similar to the data weighting algorithm by Pipe 37 to exclude misregistered data which would lead to image artifacts.
The quality of the reference bulk motion state BM Ref also influences the accuracy of the obtained motion fields. Here, BM Ref was chosen as the bulk motion state with the largest amount of acquired data. Nevertheless, other selection criteria such as choosing the bulk motion state with the sharpest image features will be investigated in the future.
So far we have assumed that bulk motion shifts occur as distinct changes of the body in a time interval, which is short compared to the overall scan time. Nevertheless, the presented approach could also be used to detect and compensate for bulk motion drifts, i.e., a small but continuous change of body position during the entire scan.
As already discussed, bulk motion affects MR and PET image quality differently. Depending on the data acquisition scheme, even small and fast bulk motion can severely impair MR image quality. For PET, the effect of motion is highly dependent on the length the subject remaining in FIG. 8. Accuracy of motion detection and estimation. (a) Error E Det of motion detection for different durations of bulk motion states. E Det is given for two different cases: splitting the obtained data into 55 time frames (blue crosses, 55 dyn) and splitting the data into 115 time frames (red circles, 115 dyn). Motion detection fails for a motion state shorter than 9% (equivalent to 12 RPE lines) and 5% (equivalent to 6 RPE lines), respectively. (b) Error E Est (blue crosses) of motion estimation for different durations of bulk motion states. In addition the cross-correlation value (XCorr, red circles) between the reference image and the images transformed to the reference bulk motion state is shown. Accurate motion estimation can be carried out even for images obtained from less than 6% (equivalent to 8 RPE lines) of the total scan time (equivalent to 128 RPE lines). E Det and the duration of the bulk motion states are given in percent of the total scan time and in numbers of RPE lines. shows the best agreement compared to the MF reference. The negative effect of the mismatch between AC and emission data is highest for L1 which is close to a strong local change of the AC values (i.e., interface between bone and soft tissue). The difference in SUV ( SUV) was calculated as the relative difference of the full width at half maximum profile area between MF and noMC, MCS, and MCSA. each bulk motion state. Nevertheless, we believe that even for short movements which do not lead to immediately visible image degradation motion detection and correction is important to ensure reliable and accurate quantitative PET imaging. 38 Respiratory motion artifacts in the MR images were minimized using navigator-based gating. This ensures excellent image quality but leads to longer scan times. Buerger et al. have presented an approach to correct for respiratory motion artifacts 39 and future studies will focus on combining respiratory and bulk motion compensation to allow for a highly efficient PET-MR acquisition of the abdominal area. In the presented PET simulation, respiratory motion was not considered to ensure the obtained results demonstrate clearly the effect of bulk motion.
The SUV values and lesion profiles shown in Fig. 9 are lower than the values assigned to the segmented emission maps because the emission values were assigned in the segmentations of the high-resolution MR images. Prior to the PET simulation, these high-resolution emission maps are filtered and down sampled to achieve a realistic PET image resolution. In addition, a Gaussian filter is applied to the reconstructed PET images as a post processing step. Figure 10(b) shows that if the AC maps are not adapted to the different bulk motion states, SUV max is overestimated. This effect is due to the fact that the lesion should be corrected with a soft tissue AC value. Nevertheless, due to the bulk motion, the AC map is misaligned and the value of bone is used which leads to an overestimation of the SUV value of the lesion. If the lesion should actually be corrected using the AC value of bone but instead the AC value of soft tissue is applied, SUV max is underestimated [ Fig. 10(c)]. For both cases, using misaligned AC maps leads to artifacts around the lesion which increase the background signal. The higher background signal decreases the measured CRC. A partial or complete metabolic response to cancer treatment is defined as a change in SUV between 20% and 50% depending on tumor and study type. 2,40 The presented results of the PET simulation indicated that even bulk motion shifts less than 20 mm can lead to an error in SUV assessment of more than 40%. If the sinograms are motion corrected but the attenuation maps are not adapted to the different bulk motion states measured SUVs can deviate from the true value by almost 25%. Both errors can thus have serious impact on assessment of cancer treatment which shows the importance of detecting and correcting for even small bulk motion shifts.
While the effect of MR based nonrigid bulk motion correction has been demonstrated for MRI in healthy volunteers the impact on PET image quality was assessed using PET simulations only, which is a limitation of the presented study. Further studies are required to verify the obtained results for in vivo simultaneous PET-MR imaging. The scatter correction applied in the image reconstruction of the PET simulation is based on segmented emission and attenuation maps. For in vivo scans, true emission maps are not available and scatter correction has to be estimated from the acquired PET data and AC map. If sinograms are corrupted and the AC maps are misaligned due to motion, this approach might lead to much more severe artifacts than simulated here.
The presented technique requires a simultaneous acquisition of PET and MR data and cannot be applied to sequential PET-MR systems.
We have shown the effect of our proposed approach in bulk motion shifts of the abdominal area. Nevertheless, this method can be applied to the detection and correction of other nonrigid motion. Especially, patients who are more likely to move during a PET-MR scan, such as newborns, could benefit strongly from this technique. 41

CONCLUSION
We have presented a novel technique which detects and corrects for nonrigid bulk motion shifts ranging from 10 to 30 mm with an accuracy of 1.71 ± 0.29 mm. In addition, the obtained bulk motion fields were used to successfully compensate for the same bulk motion in a simulated simultaneous PET-MR acquisition. Our proposed method does not require any additional data acquisition but all necessary motion information is obtained directly from the acquired 3D high resolution MR image data. Furthermore, it is applicable to both T1-and T2-weighted MR imaging and its accuracy is not affected by the different image contrasts. The results of our PET simulation suggest that bulk motion can lead to an error in SUV assessment of up to 67%. Motion compensating both PET emission data and the corresponding AC maps reduced this error to below 10% compared to a motion free reference.