MotionTTT: 2D Test-Time-Training Motion Estimation for 3D Motion Corrected MRI

Authors: tobit  klug
kun  wang
stefan  ruschke
reinhard  heckel
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 19:48:27

Abstract: A major challenge of the long measurement times in magnetic resonance imaging (MRI), an important medical imaging technology, is that patients may move during data acquisition. This leads to severe motion artifacts in the reconstructed images and volumes. In this paper, we propose MotionTTT a deep learning-based test-time-training (TTT) method for accurate motion estimation. The key idea is that a neural network trained for motion-free reconstruction has a small loss if there is no motion, thus optimizing over motion parameters passed through the reconstruction network enables accurate estimation of motion. The estimated motion parameters enable to correct for the motion and to reconstruct accurate motion-corrected images. Our method uses 2D reconstruction networks to estimate rigid motion in 3D, and constitutes the first deep learning based method for 3D rigid motion estimation towards 3D-motion-corrected MRI. We show that our method can provably reconstruct motion parameters for a simple signal and neural network model. We demonstrate the effectiveness of our method for both retrospectively simulated motion and prospectively collected real motion-corrupted data. Code is available at \url{https://github.com/MLI-lab/MRI_MotionTTT}.

Main paper: 
 introduction: 
 Magnetic resonance imaging (MRI) is one of the most important medical imaging technologies due to its non-invasiveness and ability to foster diagnosis of a wide range of diseases. However, its inherently long scan times make MRI susceptible to motion artifacts caused by patient movement during the scan. Repeating scans corrupted by motion artifacts causes additional costs and reduces patient through-put, and unnoticed artifacts can lead to misdiagnosis [2, 36] .We consider the problem of imaging under motion and propose an algorithmic solution to correct for the motion based only on the measurements acquired during the scan, without requiring additional hardware or changing the measurement process, or interrupting the clinical workflow.A traditional approach to algorithmic motion reconstruction is to jointly estimate the motion parameters and motion-corrected reconstruction [7, 6, 15] , but those methods are slow and can be inaccurate in particular for severe motion.Deep learning based approaches have been proposed to accelerate reconstruction and potentially allow to account for more severe motion. However, most existing data driven methods correct for in-plane motion within 2D MRI (see the review Spieker et al. [38] ), as 3D data for training 3D models is only scarcely available and computationally expensive to handle [18] . In practice, however, motion occurs in 3D and not in-plane. Moreover, the duration of a scan in 3D is significantly longer than in 2D making motion more likely and thus motion reconstruction more important. Finally, data-driven • We use retrospectively simulated motion to demonstrate the ability of our method to accurately estimate motion over a wide range of motion severities. Combined with a L1minimization reconstruction module we achieve effective 3D imaging under motion and outperform a classical alternating optimization baseline [7] in terms of estimation speed and estimation performance under severe motion and a deep learning based end-to-end motion correction baseline Al-Masni et al. [1] in terms of reconstructed image quality.• We demonstrate the potential of our method on prospectively acquired real motion-corrupted data achieving significant improvements in terms of visual image quality. 
related work: 
 Approaches for retrospective rigid-motion correction for MRI can be categorized into supervised deep learning-based approaches, model-based optimization approaches, and combinations thereof.Supervised deep-learning based approaches learn a mapping from undersampled measurements corrupted with simulated motion to the motion-free reference image and have been proposed for 2D [21, 24, 42, 35] and 3D MRI [18, 10, 1] . However, as demonstrated in our work the images produced by those end-to-end approaches are often blurry and the methods pertain to the type of motion simulated during training [23, 14, 16] .Alternating optimization [7] is a classical model-based approach for joint motion estimation and correction in 3D MRI, where every iteration alternates between optimizing over the motion parameters while fixing the estimated reconstruction and vice versa. However, jointly optimizing over both unknowns without any prior information is a highly complex optimization problem resulting in long run times and errors in the presence of more severe motion as demonstrated in our work.The speed and robustness of alternating optimization can be improved through augmenting either the reconstruction step or the motion estimation step with deep learning [14, 16] . However, both methods have so far only been proposed for 2D in-plane motion correction and rely on synthetic motion simulation during training, which makes them specific to the type of simulated motion. In practice, motion always occurs in 3D.For motion estimation, Singh et al. [34] pre-trains a neural network to predict a motion-free image from the motion-corrupted measurements conditioned on the true motion parameter. At inference the data consistency loss is optimized only over the motion parameters. However, the motion estimation relies on learning the characteristic of motion-corruptions. Hence, the model operates in the measurement space in which the corruptions occur, which in multi-coil MRI is of much higher dimension than the corresponding image space, thus making an extension of this approach to 3D very challenging. Our method relies on learning the characteristics of motion-free data in the image space and thus is efficent for 3D imaging as we show.Levac et al. [22] proposes a method based on diffusion models trained on generating motion-free images for 2D motion reconstruction, that also does not rely on motion simulation during training. At inference, sampling based reconstruction [17, 5] with joint motion trajectory estimation is performed.While our approach as well as the aforementioned works correct motion artifacts solely based on the acquired MRI measurements, another line of research corrects for motion prospectively or retrospectively based on additional data collected during the scan via, e.g., external detectors [30, 27] or navigator sequences [40, 43, 11] . However, those methods usually require interference with the standard clinical measurement process and are often tailored to a specific measurement sequence or setup limiting their broad applicability in practice.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] Justification: Limitations regarding the performance of the reconstruction module as well as the dependency on fully-sampled training data are discussed in Section 7.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
problem statement: 3d mri imaging under motion: 
 A 3D multi-coil accelerated MRI measurement y ∈ C C×kx×ky×kz is obtained byEQUATIONwhere A = MFE is the forward model, x ∈ C rx×ry×rz the object of interest, and z measurement noise. The measurement y consists of C-many k-space measurements collected by C coils. The expand operator is defined as Ex = [S 1 x, . . . , S C x] with coil sensitivity maps S j . The 3D Fourier transform F and undersampling mask M are applied coil-wise. We consider 3D Cartesian undersampling, where undersampling takes place in the plane of the two phase encoding dimensions k x × k y and the frequency encoding dimension k z is fully sampled (see Figure 1 ).We focus on rigid motion, where the i-th motion state is defined by three translation and three rotation parametersm i = [t i 1 , t i 2 , t i 3 , ϕ i 1 , ϕ i 2 , ϕ i 3 ]. MRI acquisition under motion can be described asEQUATIONwhere the forward model A is a function of the unknown motion states m = [m 1 , . . . , m b ], where b is the number of motion states, and of the known MRI sequence specific sampling trajectory T that specifies which part of the k-space data y is acquired when.The goal of this work is to reconstruct the volume x from the undersampled, motion-corrupted measurement y without knowledge of the motion states. We do so by estimating the motion states m from the measurement y and sampling trajectory T and use the estimated parameters to reconstruct a motion-corrected volume x.Parameterization of the forward model under motion. In practice, a measurement is acquired in batches of lines in the k-space along the frequency encoding dimension k z . A batch, referred to as shot, is acquired within a short time window followed by a pause before the next subsequent shot. It is popular to assume that a subject's position is constant during one shot and motion happens during the pause between shots. This is known as inter-shot motion [7, 6, 18, 22, 34, 16] . For inter-shot motion, the number of motions states b is equal to the number of shots and the sampling trajectory T maps the lines in the k-space acquired during the i-th shot to the i-th motion state. See Figure 1 (c, d ) for examples of sampling trajectories used in practice.However, in practice, motion can occur anytime, and thus the inter-shot motion introduces an approximation error. Motion during the acquisition of one shot is referred to as intra-shot motion [14] .Then, each k-space line acquired during such a shot can have a a distinct motion state. In this work, we investigate the capabilities of our method under both inter-and intra-shot simulated motion.Within the forward model A(T , m), motion corruption can be applied in the image or in the kspace [25] , since rotations and translations in the image space translate to rotations and linear phase shifts in the Fourier space. We apply motion corruption in the k-space, because it is computationally more efficient for our setup (see Appendix A). We use the non-uniform FFT (NUFFT) N(T , ϕ) to sample fully-sampled k-space data at the rotated coordinates for each shot. Translations are applied via linear phase shifts L(T , t) to obtain motion-corrupted k-space data.As a starting point for reconstructing the volume x from an undersampled measurement y with zero-filled (ZF) missing entries, it is common to compute a ZF reconstruction x † = A † y, whereA † y = C j=1 S * j F -1 y j .If the measurement y is motion corrupted, a corrected ZF reconstruction based on motion parameters m is x † = A † (T , -m)y. First, translations are reverted with a phase shift in the opposite direction L(T , -t). Then, rotations are corrected for via the adjoint NUFFT N adj (T , -ϕ). See Figure 2 for an illustration of the forward model and ZF reconstruction.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Justification:Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
motionttt: 
 The proposed MotionTTT consists of: 1) pre-training a neural network for 2D motion-free image reconstruction, 2) test-time-training to estimate motion parameters of the motion-corrupted 3D k-space data and 3) reconstructing the motion-corrected 3D image based on the estimated motion.Step 1: Pre-training. Given motion-free data {(x 1 , y 1 ), . . . , (x N , y N )} consisting of pairs of 2D reference images x ∈ C rx×ry and undersampled k-space data y ∈ C C×kx×ky , we train a U-net [32] f θ with weights θ to map a ZF reconstruction A † y i to the image x i by minimizing the lossL train (θ) = N i=1 |f θ (A † y i )| -|x i | 1 /∥x i ∥ 1 + FEf θ (A † y i ) -FEx i 1 /∥FEx i ∥ 1 . (3)We use this combined training loss between magnitude images and k-space data since it leads to better performance for motion-free reconstruction than using one of the individual losses (see Appendix C).Step 2: Test-time-training for motion estimation. Given an undersampled and potentially motion corrupted 3D measurement y and a sampling trajectory T we freeze the weights θ of the trained network f θ and estimate the motion parameters m by minimizing the data consistency lossEQUATION) with Adam [19] . The idea behind minimizing this loss is as follows. If applied to motion-corrupted data, at initialization with m = 0 the motion correction A † (T , -m) has no effect and the motioncorrupted network input results in a large loss as the network was trained on reconstructing motion-free data. Contrary, when the motion parameters are chosen correctly the network input is a motioncorrected ZF image, which is similar to a motion-free ZF image and results in a small loss. See Figure 2 for example images. In Section 5 below we study this loss theoretically.We call this approach test-time-training, since the adjoint A † (T , -m) can be considered to be part of the network, and by optimizing over the motion states m we are optimizing over part of the network's parameters. Methods that optimize a network at inference are referred to as test-time-training methods, and are successful at prediction under distribution shifts [39, 9] .In every iteration, the network reconstructs the entire 3D input volume slice-wise, where the slicing direction is sampled uniformly at random to be either in the r x × r y , r x × r z or r y × r z image plane. We compute gradients only for a subset (of size 5, limited by GPU memory) of slices sampled independently in every iteration. While motion has a local effect in the k-space, the artifacts spread globally in the image space hence computing gradients with respect to a single slice can contain signal about all motion parameters.In order to minimize the loss (4) reliably over different levels of motion severity the optimization scheme is important. We take a three-phase optimization approach.Phase 1 optimizes over one motion state per acquired shot. We start with a large initial learning rate in order to explore the non-convex loss landscape. Especially for strong motion initializing all parameters with 0 can lead to a large distance to the true motion parameters. During phase 1 the learning rate is decayed twice in order to converge to a stable first estimate of the motion parameters.At the start of phase 2 we compute the DC loss (4) for every estimated motion state miEQUATIONwhere the mask M mi keeps only the part of the k-space acquired during the i-th state. Motion states with a loss larger than a certain threshold are likely estimated poorly and we reset them to the average between the previous and next motion state that fall below the threshold. Phase 2 then only optimizes over the motion states that have been reset. For intra-shot motion, a thresholded motion state is not only reset, but N splits additional motion states up to the number of acquired k-space lines per shot can be introduced in order to estimate a more resolved motion trajectory for this shot.Phase 3 again optimizes jointly over all motion states with a small learning rate to converge to a final estimate of the motion trajectory.Step 3: Reconstruction. Using the estimated motion parameters m, we can obtain an estimate of the motion-corrected image directly from the network output f θ (A † (T , -m)y), and apply a data-consistency step to improve performance (U-net-DCLayer). This step moves the frequencies of the reconstructed image closer to the given frequencies, and is used for example by Chen et al. [4] . Alternatively, we can use any reconstruction method for motion-free data, such as a classical L1-minimization-based reconstruction [26] . We refer to L1 or U-net reconstruction based on motion parameters estimated by MotionTTT as MotionTTT-L1 and MotionTTT-U-net-DCLayer and their respective performance with the oracle known motion parameters as KnwonMotion-L1 and KnownMotion-U-net-DCLayer. We refer to the reconstruction after DC loss thresholding that excludes motion states with a large DC loss from the reconstruction as MotionTTT+Th-L1. It can be seen that there is a sharp minima around m = m * . This minima turns out to be unique under certain conditions. In our theory we consider discrete shifts, indicated by crosses.a = 0 a = 0 a = 1 a = 1 a = 2 a = 2 a = 3 a = 3 
theory for motion ttt: 
 We consider the following model to illustrate the principle of our method. We consider a signal x ∈ R n the lies in a d dimensional subspace described by the matrix U ∈ R n×d , i.e., there is a coefficient vector c ∈ R d so that x = Uc. We take the matrix U as a random Gaussian matrix with iid N (0, 1/ √ n) entries, so that the columns of the matrix are approximately orthonormal.Let F T be the Fourier matrix with rows chosen in the set T ⊆ {0, . . . , n -1}. We assume a measurement model where the signal x is shifted by unknown discrete integer parameters m * 1 , . . . , m * b ∈ Z, and for each shifted version of the signal, a set of measurements is collected according toEQUATIONwhere D m,T ℓ is a diagonal matrix with e i2πmj/n , j ∈ T ℓ on its diagonal. Note that this multiplication with complex exponentials in the frequency domain implements a circular shift in the time domain.In this measurement model, the signal is assumed motion-free while the measurements in the set F T ℓ are collected. The frequencies in the set T ℓ are chosen by sampling each frequency independently with probability k/n. So in expectation, k frequencies are included in the set T ℓ .We consider the network f (x) = n bk UU T x for reconstructing a clean signal, where bk is the total number of measurements collected. This choice of network is motivated by the fact that if the measurement is not motion corrupted, then we have that (see Appendix B.3)EQUATIONwhere T = T 1 ∪ . . . ∪ T b is the set of all measurements collected and (•) * denotes the complex conjugate. We consider our test-time-training loss for this model, which isEQUATIONBelow, we show that for our model, under certain conditions and with high probability, the loss has a unique minimum at L(m * ). Before stating our result, we visualize the loss in Figure 3 . It can be seen that the loss is not convex in m which makes it difficult to optimize.Theorem 1 Consider the model introduced above, and assume that the signal x = Uc is chosen randomly by drawing the entries of c iid from a zero-mean unit-variance Gaussian distribution. Let a(m) be the number of values of m 1 , . . . , m b that are non-equal to m * 1 , . . . , m * b . The following statement holds for all m ∈ {0, . . . , n -1} \ {m * } simultaneously with high probability: IfEQUATIONthen L(m) > L(m * ), where c is a numerical constant.The theorem implies that if the subspace dimension, d, and the number of shifts, b, are sufficiently small relative to the number of measurements, bk, then the loss has a global minimum at the true shift m * . 
experiments: 
 We demonstrate quantitatively and qualitatively on simulated data the ability of our method to accurately reconstruct images for a wide range of levels of motion severity in the presence of inter-and intra-shot motion. Moreover, on prospectively acquired real motion-corrupted data, we demonstrate that our method achieves significant gains in terms of visual reconstruction quality.We start with experiments with simulated inter-shot motion as described in Section 3. The model, data, and baselines considered for both inter-and intra-shot motion experiments are as follows.Model. We use a 17.5M parameter U-net [32] f θ , a common baseline with good performance on image-to-image tasks [46, 3, 13] , training details are in Appendix D.Data. We train the U-net f θ on motion-free data and evaluate MotionTTT with simulated motioncorrupted data sourced from the Calgary Campinas Brain MRI Dataset [37] (license CC BY-ND) consisting of 3D scans of size k x × k y × k z = 218 × 170 × 256 and C = 12 receiver coils. We select 40 subjects for training from the training set, 4 subjects for validation and hyperparameter tuning, and 10 subjects for testing from the validation set. We train and test with an undersampling factor of 4 using the mask from Figure 1 (c ). For training, we slice each 3D zero-filled reconstruction A † y and the corresponding 3D reference volume x along all three dimensions resulting in about 25k pairs of 2D network inputs and targets. We compute sensitivity maps from 24 × 24 × 24 auto-calibration lines of the originally fully-sampled motion-free k-space with ESPIRiT [41] .Motion simulation. We set the number of shots to B = 50 similar to how our own real data was acquired in the upcoming Section 6.3, and we use an interleaved sampling trajectory T (Figure 1 (c)), where every 50-th line in the k-space is acquired within one shot and the 3 × 3 center of the k-space containing the largest energy is sampled in the first shot. Without loss of generality, we assume the subject to be in zero-motion state at the first shot and hence do not simulate and estimate motion parameters for the first shot.As the characteristics of patient motion vary widely, synthetic rigid body motion is often simulated as random motion with rotations and translations drawn uniformly from some range, or from a Gaussian [7, 22, 16, 34] .We simulate random motion with different levels of severity by varying the number of motion events N e ∈ {1, 5, 10} per scan and the maximum possible rotations/translations M max ∈ {2, 5, 10} in degrees/mm. The shots between which a motion event occurs are sampled uniformly at random, and for each event translation and rotation parameters are sampled uniformly from [-M max , M max ]. After a motion event the subject stays in this position until the next motion event occurs. This yields 10 different levels of motion severity including the motion free case.Baselines. We compare to alternating optimization by Cordero-Grande et al. [7] , which is one of very few approaches for retrospective motion estimation in 3D MRI. The method alternates between two steps of L1-minimization reconstruction with wavelet regularization while fixing the motion parameters and four steps of motion parameter estimation while fixing the reconstruction. For the final reconstruction we perform L1-minimization from scratch based on the estimated motion parameters (AltOpt-L1) and with additional DC loss thresholding based on the estimated motion parameters (AltOpt+Th-L1), to ensure a fair comparison to our method. We also perform L1-minimization without any motion estimation (L1). We further compare to the E2E stacked U-net with self-assisted priors [1] , a method for end-to-end (E2E) motion artifact reduction in 3D MRI. The model from Al-Masni et al. [1] is trained with pairs of motion-free images and undersampled images corrupted by motion simulated as described above to predict the motion-corrected image with a single pass trough the network at inference.Hyperparameters for MotionTTT and the baselines are in Appendix D.MotionTTT accurately estimates motion over a wide range of motion severities. The results in Figure 4 show the reconstruction performance in PSNR as a function of motion severity averaged over the test set and over two independently sampled motion trajectories per example. Figure 5 and Appendix E.1.1 contain example reconstructions. For both small and large motion severities, the PSNRs for reconstruction with known and estimated motion parameters is the same for MotionTTT, indicating that the motion parameters are estimated very well. Reconstruction results for the motion parameters itself are in Appendix E.1.1.For AltOpt this is only true for small motion severities, for large ones MotionTTT significantly outperforms AltOpt. In addition, MotionTTT is about 6x faster than AltOpt for this problem, see Appendix E.1.2.The performance of E2E stacked U-net lies significantly below the performance of methods that obtain a reconstruction based on explicitly estimating motion like MotionTTT or AltOpt. With increasing motion severity E2E methods suffer from the severe degradations in their network inputs, which can not be resolved without knowledge of the motion resulting in overly smooth reconstructions.For the results in Figure 4 we used the fixed undersampling mask from Figure 1 (c) with acceleration factor 4. In Appendix E.1.3 we additionally ablate over acceleration factors 2 and 8 showing robust performance across acceleration factors and levels of motion severity.Reconstruction module. Figure 4 shows that for no motion, U-net based reconstruction (KnownMotion-U-net-DCLayer) outperforms L1-based reconstruction (KnownMotion-L1), as expected, but for larger levels of motion severity L1-based reconstruction performs slightly better. A likely reason for this is a distribution shift, which is known to result in worse performance [8] : Acquiring MRI measurements under motion typically leads to some parts of the k-space being sampled more than once, while parts that would be sampled if there were no motion are not sampled. This results in a change of effective undersampling mask and a motion-specific increase of the effective undersampling factor, which suggests why the reconstruction quality decays for all methods with increasing levels of motion even if motion parameters are known. This problem, illustrated in Figure 2 , is well known [45, 12] . However, L1-based reconstruction is not as sensitive to changes in the mask as the U-net based reconstruction, which is trained with a fixed mask. Due to this distribution shift, the performance of U-net might degrade stronger than for L1-based reconstruction.Finally, we note that, comparing MotionTTT-L1 and MotionTTT+Th-L1 Figure 4 shows that for strong motion, DC loss thresholding essentially closes the gap to KnownMotion-L1. This indicates that the thresholding does not unnecessarily exclude motion states for moderate motion, but reliably detects incorrectly estimated motion states for severe motion.Next, we study the performance of our method for intra-shot motion. As mentioned before, the model, data, and baselines are the same as for the inter-shot experiments in Section 6.1.Motion simulation. We again set the number of shots to B = 50. We use a random sampling trajectory (Figure 1 (d)) , where the 3 × 3 center is acquired first and all other k-space lines are acquired at a random order, and again assume the subject to be in zero-motion state at the first shot. We randomly select ⌈N e /2⌉ of the motion events to take place during the acquisition of one shot. Intra-shot motion is simulated by assigning a distinct motion state to each of the 182 k-space lines acquired during this shot such that the intra-shot motion trajectory connects the motion parameters from the previous to the next shot, where start and end point as well as the presence of up to two peaks due to over-and/or under-shooting within the intra-shot trajectory is randomized.Choice of the parameter N splits . As explained in Section 4, MotionTTT estimates intra-shot motion by splitting motion states that exhibit a large DC loss after phase 1 of the optimization scheme into N splits motion states. Since ground truth intra-shot motion is simulated with a distinct motion state for each of the 182 k-space lines acquired during one shot, choosing the parameter N splits < 182 results in a discretization error. On the other hand, choosing N splits large increases the difficulty of the estimation problem because the number of k-space lines corresponding to one motion state, decreases. We found N splits = 10, i.e., about 18 k-space lines per motion state, to be a good trade-off, see Appendix E.2.1 for an ablation study.Sampling order. The sampling order within a shot is crucial for the ability to estimate intra-shot motion as estimating a motion state corresponding to a batch consisting of only high-frequency and low-energy components is difficult. In Appendix E.2.2 we ablate over different sampling orders and find that acquiring all k-space lines at a random order works well.MotionTTT with intra-shot motion estimation improves performance over discarding measurements corrupted by intra-shot motion. Figure 6 shows the reconstruction performance as a function of motion severity, where half of the motion events exhibit intra-shot motion. Recall that after phase 1 of the optimization scheme outlined in Section 4 MotionTTT converged to one estimated motion state per shot. MotionTTT-L1 (Phase 1) is based on those motion states and hence does not estimate intra-shot motion. Consequently, even for the lowest level of motion considered here, MotionTTT+Th-L1 (Phase 1) improves the performance as shots corrupted by intra-shot motion are discarded during DC loss thresholding before reconstruction.In contrast, MotionTTT-L1 (Phase 3) achieves on par performance with MotionTTT+Th-L1 (Phase 3) for moderate motion levels indicating that intra-shot motion has been estimated successfully so that motion states are below the threshold we set for DC loss thresholding. Nevertheless, a performance gap remains relative to KnownMotion-L1, which results from the irreducible discretization error.For severe motion the gap between MotionTTT-L1 (Phase 3) and MotionTTT+Th-L1 (Phase 3) increases as more motion states are estimated incorrectly and have to be discarded. However, MotionTTT with only Phase 1 continues to be outperformed indicating the benefit of performing intra-shot motion estimation over discarding measurements corrupted by intra-shot motion. Figure 6 : Reconstruction performance in PSNR as a function of the simulated motion severity defined as in Figure 4 only that here half of the motion events exhibit intra-shot motion. We consider L1minimization based on known motion, motion estimated with MotionTTT after phase 1 (no intra-shot motion estimation) and after phase 3 (with intra-shot motion estimation) of the optimization scheme.Error bars are the standard deviation over test examples and randomly sampled motion trajectories.Figure 11 in Appendix E contains a qualitative comparison, where especially in the case without DC loss thresholding the difference in reconstruction quality between phase 1 and 3 is clearly visible.After thresholding the differences are less visible, however that depends on the amount of intra-shot motion that we simulate as with more intra-shot motion the difference in undersampling factor between discarding and estimating intra-shot motion becomes larger and differences more visible.We now apply MotionTTT to real motion data. We use pre-trained model f θ from Section 6.1, and the implementation details are described in Appendix D.Data. We acquired four scans from one subject. To obtain motion-free reference data, mildly motion-corrupted, and strongly motion-corrupted data the subject was instructed to not move at all or move 1-3 times at distinct time points during the acquisition. The performed motions include nodding, head rotations, and either with or without returning to the original position. Sequence parameters (see Appendix D.4) were chosen to match those in the Calgary Campinas Brain MRI Dataset [37] as close as possible. Data was acquired with an undersampling factor of 4.94 and a random sampling trajectory (similar to Figure 1 (d) ) with B = 52 shots. The acquisition of one shot lasts 1.3s followed by a pause of 1.6s resulting in a total scan duration of about 150s.Results. We find that MotionTTT achieves significantly improved visual reconstruction quality compared to no motion-correction. Applying MotionTTT results in a significant reduction of motion artifacts for both mild and strong motion (reconstructed images are in Appendix E.3). 
limitations and future work: 
 In this paper, we proposed the first deep-learning-based 3D rigid motion estimation method for 3D MRI and have demonstrated that it is effective and computationally managable at estimating motion and correcting for it.As discussed in Section 6.1 our reconstruction module (Step 3) has room for improvement. We currently use L1-minimization, but a deep-learning method, e.g., regularization with diffusion models should yield further improvements.Finally, the model used to perform MotionTTT can be improved in principle. Currently it is trained on fully-sampled data, which is scarce especially for 3D at high resolutions. It might be possible to train this module well with self-supervised training losses that require only undersampled data [44, 28, 20] . Moreover, as mentioned, there is a distribution shift in the mask and networks that work well with such distribution shifts might yield improvements.As mentioned in the Problem Statement Section 3, we simulate motion in the k-space rather than in the image domain because it is computationally more efficient if the number of motion states b is larger than the number of coils C. Here, we elaborate on this statement.In our work the MRI forward model under motion is implemented with the NUFFT, which for each shot in the sampling trajectory first computes the rotated coordinates based on the k-space data and the motion parameters of this shot. Then the k-space values at those coordinates for all shots can be obtained from a single application of the NUFFT. This requires us to compute a single interpolated version of the k-space data, which consists of C-many 3D coil volumes, thus C-many interpolated volumes need to be computed.In contrast, simulating motion in the image domain required computing a transformed 3D image volume for each motion state, which then is expanded to the coil dimension and transformed to the k-space with the forward model ( 1). Hence, b many interpolated volumes need to be computed. As in our work the number of coils C = 12 is smaller than the number of motion states b it is computationaly more efficient to simulate motion in the k-space than in the image space.In this appendix, we prove Theorem 1 from the theory Section 5. To prove the result, we upper bound the loss for the correct motion parameters, L(m * ) and lower bound the loss L(m) for all other motion parameters m ̸ = m * .Assume without loss of generality that the ground-truth shift is equal to m * = 0. We have thatL(m) = ∥D m F T f (F * T D * m y) -y∥ 2 2 = D m F T n bk UU T F * T D * m D m * F T x -D m * F T x 2 2 i = n bk D m F T UU T F * T D * m F T x -F T x 2 2 ii = n bk F T UU T F * D * m F T x -D * m F T x 2 2 = n bk F T UU T F * T D * m F T Uc -D * m F T Uc 2 2, where D * is the Hermitian transpose of the matrix D. Here, equation i follows from the assumption that the optimal motion parameters are zero, thus D m * = I, and equation ii follows from the entries of D m having absolute value one, and D m D * m = I. We first upper bound L(m * ) = L(0). We have thatEQUATIONFor the last inequality, we used that d/bk ≤ 1/4 by assumption. According to Theorem 2.6 in Rudelson and Vershynin [33] , the second to last inequality holds on the eventsEQUATIONandEQUATIONwith β = √ n. Those events hold with the probabilities, for all β > 0EQUATIONEQUATION)Next, we lower-bound L(m) for m ̸ = 0. Let a be the number of individual motion parameters in the vector m that are non-equal to the true motion parameters m * = 0. We haveL(m) = F T U n bk U T F * T D * m F T U - a b I c + F T U a b -D * m F T U c 2 2 ≥ F T U a b -D * m F T U c 2 -F T U n bk U T F * T D * m F T U - a b I c 2 2 ≥ F T U a b -D * m F T U c 2 -∥F T U∥ n bk U T F * T D * m F T U - a b I c 2 2 ≥ F T U a b -D * m F T U c 2 2 -∥F T U∥ 2 n bk U T F * T D * m F T U - a b I c 2 2 ≥ |1 -a/b| 2 bk n 1 2 - 7 8 bk n • (1 + α) β 2 n n 2 k 2 b 2 (4b + d) = |1 -a/b| 2 bk n 1 2 - 7 8 • (1 + α) β 2 n n kb (4b + d) ,where the last inequality holds with probability at least 1 -e -cα + 2d(d + b)e -cβ 2 -3e -cd which follows from the boundsEQUATIONandEQUATIONThus, we have with probability at leastEQUATIONwhich is equivalent toEQUATION)By a union bound over all motion parameter m ∈ {0, . . . , n -1} \ {m * } (there are n b many), we have thatEQUATIONwhich concludes the proof. It remains to prove the intermediate results, which we do next.Since the entries of D * m have absolute value one, and a/b ∈ [0, 1], we haveF T U a b -D * m F T U c 2 ≥ ∥(1 -a/b)F T Uc∥ 2 ≥ |1 -a/b|σ min (F T U)∥c∥ 2 ≥ |1 -a/b| bk n 1 - d bk 1 - β √ n ≥ |1 -a/b| bk n 7 2 8 2 ,where the second to last inequality holds with probability at least 1-3e -cd according to equations ( 14) and (15) , and where we used that d bk ≤ 1/64 and d n ≤ 1 64 .B.2 Lower-bounding the second term, proof of equation ( 16)Define A = n bk U T F * T D * m F T U -a b I for notational convenience. We have that P ∥Ac∥ 2 2 ≥ (1 + α) β 2 n n 2 k 2 b 2 (4b + d) ≤ P ∥A∥ 2 F ≥ d β 2 n n 2 k 2 b 2 (4b + d) + P ∥Ac∥ 2 2 ≥ (1 + α) 1 d ∥A∥ 2 F ≤ P ∥A∥ 2 F ≥ d β 2 n n 2 k 2 b 2 (4b + d) + P ∥Ac∥ 2 2 - 1 d ∥A∥ F ≥ α 1 d ∥A∥ F ≤ e -cα + 2d(d + b)e -cβ 2where the last inequality follows from the Hanson-Wright inequality as well as fromEQUATION)We next prove the bound (21) . We start with upper bounding the Frobenius norm of the matrix A. We split the squared Frobenius norm into a sum of the squared diagonal entries and squared off-diagonal entries according to∥A∥ 2 F = n kb U T F * T D * m F T U - a b I 2 F = d i=1 n kb u * i F * T D * m F T u i - a b 2 + d i=1 d j̸ =i n kb u * i F * T D * m F T u j 2 .It follows thatP ∥A∥ 2 F ≥ d β 2 n n 2 k 2 b 2 (4b + d) ≤ P ∥A∥ 2 F ≥ d β √ n 1 + n k √ b 2 + d(d -1) n kb 2 β 2 n ≤ d i=1 P n kb u * i F * T D * m F T u i - a b 2 ≥ β √ n 1 + n k √ b 2 + d i=1 d j̸ =i P n kb u * i F * T D * m F T u j 2 ≥ n kb 2 β 2 n ≤ d2(1 + b)e -cβ 2 + d 2 e -cβ 2 .The last inequality follows fromEQUATION)andEQUATION)Bounding a diagonal entry, proof of inequality (22) : First, consider a diagonal element and let Z ℓ ∈ R n×n be the mask that selects the frequencies in the set T ℓ , and recall that F ∈ R n×n is the Fourier transform. With a slight abuse of notation, we let D m be the diagonal matrix that contains the frequencies that matrix the Fourier matrix it is multiplied with, i.e., in D m F, the matrix D m is the n × n diagonal matrix with entries e i2πmℓ/n , ℓ = 0, . . . , n -1, and in D m F T , the matrix D m is the |T | × |T | diagonal matrix with entries e i2πmℓ/n , ℓ ∈ T . For convenience, we drop the index i and write u = u i . With this we haveEQUATIONHere, the entries of ũ = Fu ∈ R n are iid CN (0, 1/ √ n) distributed, since the DFT matrix F has orthonormal columns. Thus, by a union boundEQUATION)where we used that, for all β > 0,EQUATIONandEQUATIONThis concludes the proof of inequality (22) . It remains to proof equations ( 26) and (27) .Proof of equation (26): Note thatEQUATIONSince z b,i -k n is a sub-Gaussian zero mean random variable a concentration inequality for sub-Gaussians yieldsEQUATION)Here, we used that ũ * • ũ is the entrywise product, and the random variable ∥ũ * • ũ∥2 2 = n i=1 (ũ * i ũi ) 2 concentrates around its expectation n3σ 4 = n3/n 2 = 3/n.Thus we have shown that the random variables s ℓ = ũ * Z ℓ -k n I D * b ũ (conditioned on ũ) are sub-Gaussian. The random variables are also zero-mean and independent (if conditioned on ũ), and thus by concentration of sub-Gaussian random variables we getP b ℓ=1 ũ * Z ℓ - k n I D * m ℓ ũ ≥ 1 √ n β √ b ≤ 2e -cβ 2 .This concludes the proof of equation (26) .Proof of equation ( 27): Next, consider the random variable u * F * D * m ℓ Fu = ũ * D * m ℓ ũ in equation (24) . If m ℓ = 0, i.e., there is no shift, we have D 0 = I, and thus the random variable becomes ũ * ũ which is a sum of sub-exponential random variables and concentrates around 1. In case m ℓ is an integer non-equal to zero, we have that u * F * D * m ℓ Fu = u T u m ℓ , where u m ℓ is a vector circularly shifted by m ℓ . We can write u T u m ℓ as two sums of independent Gaussians, to see this consider the case m ℓ = 1 and note that for this case we haveu T u m ℓ = (u 1 u 2 + u 3 u 4 + . . .) S1 + (u 2 u 3 + u 4 u 5 + . . .) S2 .By the union bound and by Hoeffding's inequality, we get thatEQUATION)Combining this via a union bound for all summands ℓ = 1, . . . , b yields the bound in equation (27) .Bounding an off-diagonal element, proof of inequality (23) : Now consider an off-diagonal element ũ * i D * m ũj . Since ũi and ũj are independent and the entries have zero mean and variance 1/n, this is a sum of Gaussian random variables with variance 1/n. Thus, by sub-Gaussian concentration or Hoeffding's inequality, we have thatEQUATION)where we used that ∥D * m ũj ∥ In the main body, we stated Equation ( 7), i.e.,EQUATIONwhere T = T 1 ∪ . . . ∪ T b is the set of all measurements collected.To see that this approximation is accurate, note that for the noiseless case with a known shift, where y = D m Uc the network approximately reconstructs the signal sinceEQUATIONwhere the first equality follows because the matrix D m F T has orthonormal rows, and the approximation holds since U T F * T F T U concentrates around n kb if U is a random subspace, and if the number of measurements, T is sufficiently large relative to the dimension of the subspace, d.To perform motion parameter estimation as outlined in Section 4 we use gradient based optimization with the Adam [19] optimizer. Phase 1 runs for 70 iterations with an initial learning rate of 4.0, which is decayed twice by a factor of 4 at iterations 40 and 60. The DC loss threshold used to determine incorrectly estimated motion states after phase 1 is set to 0.575. If no motion states fall above the threshold, the optimization continues for another 30 iterations with one additional learning rate decay after 10 iterations. If motion states fall above the threshold, phase 2 runs for 30 iterations with a learning rate of 0.5. Finally, phase 3 runs for another 30 iterations with a leraning rate of 0.05.For severe motion we found that optimizing only over rotation parameters for the first few steps facilitates their correct estimation. To avoid single motion parameters to get stuck at a large value early during the optimization we clamp the estimated motion parameters at [5, 8, 10, 12, 15] degrees/millimeters for steps smaller than [15, 30, 45, 60, 150] .We use TorchKbNufft (MIT license), an implementation of the NUFFT from Muckley et al. [29] and use the option for density compensation based on the method of Pipe [31] when applying the adjoint NUFFT. To allow differentiation with respect to the input coordinates we build on the extension Bindings-NUFFT-pytorch from Alban Gossard (MIT license).We run MotionTTT on a Nvidia L40 GPU with 46GB memory or on a Nvidia A100 GPU with 80GB memory. Within our implementation there are three hyperparameters that control the required GPU memory. Recall that in every iteration of MotionTTT the entire 3D volume is reconstructed slice-wise. However, gradients are only computed for a subset of randomly selected slices of size 5 as described in Section 4. The size of this subset is the first hyperparameter that we can control.The second parameter is the batch size of the NUFFT. As described in Appendix A the NUFFT is applied for each of the C coils. The TorchKbNufft package allows batch wise computation at the cost of increased GPU memory utilization.As in our 3D setup we have to deal with a lot more motion states (B = 50) than previous work that studied the 2D setup, we implemented a third option to reduce GPU memory consumption. To this end, we split the estimated motion states into two batches of size of, e.g., 25 each and backpropagate the gradients subsequently before performing a single optimizer step that updates the estimated motion states.Note that hyperparameters batch size of NUFFT and of motion states per backpropagation only affect the run time but do not change the optimization problem, where the cost of the latter could be compensated by increasing the number of GPUs accordingly. Hence, at the cost of prolonged run times or a larger GPU cluster MotionTTT can be applied with any number of motion states. With our hardware and for our C = 12-many coils we could use a NUFFT batch size of 12 (4) and a batch of 50 (25) motion states for backpropagation in case of the A100 (L40) GPU.We perform DC loss thresholding for excluding the k-space data acquired during the i-th shot from the reconstruction based on its estimated motion state mi and the DC loss (5) L TTT ( mi ) > δ with a threshold of δ = 0.575.We perform L1-minimization with wavelet regularization based on the estimated motion parameters. We run 50 steps with SGD and a learning rate of 5 × 10 7 and regularization weight λ = 10 -3 .The U-Net reconstruction given the true motion parameters KnownMotion-U-net-DCLayer from Figure 4 , is obtained by fine-tuning a DCLayer as proposed in Chen et al. [4] . If the U-Net reconstruction is x U-Net = f θ A † (T , -m)y , then the reconstruction of the DCLayer is:x = arg min x ∥A(T , m)x -y∥ 1 ∥y∥ 1 + λ ∥x -x U-Net ∥ 1 ∥x U-Net ∥ 1In the experiments, we set λ = 0.1. The choice of learning rate and number of steps is critical for optimizing the DCLayer. After a grid searching, we identified the optimal learning rate and steps for various severity levels, which are summarized in the following For applying MotionTTT to the real motion-corrupted data in Section 6.3 we use the same hyperparameter configuration as described above up to choosing a slightly smaller initial learning rate of 1.0 during phase 1 to explore the space of motion parameters more slowly and set the threshold parameter to 0.70 due to a generally higher DC loss level of the scanner data compared to the data used for the simulation.For obtaining the final reconstruction with L1-minimization based on the motion parameters estimated with MotionTTT we set the hyperparameters learning rate to 10 -8 and regularization weight to λ = 3 × 10 -8 .To perform alternating optimization as described in Section 6.1 we run SGD with a learning rate of 5 × 10 7 and regularization weight λ = 10 -4 during the reconstruction steps and a learning rate of 5 × 10 -11 during the motion estimation step. In both steps the loss is the MSE between predicted and given measurement. The optimization process is capped at 500 iterations, but it terminates early if the reconstruction loss falls below the threshold of e 13 .After alternating optimization we perform L1-minimization from scratch based on the estimated motion parameters. We run 50 steps with SGD and a learning rate of 5 × 10 7 and regularization weight λ = 10 -3 .We train the stacked U-net with self-assisted priors proposed in Al-Masni et al. [1] as an example from the class of methods for end-to-end motion correction for 3D MRI. Their method consists of a stack of two U-nets that are trained to map the absolute value of a zero-filled (ZF) reconstruction of the motion corrupted and undersampled measurement to the absolute value of the motion-free image enabling end-to-end motion correction with a single forward pass through the network. The network operates slice-wise however neighboring slices are provided as context at the network input to account for the 3D nature of the problem.We train the E2E stacked U-net on the same dataset and in a similar way as we trained the U-net on motion-free data in Section D.1.1. In every step one of the 40 3D training volumes is loaded from which the fully sampled target volume is computed as well as three different input volumesthe motion-free undersampled ZF reconstruction and two different motion-corrupted undersampling ZF reconstructions generated with the inter-shot motion simulation described in Section 6.1 with a random level of motion severity and a fixed motion severity of level 1 respectively. For each training example and level of motion severity a fixed set of three different motion trajectories is considered during training.We include the task of motion-free and motion severity of level 1 reconstruction in every training step to ensure that the model does not loose too much performance when reconstructing motion-free data compared to a model trained only for motion-free reconstruction and to bias the model towards performing well for very mild motion artifacts as it is the regime in which end-to-end methods can still provide reasonably good results, whereas for severe motion the reconstruction results are far from a potential use for medical diagnosis. See Figure 7 for examples of mild and very severe motion.For each of the three degraded volumes 20 random slices are selected in each plane r x × r y , r We train the model for 230 epochs with the L1 loss from equation ( 37) and the Adam optimizer with a learning rate of 3 × 10 -4 which is decayed twice by a factor of 10 at epochs 200 and 220 and pick the model from the epoch with the highest reconstruction PSNR for motion-free reconstruction on the validation set.We adopt the network design from Al-Masni et al. [1] , where we set the number of channels in first layer of both U-nets to 64 resulting in a total of 15.9M network parameters. We used instance norm instead of batch norm in our network as we found it to give more stable results across levels of motion severity.The final motion-corrected reconstructions are obtained by reconstructing the 3D volume slice-wise in the axial r x × r y plane. Before the PSNR scores shown in Figure 4 are computed the reconstruction is aligned with the motion-free reference volume.In this section we provide additional information regarding the acquisition of our own data used in Section 6.3. This study was exempt from Institutional Review Board (IRB), but potential risks were disclosed to the subject and experiments were conducted with informed consent. Data was acquired on a Ingenia Elition 3.0T X scanner (Philips Healthcare, Best, The Netherlands) using the standard 16-channel dStream HeadSpine coil array, where C = 13 channels were used during acquisition. We perform 3D T1-weighted Ultra-fast Gradient-echo (TFE) imaging with a 1mm isotropic resolution and a matrix-size of k x × k y × k z = 222 × 236 × 512, an undersampling factor of 4.94 and a linear sampling trajectory illustrated in Figure 1 (d) . We subsample the data by a factor of two along the fully-sampled frequency encoding dimension k z to obtain a similar field of view as in our training data with k z = 256. In one shot 204 lines in the k-space are acquired resulting in a total number of 52 shots. See Table 1 for an overview of all sequence parameters. In this section we provide additional ablation studies and further qualitative examples to complement the experimental results presented in Sections 6.1, 6.2, and 6.3 on simulated inter-/intra-shot motion and real-motion respectively. We present additional results on inter-shot motion estimation, analysing reconstructions at different levels of motion severity. Qualitative comparisons for mild and severe cases highlight the strengths and limitations of the motion estimation methods used.Figure 5 in Section 6.1 in the main body shows reconstructed images for a simulated motion severity level 9.Figure 7 shows additional reconstruction results at a lower severity level (level 3); for this lower motion severity level both MotionTTT+Th-L1 and AltOpt+Th-L1 achieve results comparable to KnownMotion-L1.Figure 8 shows an example of predicted motion and corresponding DC loss for a simulated inter-shot motion scenario at severity level 9. The DC loss effectively detects incorrectly estimated motion states, highlighting their locations. This capability is particularly useful for DC thresholding, which improves robustness, as discussed in Section 6.1. In this scenario, MotionTTT failed in only one shot, whereas the AltOpt method failed in 12 shots, further demonstrating the effectiveness of the MotionTTT approach.Beyond reconstruction performance, computational time is an important factor for real-world applications. If the reconstruction process is too slow, the algorithm may be impractical for clinical use. The following table summarizes the average computational times for the AltOpt and MotionTTT methods, tested on an Nvidia A100 GPU. As shown in the table, our MotionTTT method is approximately 6 times faster than AltOpt when early stopping is applied. Without early stopping, AltOpt is even 10 times slower than MotionTTT. In this section we investigate the role of the undersampling factor on the ability of MotionTTT to estimate inter-shot motion. We re-train the U-net on two additional Cartesian undersampling masks with acceleration factors R = 2, 8 in addition to the existing results with acceleration factor R = 4 (see Figure 9 a,b,c for the masks). Figure 10 shows the reconstruction performance in PSNR based on motion parameters estimated by our MotionTTT compared to ground truth motion over three levels of motion severity and the three acceleration factors.As expected, the overall performance decays with increasing acceleration factors and motion severities.For mild and moderate motion, MotionTTT achieves highly accurate motion estimation for all acceleration factors indicated by the vanishing performance gap relative to using ground truth motion.For the most severe motion, a small performance gap exists for all acceleration factors due to incorrectly estimated motion states that are discarded from the final reconstruction via DC loss thresholding. In fact, under severe motion an average of 2.5/100, 2.0/50 and 0.12/25 shots have to be discarded for acceleration factors 2,4 and 8.We conclude that MotionTTT can achive highly accurate motion parameter estimation robustly across different acceleration factors. We attribute the slight increase in discarded motion states for smaller acceleration factors to the increased complexity of the optimization problem as the number of unknown motion states to be estimated increases linearly in the number of acquired shots.In this section we present addition ablation studies for the choice of the number of motion states N splits estimated per shot and the sampling order used in our experimental results on simulated intra-shot motion estimation in Section 6.2. We also show reconstruction results from the experiments in the main body in Figure 11 .We start with an ablation study on the choice of the hyperparameter N splits that determines the number of motion states that are introduced at the end of phase 1 of MotionTTT's optimization scheme outlined in Section 4 for each shot that exhibits a data consistency loss larger than a certain threshold. As discussed in Section 6.2, the choice of this hyperparameter trades off the irreducible discretization error versus the available signal per estimated motion state and the computational complexity.The discretization error results from estimating only N splits -many motion states for a shot that when affected by intra-shot motion exhibits a distinct motion state for each k-space line (182 in our setup) acquired within this shot. The error decreases with increasing N splits and depends on severity of motion as fast movements with a large amplitude result into a large discretization error.On the other hand, increasing N splits decreases the amount of k-space signal per estimated motion state potentially leading to more incorrectly estimated motion states. Additionally, doubling the number of motion states N splits doubles the computational complexity of phase 2 if the available hardware does not allow for parallelization (see Appendix D.1.3 for a discussion of computational aspects).To find the value for N splits that trades off those two effects, we conduct the following experiment. We simulate motion trajectories containing 5 intra-shot motion events following Section 6.2 with maximal motion M max = 5. We assume that all motion parameters are known except the ones during the intra-shot events. We then apply MotionTTT+Th-L1 with a random sampling order to estimate the motion states during intra-shot motion for different levels of discretization defined by the number of motion states per shot N splits ∈ {5, 10, 20}. This simulates the case where phase 1 of the intra-shot motion estimation scheme perfectly estimates the inter-shot motion parameters and detects all shots corrupted by intra-shot motion for which then motion parameters are estimated during phase 2. This allows us to focus the evaluation solely on the ability of the method to perform intra-shot motion estimation independently from the performance in phase 1.We obtain PSNRs averaged over the four examples in the validation set and five simulated motion trajectories per example as 35.84, 36.05 and 36.07 for N splits = 5, 10 and 20. Further, the average numbers of k-space lines that are discarded due to DC loss thresholding before the final reconstruction are 289, 197 and 302 out of 5 * 182 = 910 lines affected by intra-shot motion.For N splits = 5 we obtain the lowest performance, which we attribute to a large discretization error, which even results into more motion states with a DC loss above the threshold than for N splits = 10.For N splits = 10, 20 we obtain the same performance, however due to different reasons. While for N splits = 20 a smaller discretization error can be achieved, more motion states are estimated incorrectly and are discarded due to a more difficult optimization problem as more motion states are estimated with less signal per motion state. We use N splits = 10 for our experiments in Section 6.2 as it corresponds to smaller computational costs for the same performance.Next, we present an ablation study on the role of the sampling order for the ability of MotionTTT to estimate intra-shot motion as in Section 6.2.In order to focus our evaluation on the performance for intra-shot motion estimation, we consider the same experimental setup as in Appendix E.2.1, where we simulate motion trajectories containing 5 intra-shot events with maximal motion M max = 5, apply MotionTTT only to the motion states affected by intra-shot motion and assume all other motion states to be known. We fix the number of motion states per shot to N splits = 10 and report the performance for different sampling orders including the interleaved and random sampling orders used for the inter-and intra-shot motion experiments, as well as a linear and a deterministic order. See Figure 9 masks c)-e) for visualizations.The linear order acquires k-space lines according to their k y index from low to high. Hence, all shots towards the start and end of the scan only contain high-frequency components. The interleaved order distributes high-/low-frequency components evenly across shots, but within one shot the order follows the k x index from low to high. Hence, all k-space lines acquired towards the start and end of one shot are high-frequency components. The random order acquires all k-space lines at a random order. The deterministic order we designed such that the average distance to the center of the k-space of two k-space lines acquired after each other is approximately constant, which ensures that after the acquisition of a high-frequency components a low-frequency component is acquired. The random, interleaved and deterministic orders all start by sampling the 3 × 3 center of the k-space.Table 2 shows the average performance of MotionTTT+Th-L1 in PSNR for each sampling order as well as the average number of k-space lines that are discarded due to the DC loss thresholding before the final reconstruction.While random and deterministic orders perform equally, the interleaved order results in more incorrectly estimated motion states. For the linear order we obtain a significant loss in performance, and over a third of the 9078 lines acquired in total are discarded including many lines that are not affected by intra-shot motion and for which we assume knowledge of the true motion state in this experiments. In fact, under the linear order entire shots at the start and end of acquisition are discarded due to their large DC loss despite assuming knowledge of the motion states. We attribute this to the U-net used during MotionTTT reconstructing high-frequency components not as faithfully as low-frequency components, and conclude that in order to reliably estimate the motion state of a batch of measurements the batch needs to contain low-frequency components.The deterministic order fulfills this requirement by design and the random order with high probability as for the number of motion states per shot N splits = 10 the smallest batch of measurements we consider consists of about 18 k-space lines, which is large enough to contain both high-and lowfrequency components especially as the center of the k-space is much more densely sampled than towards the edge. In our experiments we hence use a random sampling scheme due to its ease of implementation on top of any arbitrary design of the undersampling mask.We now discuss the real motion experiments from Section 6.3 in more detail.Visual reconstruction analysis. Figure 12 shows the reconstruction of the volume in three orthogonal planes (axial, sagittal, and coronal). For mild motion, artifacts are still noticeable, particularly in the k y -k z plane. However, the MotionTTT-L1 method effectively mitigates most of these artifacts, resulting in a much clearer reconstruction than the other methods. For strong motion, many details across all three planes are entirely occluded, yet the MotionTTT-L1 method significantly improves image quality, restoring visibility to numerous details that would otherwise be lost due to severe motion.Estimated motion parameters and DC loss. Figure 13 shows the estimated motion parameters as well as the DC loss (5) per state for the strong motion example. The number of motion states is 69, which results from starting phase 1 of MotionTTT with 52 motion states (one per shot), of which four exceeded the threshold after phase 1. As the number of additional motion states per corrupted shot is N splits = 10, phase 2 and 3 continue with 92 motion states from which 23 were discarded during the final DC loss thresholding resulting in the 69 motion states depicted in Figure 13 . The areas shaded in gray correspond to the motion states that were added after phase 1, i.e., intra-shot motion states, whereas the white areas correspond to motion states each pertaining a single shot. Consequently, the predicted motion and the DC loss of MotionTTT (phase 1) is constant within the gray areas as intra-shot motion estimation is performed starting from phase 2. As we can see intra-shot estimation, i.e., MotionTTT (phase 3) significantly reduced the DC loss compared to phase 1 and the estimated motion is reasonably smooth, indicating the benefit of estimating intra-shot motion in practice. Justification: The abstracts claims to demonstrate the capabilities of our new method with simulated and real motion data. Experimental results for both are described in Section 6.Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [Yes] Justification: The proposed method is targeted at the specific problem of correcting motion artifacts in MRI, which potentially improves the quality of medical diagnosis and increases the effective patient throughput of the MRI machine as discussed in Section 1. We can not see any negative societal impacts arising from the proposed approach. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [Yes] Justification: The license of the used dataset [37] is mentioned in Section 6.1. Specialized software packages used are cited with license in Appendix D. Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 13. New Assets Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets? Answer: [Yes] Justification: We release the code for our method at https://github.com/MLI-lab/ MRI_MotionTTT. The method itself is described in detail in Section 4 and an instruction on how to use the code can be found in the repository. Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)? Answer: [Yes] Justification: One subject participated in the real motions experiments; a compensation was not provided due to the negligible effort and risk of harm; details about instructions given to the participant are provided in manuscript Section 6.3. Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained? 
 
 Appendix:  


Decision: The manuscript is on the verge of rejection, yet the authors have presented a compelling rebuttal. This AC acknowledges the complexity of the problem and appreciate this first work towards 3D motion-compensated MRI reconstruction. This is why the AC believes this could be an interesting paper to be discussed at the conference. While this AC is inclined to support weak acceptance, the AC also understands it may not meet the priority threshold given the overall acceptance criteria/ranking. If the paper is accepted, it is with the expectation that the key points from the rebuttal will be effectively integrated into the manuscript to enhance its overall impact.

Reviews: review: paper_summary: paper_summary: The paper presents MotionTTT, a deep learning-based method for estimating and correcting rigid motion in 3D MRI images. The approach leverages a neural network pre-trained for 2D motion-free image reconstruction and employs test-time-training (TTT) to estimate motion parameters from motion-corrupted 3D measurements. The effectiveness of the method is demonstrated through evaluations on both simulated and real datasets.
strength_weakness: strength_weakness: The theoretical foundations of MotionTTT are generally robust, but there are notable gaps. The method builds on existing neural networks for 2D image reconstruction, extending them to handle 3D rigid motion estimation. 

During the pre-training step, $f_\theta$ is trained to map under-sampled k-space data back to fully sampled k-space data. The authors show that under-sampled k-space data with motion (where certain regions of the k-space undergo rotation and phase shifting) results in higher reconstruction loss for $f_\theta$. However, this indirect optimization lacks theoretical proof and needs more robust experimental validation. 

For instance, the performance of this method may heavily depend on the under-sampling pattern of the k-space (such as the undersampling factor and the linear sampling trajectory) and the specific regions in k-space that are affected by motion. However, in the experiments, the setup of the undersampling of the k-space was fixed. 

A more thorough theoretical and experimental investigation into why this optimization works and under what conditions it is most effective would strengthen the soundness of the method.
questions: questions: 1. The idea behind minimizing loss (4) is not intuitive and clear. How do you prove the effectiveness of using loss (4) to recover and remove the motion corruption? Will the effectiveness be affected by the undersampling factor and trajectory? What if you have a higher or lower factor? or maybe a spiral trajectory? Can you find the optimal factor for your method? 
2. Based on the assumption that rigid motion leads to rotation and phase shift in k-space, is the linear trajectory the optimal undersampling trajectory? 
3. Do you really need the fully sampled k-space data to learn $f_\theta$?
limitations: limitations: The indirect TTT optimization method lacks theoretical justification. The authors should investigate both theoretically and experimentally why this approach works and under what conditions it is most effective.

The performance of MotionTTT may be heavily dependent on the under-sampling pattern of the k-space and the specific regions affected by motion. The method’s effectiveness might vary significantly with different sampling patterns and motion artifacts.

The method relies on a pre-trained network using fully-sampled data, which may not always be available in practice.
paper_summary: paper_summary: This paper proposes a method for estimating the motion of patients, such that accurate motion-corrected images could be reconstructed. The key idea is that a neural network trained for motion-free reconstruction has a small loss if there is no motion, thus optimizing over motion parameters passed through the reconstruction network enables estimation of the motion.
strength_weakness: strength_weakness: - There is no technical novelty in this work. It introduces some existing modules to address the proposed problem.
- The experiments are weak to demonstrate the efficacy of the proposed method. 
  1）There lacks strong baselines and there are no compared methods. Although according to the authors, there are no other methods designed for 3D rigid motion estimation for 3D motion-corrected MRI, there should be some methods that could reconstruct the images.
  2) The data are little for evaluating the methods, and more averaged quantitative results are suggested to report.
questions: questions: See the weaknesses.
limitations: limitations: Yes, limitations are discussed.
paper_summary: paper_summary: This paper proposes a motion correction MRI reconstruction algorithm for 3D brain MRI. The proposed technique consists in a deep learning-based estimation of rigid motion parameters, which allows to correct the k-space before a final reconstruction.
Estimation of motion parameters are based on a single optimization step freezing the reconstruction UNet parameters, and assuming loss will only (or mainly) be influenced by motion.
strength_weakness: strength_weakness: The overall structure of the paper is quite confusing, which makes it hard to follow. Experiments and alternative techniques seem to be presented throughout the result section.
A proper ablation study seems required, in order to assess the added-value of each step. Proper baseline also seems necessary especially regarding the reconstruction step or the removal of motion artifacted lines, as more recent techniques could have been used as SOTA.
The overall motion estimation relies on a first 2D reconstruction pipeline however the performance of the technique seems to be lower than a standard L1 reconstruction technique. If so, how could the authors be sure that motion parameters are not biased or affected by the low performance of the 2D reconstruction pipeline.
Why did the authors not try to directly reconstruct from the 3D k-space domain, in order to estimate motion parameters in 3D? Moreover, I believe the authors should also provide further details on they are switching back and forth from 3D k-space domains to 2D, as a simple slicing in a given direction of the 3D k-space is not appropriate.
questions: questions: Could the authors find another database of 3D data in order to train a 3D reconstruction UNet directly?
Given the proposed technique aims also at suppressing corrupted k-spaces lines, would it be sensible to train a reconstruction technique with higher a acceleration factor?
Would the proposed technique be applicable for other sampling strategies, especially would it be possible to apply on stack of stars acquisitions or other 3D sampling strategy.
Why did the authors not apply other reconstruction (standard L1) reconstruction for the real dataset for comparison? Especially since it was outperforming the UNet approach on simulated motion data.
limitations: limitations: The authors should discuss the risks of hallucinations and bias of the proposed technique.
paper_summary: paper_summary: This work proposes a novel deep learning-based framework for recovering high-quality 3D MR images from undersampled and motion-corrupted k-data. The proposed approach is well motivated and technically sound. The authors perform extensive experiments on simulated and real MR datasets, which confirm the effectiveness of the proposed framework.
strength_weakness: strength_weakness: For this submission, I have a few minor suggestions as follows.

**Various types of rigid motion**

In line 205, the authors show that random rigid motions ($M_\text{max}=[2,5,10]$) are simulated. However, the rigid motion in the real world could follow some patterns, such as involuntary motion and abrupt motion. I think it is better to test these different movements.

**Compared methods**

For supervised methods, only U-net is used as a baseline. Advanced supervised models [1][2] are not discussed. Furthermore, motion correction methods [3] based on diffusion models are not compared. 

> [1] Han Y, Yoo J, Kim H H, et al. Deep learning with domain adaptation for accelerated projection‐reconstruction MR[J]. Magnetic resonance in medicine, 2018, 80(3): 1189-1205.

> [2] Liu J, Kocak M, Supanich M, et al. Motion artifacts reduction in brain MRI by means of a deep residual network with densely connected multi-resolution blocks (DRN-DCMB)[J]. Magnetic resonance imaging, 2020, 71: 69-79.

> [3] B. Levac, A. Jalal, and J. I. Tamir. “Accelerated Motion Correction for MRI Using Score-Based 388 Generative Models”. ISBI 2023.

**Some typos**

For example, the symbol $\mathcal{E}$ in Eq. 3 is not defined.
questions: questions: See the section of Strengths and Weaknesses, please.
limitations: limitations: The authors mention enough limitations of their work.

Comments: title: Discussion period is coming to an end
comment: Can you be specific which concerns remain? You point out that “this paper is a good trial to address this problem” but at the same time argue that the work lacks novelty and that there “should be other methods”, without pointing to any specific work in the literature. If you have any concrete concerns that remain, we would be happy to clarify.
Thanks for noting that our paper makes a significant contribution to 3D motion correction, and many thanks for the reply.

Regarding comparing to more methods, we would have been happy to compare to more methods, but found no suitable 3D methods other than the one we compare to (alternating minimization). Please not that all methods you mentioned are for 2D motion correction, and they do not extend to 3D in a straightforward manner. We think it can’t be expected to extend other work significantly in order to generate new baselines, this would be independent contributions. Motion estimation and correction in 3D is extremely difficult, this is testified by the fact that the vast majority or papers only studies 2D motion reconstruction, even though motion in MRI always occurs in 3D space. Please let us know if you have any baseline in mind that has been applied to 3D motion estimation or correction by the papers that proposed the method.
Thanks to the authors for the rebuttals. Some of the concerns are addressed but some remain. After reading the comments of all reviewers, I keep my original rating.
Thank you for the detailed response. I have also reviewed the comments from other reviewers. While I believe this paper makes a significant contribution to 3D motion correction in MRI, there are some limitations, such as the limited comparison of methods. Therefore, I recommend a weak acceptance for this submission.
We would like to thank the reviewers once again for their insightful comments.
As the open discussion period is coming to an end, we would very much appreciate if the reviewers could kindly respond to our rebuttal, this would provide us with the opportunity to address or clarify any remaining concerns.