Combining prospective and retrospective motion correction based on a model for fast continuous motion

Prospective motion correction (PMC) and retrospective motion correction (RMC) have different advantages and limitations. The present work aims to combine the advantages of both for rigid body motion, aiming at correcting for faster motions than was previously achievable. Additionally, it provides insights into the effects of motion on pulse sequences and MR signals with a goal of further improving motion correction in the future.


| INTRODUCTION
Motion artifacts in magnetic resonance imaging (MRI) lead to repeated examinations or reduce diagnostic confidence, with negative effects on both patient comfort and economic performance of radiology centers. 1,2 Numerous remedial strategies have been proposed for selected examination types. 3 Two sequence-independent approaches for reducing motion artifacts are the prospective motion correction 4 (PMC) with external tracking devices [5][6][7] and retrospective motion correction (RMC). 7 The latter may rely either on additional motion information 8,9 or redundancies in the MRI raw data. 10,11 In a typical MRI pulse sequence with PMC, a correction of the slice or slab position and orientation is applied prior | 1285 HUCKER Et al. to the radio frequency (RF) excitation pulse, once per repetition time (TR), which for Cartesian trajectories means for every k-space line 5 or echo train. 12 Although this strategy is effective in many cases, artifacts may still appear if substantial motion occurs within a TR period. Additionally, the latency between the measurement of the object position (eg, with an external tracking device or a navigator) and the availability of corresponding quantitative motion information at the pulse sequence controller lies typically between 15 and 150 ms, 5,13,14 which is comparable to typical echo time (TE) or TR, and may even exceed both substantially for fast sequences. Figure 1 shows possible temporal relations within a PMC-enabled sequence.
Considering sequences with prolonged contrast preparation modules, it may be advantageous to apply corrections during preparation periods. For example, long diffusionweighting (DW) gradients may be split up into segments that are corrected individually. 15 However, even this procedure cannot entirely compensate for motion during an active gradient due to velocity effects and the nonzero latency mentioned above. Therefore, multishot DW EPI required additional post-processing. 16 Compared to PMC, the retrospective correction is less influenced by the latency of the motion data. Some data-driven methods 11,17,18 model the motion directly from the MR data. The RMC techniques that rely on navigator or tracking systems may apply kinematic motion models which simultaneously account for the known latency when precise time stamps are attached to the motion tracking data, synchronized with the MR scanner. Recently, several RMC approaches based on neural networks, deep learning, and other artificial intelligence (AI) techniques [19][20][21] have been proposed that provide good correction results for specific pulse sequences without requiring external insights concerning the underlying subject motion. However, RMC has no means of adapting the k-space trajectory, which remains constant in the device coordinate system, fixed at the beginning of the measurement. Therefore all rotations of the object during the measurement result in so-called "pie-slice" sampling artifacts, a term describing regions in the k-space represented in object coordinates, where no data were acquired, leading to a violation of the Nyquist criterion. 7,8,10,22 Additionally, the majority of existing RMC techniques are based either on specific pulse sequences (eg, extended to acquire navigator data) or on specific trajectories with data redundancy providing intrinsic navigation information. [7][8][9][10][22][23][24] The goal of the present study is to combine both PMC and RMC strategies, utilizing advantages of the two for enabling correction for faster motions than was previously possible. In the proposed approach, PMC is used to adapt the position of the slice and orientation of the encoding gradients in real time, whereas the residual motion between these adaptions is corrected for retrospectively. We further calculate and compensate for velocity effects during all gradient pulses, also those that do not contribute to the k-space trajectory (eg, bipolar gradients, slice rephaser, etc.) when no motion is present. The proposed framework is independent of pulse sequences and protocols; it solely requires PMC-enabled sequences and an access to pulse sequence waveforms and timing.
To enable fine temporal resolution, the tracking information is interpolated using the motion model described below. The latency between the optical acquisition and the application of the corresponding tracking data by the scanner is carefully characterized to align the clocks of both sub-systems. An accurate representation of the motion effects in the underlying MR experiment is achieved by exporting all gradient pulses, RF, and readout (ADC) events from the actual MR protocol using the vendor-provided simulator. The data are then used to calculate the effective encoding trajectory, which in turn is employed in the respectively adapted iterative image reconstruction.

| Motion and MR signal
The motion model in this study assumes a continuous steady motion of a rigid body between the subsequent tracking poses with constant linear and angular velocities and is based on Chasles' theorem. 25 Each pose comprises independent measurements of position and orientation of the body extracted from a single optical frame (or several time-synchronized shots from different cameras in case of stereoscopic tracking). Chasles' theorem states that any rigid body motion can be described by a combination of a rotation about an axis and a translation along that axis. Although translations and rotations in 3 dimensions generally do not commute, in this particular case, the translation may precede or follow the rotation or take place simultaneously. This theorem forms a basis of the so-called screw interpolation, which we use to describe the trajectory of an object between subsequently measured poses.
Translation and orientation of the object can be used to determine the trajectory of every voxel in the object assuming the imaging coordinates remain locked to the object, which is the goal of the motion correction. We define the position of each voxel at its center and denote it as � ⃗ p i (t) for the running voxel index i.
The phase of the MR signal of a voxel is calculated by Equation (1) is discretized in the time domain with intervals defined by centers of RF pulses, time points of ADC samples, and times at which gradient slew rate changes (eg, 4 points for trapezoid pulses), and for PMC-enabled protocols time points of the pose measurements. These events define the variable time boundaries for solving the integral in Equation (1). Efficient integration is achieved by the optimization steps described in the following section. Conventionally, in MRI, the encoding vector � ⃗ k is defined as the time integral of the gradients applied since the excitation of the spin system. Here, we adhere to an alternative definition inspired by the local k-space concept proposed by Noll 26 and developed further by Gallichan et al. 27 We define the effective encoding vector as the phase gradient between the neighboring voxels: where the spatial gradient operator is applied in the object coordinates. As discussed below and proven in Supporting Information, phase distributions due to a combined action of rigid body motion and pulsed gradients are limited to the first spatial order, meaning that the effective encoding vector � ⃗ k is independent of the spatial location. This means that for an arbitrary rigid object motion during MR pulse sequences, the effective k-space trajectory in the object coordinate frame may deviate from the nominal trajectory, in a matter which is identical within the entire object.

| Efficient computation of the voxel phase evolution
The gradient waveform �� ⃗ G (t) in Equation (1) can be exported from a vendor sequence simulation tool and is assumed to be known. The continuous trajectory of the given voxel � ⃗ p i (t) is interpolated based on the discrete motion tracking data using the screw interpolation. The variable time stepping defined in the previous section allows one to assume that within each time step the gradient slew rate vector remains constant and the object is moving with constant linear and angular velocities.
The screw interpolation assumes a rotation about an axis and a translation along that axis to describe rigid body motion between 2 poses. In the PMC terminology, these poses are expressed in the scanner (device) coordinate system C D with axes x, y, and z. If the screw rotation is handled in a separate coordinate system (C I ), which axes are denoted X, Y, and Z, the integration of Equation (1) can be handled efficiently. Coordinate system C I is defined such that its Z axis is aligned parallel to the screw rotation axis, then the (screw) rotation takes place entirely in the X-Y plane. C I is only rotated but not translated with respect to C D , which means that the origins of the 2 systems coincide and gradients induce no frequency offset at (X, Y, Z) = (0, 0, 0). To simplify and improve computational performance, C I is rotated about Z until the Y component of the applied gradient vanishes.
Transformation to C I can be described as C I {� ⃗ p } or as rotation R of the i-th voxel at the position � ⃗ Using the screw motion model and the new orientation of C I is advantageous: � ⃗ x i (t) can be separated into a component � ⃗ x XY,i (t) rotating within the X-Y plane and a linear motion along the Z axis, � ⃗ x Z,i (t): The rotation of the voxel within the X-Y plane is expressed as follows: .
The axis of rotation is defined unambiguously by the screw model. The representation as rotation matrix R t with the start positions X 0,i , Y 0,i in C I requires shifting the slice position to the origin and back. Therefore, it is more convenient to calculate in polar coordinates. There 2 types of variables: on the one hand the time t and the angular velocity , and the position of the rotation axis (r X0 ,r Y0 0), in C I , which are defined solely by the motion itself, on the other hand the radius of the rotation r i and the angular start offset of the defined rotation i depend on the position of the voxel inside the object.
The motion parallel to the Z axis is described by the velocity v Z and the start position along the axis Z 0,i : To model the gradient �� ⃗ G (t), a first-order description in time domain is used, which allows one to incorporate a (constant) slew rate (g x , g y , g z ) and constant gradient strength (G x , G y , G z ): The gradient vector is rotated into the integration coordinate system C I and consists of the slew rate with components g X , g Y , and g Z , and the constant gradient components G X , G Y , and G Z , respectively.
Based on Equations (1), (4), and (7), the phase accrued by the voxel i within the time step j is: Solving the scalar product results in which can be solved analytically within each time step.
Performance and accuracy of solving Equation (1) is improved by the introduced variable time stepping combined with analytical calculations within these discretized steps (compared to a trivial fixed-step numerical computation).
The resulting phase contributions ij according to Equation (10) are accumulated from the center of the excitation RF pulse to every ADC sampling time point.
For the special case when the object is not rotated between poses and therefore the point of rotation is arbitrary, only the translational motion of the screw model is used and � ⃗ x XY remains constant.

| From voxel phase to encoding trajectory
For an arbitrary motion of a rigid object according to Equations (4)(5)(6) in the presence of pulsed gradients described in Equation (7), Equation (10) is used to calculate the phase evolution of the MR signal for every voxel. As noted in the discussion of Equation (2), the resulting spatial phase distribution consists of a linear phase gradient and a global constant phase contribution. From this phase distribution, an effective encoding trajectory � ⃗ k t can be extracted according to Equation (2) and the residual will form a global phase term.
Theoretically, the calculation of the encoding trajectory is possible by integrating Equation (10) over time and space, rotating the result from C I to object coordinates and applying a spatial derivative. But C I is defined for every time step separately and hence is time dependent. Instead, our algorithm only calculates the phase of 4 points defined in the object coordinates at the beginning of the experiment as follows: the slice/FOV center, 1 mm offset in the readout, phase encoding, and slice encoding directions, respectively. The global offset is taken directly from the phase evolution of the center point, whereas the components of � ⃗ k are calculated in units of [rad/mm] from the phase differences at the corresponding offset points to the one at the slice center. The difference between the calculated global phase offset and the phase offset which was used by the scanner (depending on the slice offset) is used to correct the sampled raw data. The effective k-space trajectory is used in the nuFFT-based 28,29 image reconstruction.

| Latency
Several components contribute to different types of latencies between the optical tracking system and the acquisitions of the image, as sketched in Figure 1. A more comprehensive explanation can be found in Supporting Information.
In summary, different latencies for slice, phase, and readout encoding directions can be observed, with these relations between the axes staying constant during the pulse sequence. Additionally, there is a jitter contribution to the latency, which is constant for all encoding axes but changes from TR to TR. This work uses the timestamps, which are generated by the tracking device in step (1) in Figure 1. Therefore, imaging effects originating from these latencies can be compensated retrospectively.

| Implementation
The algorithm was implemented in Matlab (The Mathworks Inc., Natick, MA, USA). The current implementation relies on the motion tracking data of the MPT tracking system (Metria Innovation Inc., Milwaukee, WI, USA), the logging data of libXPace (PMC library for Siemens MR systems, Freiburg, Germany), and gradient data produced by the vendor simulation tool (SIM command of the IDEA framework, Siemens Healthineers, Erlangen, Germany). The implemented framework allows one to combine prospective and retrospective motion correction and offers different options for forward simulation and image reconstruction.
Sequence export with the SIM tool takes place with the protocol parameters identical to those in the MR experiment. The sequence export can optionally include exactly the same tracking data as used for PMC, captured by libX-Pace during the MR acquisition. In this case, the proposed framework is used to provide additional retrospective correction combined with the previously applied PMC. The SIM tool can also be used to export the sequence data in the stationary device coordinate system, in which case the framework can be used for pure retrospective motion correction. Further usage scenarios include evaluation and debugging of PMC sequences either with recorded or synthetic motion data.

| Verification experiments
The performance of the method was tested by imaging a rotating phantom (Figure 2, left). The cylindrical phantom is filled with a gel and can rotate about the axis of the cylinder. This axis was aligned parallel to the y axis of the MR imager with the center of the cylinder placed close to the isocenter. A coronal slice at the isocenter was imaged to reduce both B 0 inhomogeneity and gradient distortion effects. The phantom was connected to a turbine and the air flow was controlled manually to achieve desired rotation speeds which were calculated from the motion tracking logs. The MPT tracking system has a frame rate of 85 Hz and a high angular accuracy without any temporal filtering of the tracking data, which could otherwise affect accuracy and latency for fast motions. It was connected to the MR system via Ethernet and libXPace was used to update slice positions in the sequence and log the used motion data packets.
For imaging, a gradient echo (GRE) sequence was used, comprising a slice selective excitation, a phase encoding gradient, and a readout. One k-space line per TR was acquired. PMC was implemented by incorporating libXPace into the sequence with slice position updates immediately prior to the beginning of the slice selection gradient. In order to further emphasize the velocity effects, an additional bipolar gradient was introduced. For the non-moving phantom, the phase accrual due to the bipolar gradient is expected to be zero. For a rotating phantom, an additional phase offset and an encoding vector deviation ( ���� ⃗ Δk) proportional to the rotational speed is expected (see Supporting Information for derivations). Experiments were performed with different bipolar gradient strengths.
Imaging parameters for phantom experiments were: matrix size: 128 × 128, FOV: 192 × 192 mm, slice thickness: 1.5 mm, TE = 7.4 ms, TR = 100 ms, bandwidth: 400 Hz/ pixel, and readout oversampling factor: 2. All data were acquired with a birdcage Tx/Rx coil (distributed by Siemens Healthineers, Erlangen, Germany). The duration for the bipolar gradient was: 2 x 1730 µs plateau and 4 x 470 µs for the ramps. The strength of the bipolar gradient was set to achieve a gradient strength of ≈45 mT/m (using 2 orthogonal gradients of ±33 mT/m), assuming that at a rotation speed ≈100°/s two phase wraps would be seen within the phantom.
For in vivo measurements, the image resolution was changed as follows: matrix size: 256 × 256, FOV: 230 × 230 mm, and slice thickness: 2.5 mm. One normal volunteer was recruited for the study and a written informed consent was acquired prior to the scan in accordance with the protocol approved by the Ethics Committee of the University of Freiburg. Imaging in the presence of 2 kinds of motion was performed: a periodic, larger low-frequency head motion from left to right (≈20° amplitude, period time ≈1.3 s) and fast, repeated shakes of the head (≈5.4° amplitude, period time ≈0.4 s). A transversal slice orientation was chosen so that rotation occurs mainly in-plane.
Preliminary measurements revealed that the birdcage coil used is not perfectly homogeneous. B 1 inhomogeneities resulted in substantial residual artifacts when data from a large range of orientations were included in the reconstruction. For the phantom experiments, this problem was addressed by performing retrospective gating. Therefore, the experiments were carried out with 72 or 128 repetitions, leading to every k-space line to be measured repeatedly (in total: 72 × 128 = 9216 or 128 × 128 = 16384 lines). In post-processing, every line was associated with the phantom orientation when it was acquired. The complete dataset was sought for the smallest angular interval in which every phase encoding line index was found at least once.
For better comparison instead of a straightforward standard FFT, the same iterative nuFFT implementation is used as for the RMC. For this pseudo FFT (pFFT), the nominal Cartesian grid without phase correction was used. This has 3 advantages: the same processing pipeline can be used reproducing unexpected image alternations, averaging of the repeated kspace lines is performed by the nuFFT algorithm in k-space domain, and missing k-space lines (as for the optimized kspace density images) are handled automatically. We have verified that images produced by the standard FFT and pFFT from Nyquist-sampled Cartesian k-space data are identical.

| Forward simulations
For realistic simulations, tracking data from rotation measurements and the gradient waveforms which were executed by the MRI system with PMC enabled were used. The entire sequence timing including all gradients, RF pulses and ADC events were exported to generate the k-space trajectory, which was then used for forward simulations. Image artifacts were simulated based on a dataset acquired in the stationary phantom or in a volunteer without motion after an optimal shimming. Complex images reconstructed from these measurement data ( Figures 3A and 8A) were used as a virtual phantom. The simulated k-space trajectory and the global phase offset evolution were provided to the nuFFT algorithm to generate simulated frequency-domain data.

| RESULTS
In the first experiment, the phantom spinning at a constant speed (≈97.3°/s) was imaged for 128 repetitions with PMC enabled. The bipolar gradient changed for every repetition between 5 states: "off," "+PE", "−PE", "+RO", and "−RO" with a gradient strength of 33 mT/m per axis and the sign before the axis acronym referring to the first lobe. A full nonoverdetermined set of 128 k-space lines (every PE-line with a unique index exists exactly once) was selected by the described angular gating resulting in an angular interval of 10°. Due to the periodic switching of the bipolar gradients, different bipolar gradient states entered the dataset. The different states of the bipolar gradients in the presence of monotonous motion were used to model the possible phase effects from motions with changing directions. Even with PMC-enabled, image artifacts are expected due to different velocity-induced phase accumulation effects for different k-space lines. Based on the effective encoding trajectory, a forward simulation of the expected PMC performance is presented in Figure 3B. For these simulated data and the raw data from the actual PMC experiment ( Figure 3C), the pFFT was applied for image reconstruction assuming the nominal Cartesian k-space trajectory resulting in strong image artifacts. The combination of PMC and RMC substantially improves the phantom appearance as seen in Figure 3D. To investigate the possible sources of the residual artifacts the k-space sampling density distribution was calculated by counting the number of k-space samples per rectangular bin on a Cartesian grid. This and all following sampling density diagrams use the identical color coding with black corresponding to missing samples.

F I G U R E 2 A,
The phantom is the cylindrical container (ID = 109 mm, H = 46 mm) filled with hydroxyethyl cellulose gel mixed to high viscosity (water, 11% hydroxyethyl cellulose, 0.2% copper sulfate) and suspended on ceramic ball bearings allowing it to rotate about the axis of the cylinder, belt drive, and windmill chamber. The phantom is connected via a belt transmission to a turbine and the air flow through a nozzle can be manually controlled to achieve the desired rotation speeds. The marker is fixed in the top at the rotational axes. B, The marker fixation mouthpiece for the in vivo measurements. The mouthpiece was fixated to the upper jaw for rigid body motion Figure 3E shows numerous black stripes and wedges (pieslice artifacts) where no samples are available and therefore the Nyquist condition is violated.
All encoding lines from the same angular segment (453 in total) were used to calculate the images in Figure 4. By using all lines, it was possible to reduce Nyquist violations, as evidenced by the sampling density plot in Figure 4D and a vastly improved image quality for PMC+RMC images ( Figure 4C). It is noteworthy that the increase in the number of k-space lines used for the pFFT reconstruction did not change the appearance of the magnitude images in Figures  4A,B in comparison to the corresponding images in Figures  3B,C. To evaluate whether the artifact reduction observed in Figure 4C was not due to averaging of the repeated k-space lines, another set of 128 lines was extracted from the total of 453 lines, achieving an optimal k-space coverage. Figure 5D shows that the sampling density is more evenly distributed than in Figure 3E and the number of gaps in this case is almost as low as in Figure 4D. Correspondingly the amount of ghosting in Figure 5C is substantially reduced compared to Figure 3D.
Another phantom experiment was conducted with a variable rotation speed and a fixed bipolar gradient (±33 mT/m simultaneously applied in readout and phase directions) with 72 repetitions. As in the previous experiment, gating based on the angular position was performed resulting in an interval of 16.4°. The speed of the phantom for the selected 128 lines ranged from 63.4°/s to 280.36°/s (avg: 163.4°/s; std: 45.2°/s). The images in Figure 6 were reconstructed in the same way as in the previous experiment. The overall image quality is similar to that in Figure 3 but PMC+RMC ( Figure 6C) contains in addition to ghosting, even stronger signal void artifacts. In this particular experiment, the k-space sampling pattern in Figure 6D shows obvious pie-slice effects, with a rather large sampling gap located close to the k-space center, which may explain the observed signal void in the image.
Similar to the first experiment, an optimal set of 128 lines from the same angular interval was extracted optimizing the k-space sampling density, with angular speeds ranging from F I G U R E 3 The rotating phantom was imaged with PMC, with a constant angular velocity and a variable bipolar gradient in subsequent repetitions. Angular gating was applied to collect a complete k-space set of 128 lines with the repeated lines discarded resulting in a gating opening-angle of 10.0°. Due to the variable bipolar gradient, varying velocity effects in different repetitions entered the same k-space dataset producing strong artifacts in the standard PMC imaging. A, Phantom imaged after shimming and without rotation and without bipolar gradient, used as a virtual phantom for subsequent forward simulations. B, A simulated image showing the expected PMC performance was generated by using the calculated effective encoding trajectory to generate k-space sampling points from image A using nuFFT followed by a pseudo FFT reconstruction. C, The experimental results of a prospectively motion corrected experiment. D, Reconstruction from the same data as C using the effective encoding and an iterative CG-nuFFT reconstruction. E, The sampling density of the effective encoding trajectory, with the following color coding: black=0; green=1; yellow=2; white=3+ sampling points per nominal sampling grid bin | 1291 HUCKER Et al.
69.9°/s to 246.9°/s. Figure 7D shows a more uniform sampling distribution and correspondingly the PMC+RMC reconstruction in Figure 7C displays a much better image quality.
The images from the first in vivo experiment are shown in Figure 8, acquired in the presence of left-right rotational head motion and without a bipolar gradient. Figure 8A is a reference image measured without motion, used for the forward simulation in Figure 8B. Noteworthy is the F I G U R E 4 From the same experiment as in Figure 3, all k-space lines within the same 10° gating angle interval were used (453 lines in total). A, Simulated image showing the expected PMC performance based on the virtual phantom in Figure 3. B, The corresponding experimental result with PMC. The oversampling in both cases was removed prior to the pseudo FFT by averaging the k-space lines with identical indexes. C, The result of the combined PMC and RMC applied to the experimental data in B. D, The sampling density of the effective sampling trajectory, with the same color coding as in Figure 3 F I G U R E 5 From the 453 k-space lines used in Figure 4, 128 lines were selected to produce most homogeneous k-space sampling. A, Simulated image showing the expected PMC performance. B, The results of the PMC experiment. C, The results of the combined PMC and RMC reconstruction. D, The sampling density of the effective encoding trajectory, with the same color coding as in Figure 3 similarity of this image based on the simulated k-space data to the image reconstructed from the pure PMC measurement data shown in Figure 8C. To be able to better appreciate the accuracy of forward simulations, two areas showing prominent artifacts were zoomed in: a diagonal black stripe in the prefrontal area in Figure 8E and the double folding artifact originating from the superior sagittal sinus in Figure 8F. The second in vivo experiment was performed with the same type of motion but with the bipolar gradient (±33 mT/m applied simultaneously in readout and phase direction). To achieve a sufficient k-space sampling density, we acquired 3 repetitions. For the forward simulation in Figure 9A (based on 8A) and the PMC measurement Figure 9B, the 3 repetitions were averaged by the pFFT. For the image reconstruction in Figure 9C, k-space lines from all 3 repetitions were used in the RMC.
In the last in vivo experiment, the amplitude of head rotations was reduced, but the frequency increased. Additionally, we have reduced the bipolar gradient strength to ±16 mT/m applied simultaneously in readout and phase directions. Figure 10 shows the artifact reduction achieved by the combination of the PMC and RMC if the k-space sampling density is increased. Images in Figure 10A are created from raw data of 1 repetition, B by 2 repetitions, C by 3, and D by 4, respectively. The second repetition reconstructed separately shows similar image quality as Figure  10A (see Supporting Information Figure S5). Comparing Figure 10A to B,C, or D, one can observe reduced artifacts upon the increase of the data amount, where PMC+RMC rapidly recovers the image quality close to the original, whereas PMC still suffers from substantial ghosting and prominent signal void artifacts. Sampling densities for Figures 10A-D are provided in Supporting Information F I G U R E 8 Experimental results with continuous head motion without bipolar gradients. A, nonmoving reference image. All other images correspond to repeated left to right head rotations with the peak rotation of ≈±10° with ≈0.8Hz and a peak angular velocity of ≈±50°/s. See Supporting Information Figure S3 for detailed motion plots. B, Simulation of the predicted PMC performance. C, PMC experiment. D, Combined PMC and RMC reconstruction of the same experimental data. Images (E) and (F) show zoom in views of the characteristic artifacts Figure S6. A similar experiment without bipolar gradient can be found in Supporting Information Figure S7, with weaker artifacts but a comparable correction efficacy.

| DISCUSSION
Subject motion during MRI examinations remains to be a factor substantially limiting the achievable image quality in clinical settings and image resolution and measurement time in advanced research applications. 3,30 For brain imaging in normal volunteers, both PMC and numerous retrospective correction methods offer several viable options to circumvent the adverse effects of involuntarily motion. 6,7,9,31 However, no data are available in the literature to date concerning the motion correction options for patients with neurologic motion disorders such as essential tremor (ET) and Parkinson's disease (PD). ET is a common neurological disorder. In the classic ET, the rate of tremor may vary between 4 and 12 Hz and is characterized by a kinetic or postural tremor. It can occur in all parts of the body (≈1/3 cases in the head). 32,33 PD is the second most common neurodegenerative disorder after Alzheimer's disease with a core symptom of resting tremor 34 (in 60-70% PD patients) and gradually aggravates with the progression of disease. 35

F I G U R E 9
Experimental results with continuous head motion and bipolar gradients. Motion paradigm was similar to that in Figure 8. See Supporting Information Figure S3 for detailed motion plots. The strength of the bipolar gradient was 33 mT/m applied simultaneously on both phase and read axes. The raw data of 3 repetitions were combined. A, Simulation of the predicted PMC performance based on Figure 8A. B, PMC experiment. C, Combined PMC+RMC reconstruction of the same experimental data F I G U R E 1 0 Experimental results with bipolar gradients and faster continuous head motion of reduced amplitude (full amplitude ≈5° at ≈3-3.5 Hz, see Supporting Information Figure S4 for detailed plots). The amplitude of the bipolar gradient was 16 mT/m applied simultaneously to phase and read axes. Left column shows PMC/pFFT reconstruction, right column combined PMC+RMC reconstruction. A, Data of first repetition. B, Data of first and second repetition (see Supporting Information Figure S4 for data of only second repetition). C, Data of 3 repetitions combined. D, Data of 4 repetitions combined. The corresponding k-space sampling densities are plotted in Supporting Information Figure S6 | 1295

HUCKER Et al.
A combination of the high frequency of tremor, the limited sampling rate and latency of the current motion tracking systems is expected to present a substantial challenge to PMC. Even if the head position is continuously monitored during the scan, and the MR acquisition parameters are updated as soon as the tracking data arrive to the sequence controller, the requirements to the effective head motion correction for tremor patients are expected to go beyond the trivial position updates. 36 The present study lays out an important ground stone toward extending the camera-based motion correction to such patients. It explores the 2 possible physical mechanisms leading to the data corruption in the presence of fast motion. The first is the velocity effects during the large gradient lobes such as prephaser gradients or crushers, which is addressed and investigated by means of introducing an additional bipolar gradient to the gradient echo sequence. The second group of effects is related to the limited temporal rate of sampling of motion and the unavoidable latency in the prospective correction chain.
In order to address both of the abovementioned mechanisms, we have developed a framework capable of replicating offline all gradient, RF and ADC events of the PMC-enabled sequence based on the tracking logs captured during the experiment. This exhaustive information is then combined with additional tracking information and the continuous motion model to calculate the actual phase evolution in the object's frame of reference. As a model a 2D spoiled gradient echo sequence was selected for its simplicity and robustness.
The validity of the framework is supported by the fact that the forward simulations in the Figures 3-9 produce similar artifacts as those observed in the actual measurements. Especially in Figure 8 ghosting artifacts are reproduced particularly closely. Minor differences between the simulation and measurement can still be observed in some cases, most notably in Figure 9. The possible reasons for such deviations are discussed below. In general, the forward simulation performs sufficiently well for making qualitative predictions about the artifact levels based on a given motion trajectory, the defined sequence and the used protocol. This makes it suitable to compare different sequences, or to optimize protocols based on the previously captured motion trajectory.
However, the ability to accurately predict image artifacts in the forward simulation does not necessary imply that artifactfree images can be reconstructed from the measured k-space data corrupted by motion. In addition to the local variations of the sampling density in k-space, there are several effects in the in vivo experiments, which are not accounted for by the presented method. The most prominent of these the RF coil inhomogeneity, spacial variation of residual shimming imperfections or trough-plane motion ignored in our current 2D correction implementation. When the k-space trajectory is oversampled all of these effects can be partially compensated for, making it difficult to attribute the residual artifacts to one of the possible causes. But, if the intended k-space trajectory was not overdetermined, rapid motion is likely to result both in local violations of the Nyquist sampling criterion and increased data inconsistencies due to the above-mentioned additional imperfections, leading to ghosting and signal modulation in the images. Therefore, although in some measurements combining PMC and RMC could directly improve the image quality, as in Figure 8D, in other experiments as Figures  3D or 6C the nonuniform density distribution in k-space led to persistent artifacts. In case of phantom experiments, we attribute the residual artifacts to the Nyquist sampling violations because the additional effects were effectively suppressed by pure rotational motion and angular gating. In case of Figure  10A, nonuniform density distribution is only a minor component among multiple effects. All these residual artifacts could be reduced either by selecting a different subset of k-space lines leading to a more optimal sampling density, as in Figures  5 and 7 or by including additional k-space lines, as in Figures 4,9 and 10. In particular, Figure 10 explores the influence of oversampling on the efficacy of the combined correction for strong motion. As seen, doubling the amount of the available data from Figure 10A to B leads to a strong artifact reduction. Adding another full set of k-space lines ( Figure 10C) continues to reduce artifacts. A further increase of the data amount, however, does not lead to a significant improvement. Noteworthy is the possibility of combining multiple short measurements incrementally, for example, measure 1 repetition with PMC as in Figure 10A and if the image quality is insufficient additional equivalent measurement(s) can be triggered. With the accurate forward model multiple measurements, each with insufficient image quality, can be combined retrospectively as in Figures 9C or 10B. If strong motion is expected it may be beneficial to use phase oversampling from the beginning to increase the sampling density.
Another possibility, and also a necessity for future clinical usage, would be to apply parallel imaging for reconstruction with these undersampled k-spaces. Spatially varying coil sensitivities are expected to produce artifacts in the presence of motion if ignored, but can also be turned into advantage if properly accounted for. Although techniques for parallel imaging with rigid body motion were proposed in the literature 37,38 none of them are designed for continuous motion. Also compressed sensing techniques 39 could be used to handle the sampling gaps in k-space. Gradient nonlinearities can be compensated for by adapting corresponding reconstruction algorithms. 40 In the current work, we assumed that the tracking data accurately represent head poses at any time and are only affected by temporal deviations. The proposed approach can, however, be trivially extended to account for further tracking errors such as tracking noise, 41 cross-calibration errors 14,42 or marker fixation problems. 43 But all these possible extensions are expected to increase the complexity and reduce the robustness of the presented framework. They might be added in future upon a careful debugging and extensive validation.
Another option is to combine the proposed framework with current machine learning strategies [19][20][21] to suppress residual artifacts. Alternatively, the framework can be used to create realistic and rich training data for such methods based on the simulated or measured motion trajectories. Due to the ability of the framework to reproduce both magnitude and phase of the images and its sensitivity to the gradient timing and waveforms, such AI methods can be made more aware of the underlying physical effects and the specifics of the particular sequences, potentially improving their robustness and correction efficiency.
Although presented and validated in 2D experiments, our RMC framework is expected to perform better with 3D sequences. Despite of the increased computation, the overall effort remains tractable (see Supporting Information on computational performance). Sole image reconstruction times comparable to undoing PMC 44 are expected, that can be further accelerated by graphic processing units (GPUs). 45 For 2D imaging, through-plane residual motion may cause a shift of the excited or pre-saturated volume within the object leading to spin-history effects. 46 Residual through-plane motion becomes relevant when the local misalignment of the slice profiles becomes comparable to a fraction of their thickness. A slice thickness in clinical MRI typically substantially exceeds in-plane voxel dimensions. Although in our experiments the subjects were instructed to mainly perform in-plane motions, also substantial through-plane rotational components were present (see Supporting Information Figures S3  and S4) which the proposed RMC was able to suppress, probably due to the above-mentioned reasons. Further research may be needed for optimizing through-plane effects for 2D imaging for target subject populations and specific protocols.
For the current state of knowledge, it is difficult to define a clear criterion concerning the patient groups and motion types, for which the described framework may be necessary. Based on the results in Figures 8-10 and Supporting Information Figures S4, S5 and S7, for repeated oscillatory motions with peak angular velocities of approx. 50°/s both substantial residual artifacts can be seen in PMC images and the proposed RMC is effective.
Multi-spin echo-train sequences are of high clinical relevance, but accumulated motion-induced phase and changes of the encoding trajectory are expected to violate the Carr-Purcell-Meiboom-Gill (CPMG) condition. By adding the concept of extended phase graph (EPG 47 ) to the presented work this violation may potentially be quantified for every k-space line optionally followed by appropriate compensation strategies.
From our observations (see Supporting Information), calculating the k-space trajectory on the fly directly within the scanner should be achievable and would greatly simplify image reconstruction, especially for non-Cartesian encoding trajectories. Once this basic infrastructure is available it can be trivially extended to account for bulk motion in flow or diffusion weighted imaging.

| CONCLUSION
The presented framework enables to calculate deterministic effective encoding trajectories in the presence of rigid body rotations and translations, shows good results both for forward simulation of the motion artifacts and as a input for the retrospective motion correction. It can further provide a powerful debugging platform for simulating the motion artifacts in various sequences under simulated or realistic motion paradigms.
Supporting Information is available as separate document. It contains additional sections about the detailed latency contributions, experimental limitations and requirements and the computational performance, as well as Supporting Information Equations and Supporting Information Figures S1 to S7 FIGURE S1 Simulation of the influence of the pie-slice effects. Motion data from the PMC experiment with variable speed and a bipolar gradient. Every k-space line was used exactly once (128 lines in total). A, The measured phantom without motion was used as image source for the forward simulation. B, Simulated data were used in pFFT reconstruction. C, Iterative reconstruction simulating the ultimate PMC+RMC performance. As seen, local violations of the Nyquist condition lead to severe artifacts FIGURE S2 Simulation of the compensation of the pieslice effects by using optimized k-space density distribution. Motion data from the PMC experiment with variable speed and a bipolar gradient. 128 k-space lines selected based on optimized density distribution. A, The measured phantom without motion was used as image source for the forward simulation. B, Simulated data were used in pFFT reconstruction. C, Iterative reconstruction simulating the ultimate PMC+RMC performance. Compared to Supporting Information Figure  S1 artifacts are strongly reduced.

FIGURE S3
Representative fragment of the motion logs from the first in vivo experiment plotted in MRI device coordinates, (see Figure 8 for corresponding MR images) FIGURE S4 Representative fragment of the motion logs from the third in vivo experiment plotted in MRI coordinates, (see Figure 10 for corresponding images). FIGURE S5 Images reconstructed solely from the data of the second repetition. The corresponding k-space data were combined with the k-space data of Figure 10A to generate Figure 10B FIGURE S6Density plots of Figure 10.

FIGURE S7
The experiment shown in of Figure 10 was also performed without a bipolar gradient. A similar head motion paradigm was applied-as far as it is possible for volunteers to replicate the same motion twice. A, Data of first repetition; B, data of first and second repetition; C, data of 3 repetitions combined; D, data of 4 repetitions combined