naranjal white

# Motion corrected MRI with DISORDER: Distributed and Incoherent Sample Orders
for Reconstruction Deblurring using Encoding Redundancy

L Cordero-Grande  Centre for the Developing Brain, School of Biomedical
Engineering and Imaging Sciences, King’s College London, King’s Health
Partners, St Thomas’ Hospital, London, SE1 7EH, UK  Biomedical Engineering
Department, School of Biomedical Engineering and Imaging Sciences, King’s
College London, King’s Health Partners, St Thomas’ Hospital, London, SE1 7EH,
UK  G Ferrazzi  Biomedical Engineering Department, School of Biomedical
Engineering and Imaging Sciences, King’s College London, King’s Health
Partners, St Thomas’ Hospital, London, SE1 7EH, UK  RPAG Teixeira  Centre for
the Developing Brain, School of Biomedical Engineering and Imaging Sciences,
King’s College London, King’s Health Partners, St Thomas’ Hospital, London,
SE1 7EH, UK  Biomedical Engineering Department, School of Biomedical
Engineering and Imaging Sciences, King’s College London, King’s Health
Partners, St Thomas’ Hospital, London, SE1 7EH, UK
J O’Muircheartaigh  Centre for the Developing Brain, School of Biomedical
Engineering and Imaging Sciences, King’s College London, King’s Health
Partners, St Thomas’ Hospital, London, SE1 7EH, UK  Department of Forensic and
Neurodevelopmental Sciences, Institute of Psychiatry, Psychology and
Neuroscience, King’s College London, King’s Health Partners, 16 De Crespigny
Park, London, SE5 8AF, UK  AN Price  Centre for the Developing Brain, School
of Biomedical Engineering and Imaging Sciences, King’s College London, King’s
Health Partners, St Thomas’ Hospital, London, SE1 7EH, UK  Biomedical
Engineering Department, School of Biomedical Engineering and Imaging Sciences,
King’s College London, King’s Health Partners, St Thomas’ Hospital, London,
SE1 7EH, UK  JV Hajnal  Centre for the Developing Brain, School of Biomedical
Engineering and Imaging Sciences, King’s College London, King’s Health
Partners, St Thomas’ Hospital, London, SE1 7EH, UK  Biomedical Engineering
Department, School of Biomedical Engineering and Imaging Sciences, King’s
College London, King’s Health Partners, St Thomas’ Hospital, London, SE1 7EH,
UK

###### Abstract

Purpose: To enable rigid-body motion tolerant parallel volumetric magnetic
resonance imaging by retrospective head motion correction on a variety of
spatio-temporal scales and imaging sequences.

Theory and methods: Tolerance against rigid-body motion is based on
distributed and incoherent sampling orders for boosting a joint retrospective
motion estimation and reconstruction framework. Motion resilience stems from
the encoding redundancy in the data, as generally provided by the coil array.
Hence, it does not require external sensors, navigators or training data, so
the methodology is readily applicable to sequences using 3D encodings.

Results: Simulations are performed showing full inter-shot corrections for
usual levels of in-vivo motion, large number of shots, standard levels of
noise and moderate acceleration factors. Feasibility of inter- and intra-shot
corrections is shown under controlled motion in-vivo. Practical efficacy is
illustrated by high quality results in most corrupted of $208$ volumes from a
series of $26$ clinical pediatric examinations collected using standard
protocols.

Conclusion: The proposed framework addresses the rigid motion problem in
volumetric anatomical brain scans with sufficient encoding redundancy which
has enabled reliable pediatric examinations without sedation.

Index terms— magnetic resonance, image reconstruction, motion correction,
parallel imaging, distributed and incoherent sampling

##  1 Introduction

Tolerance against motion is desirable in magnetic resonance imaging (MRI).
This includes brain MRI, where significant motion-induced image degradation
prevalence has been documented (Andre et al., , 2015) and high resolution
imaging quality may be compromised by head motion (Budde et al., , 2014).
Rigid-body MRI motion correction (Zaitsev et al., , 2015; Godenschweger et
al., , 2016) can be tackled via prospective or retrospective techniques.
Prospective techniques (Maclaren et al., , 2013) compare advantageously in
terms of spin-history, dephasing confounders or k-space density guarantees.
Particularly, optical tracking systems have been proposed for head motion
estimation, with corrections showing impressive accuracy and latency (Schulz
et al., , 2012). However, prospective methods require additional hardware
and/or scanner modifications, often involving intrusive markers attached to
the subject. In addition, satisfactory corrections may not always be possible
due to unpredictability or complexity of motion, or maximized sampling
efficiency requirements. Retrospective techniques may facilitate scanning or
improve prospective results (Aksoy et al., , 2012), particularly for 3D
encodings, where spin-history is less of a problem, and when using non-linear
reconstruction paradigms (Adcock et al., , 2014) to deal with non-homogeneous
sampling density after motion.

Motion compensation is strongly dependent on motion estimation from the
measured information. Some methods have proposed the use of navigators, where
surrogate motion-sensitive information is interleaved with the main
acquisition and correction is applied either prospective or retrospectively
(Tisdall et al., , 2016; Johnson et al., , 2016; Gallichan and Marques, ,
2017). Due to variability in time requirements for different MRI sequences,
application of a given navigator is usually limited to a specific sequence
type. Furthermore, particular care has to be taken to prevent spin-history or
saturation and, sometimes, scanning efficiency may be compromised.
Alternatively, sequences can be constructed with relative resilience to motion
or, similarly, sampling schemes can be designed to function as implicit
navigators. This is the case for spiral and radial trajectories (Bammer et
al., , 2007; Anderson III et al., , 2013; Pipe et al., , 2014), where
temporally distributed low resolution information is used for motion
estimation, with retrospective corrections usually grounded on an intermediate
reconstruction of fully formed images for each motion state, often involving
non-linear methods. Finally, other approaches have explored the redundancy of
the information sensed by parallel MRI to detect and discard localized
inconsistencies in k-space measurements (Bydder et al., , 2002), usually
requiring prior image models to limit noise amplification and improve
inconsistency detection (Samsonov et al., , 2010).

Building on models of MRI acquisition in the presence of motion (Batchelor et
al., , 2005; Bammer et al., , 2007), some methods have proposed formulations
for motion estimation from the k-space that do not require navigators (Odille
et al., , 2008; Loktyushin et al., , 2013; Cordero-Grande et al., , 2016). Our
previous work (Cordero-Grande et al., , 2016) introduced a data-driven
reconstruction method for retrospective multi-shot rigid-body motion
correction or _aligned reconstruction_ taking advantage of the encoding
redundancy in the measured data. Performed simulations showed that the ability
to solve the aligned reconstruction problem is strongly sensitive to the
k-space encoding order, which suggested that opportunities exist to maximize
the sensitivity to motion by appropriate sampling order designs. Consequently,
in this paper we introduce the Distributed and Incoherent Sample Orders for
Reconstruction Deblurring using Encoding Redundancy (DISORDER) framework as a
flexible way to correct for head motion on a variety of spatio-temporal scales
and imaging contrasts by optimizing the _sample orders_ for k-space coverage.
In addition, we propose some technical refinements to the aligned
reconstruction formulation and extend the simulation domain. The technique is
implemented on a $3\,\mbox{T}$ scanner and tested on controlled motion scans
and pediatric examinations including magnetization-prepared rapid acquisition
gradient echo (MP-RAGE), fast spin echo (FSE), fluid attenuated inversion
recovery (FLAIR), spoiled gradient echo (SPGR), and balanced steady-state free
precession (bSSFP) sequences. A MATLAB implementation to reproduce the
experiments is made available at
https://github.com/mriphysics/DISORDER/releases/tag/1.1.0.

##  2 Theory

###  2.1 Aligned reconstruction

Assuming whitened measurement noise (Pruessmann et al., , 2001), the aligned
reconstruction for parallel volumetric imaging can be formulated as:

$(\hat{\mathbf{x}},\hat{\boldsymbol{\theta}})=\displaystyle\operatornamewithlimits{argmin}_{\mathbf{x},\boldsymbol{\theta}}r_{\mathbf{x},{\boldsymbol{\theta}}}=\displaystyle\operatornamewithlimits{argmin}_{\mathbf{x},\boldsymbol{\theta}}\|\mathbf{A}\boldsymbol{\mathcal{F}}\mathbf{S}\mathbf{T}_{\boldsymbol{\theta}}\mathbf{x}-\mathbf{y}\|_{2}^{2},$
(1)

where $\mathbf{x}$ is the image to be reconstructed, $\boldsymbol{\theta}$ are
the motion parameters, $r$ is the loss function, $\mathbf{y}$ is the measured
k-space data, $\mathbf{T}$ is a set of rigid motion transformations,
$\mathbf{S}$ are the coil sensitivities, $\boldsymbol{\mathcal{F}}$ is the
discrete Fourier transform (DFT), and $\mathbf{A}$ is a sampling mask. We are
interested in reconstructing a 3D image of size $V=V_{1}V_{2}V_{3}$ with
$V_{d}$ the number of voxels along dimension $d$ from $\displaystyle
N=C\sum_{m=1}^{M}E_{m}$ $C$-element coil array samples of a discretized
k-space grid of size $K$. $E_{m}$ denotes the number of samples within
_segment_ $m$ and $M$ is the number of segments in the sequence, with each
segment associated to a specific motion state. Detailed information about the
terms in Eq. (1) can be found in Cordero-Grande et al.,  (2016). Here we
provide a brief description of their structure:

  * •

$\mathbf{y}$ is a $N\times 1$ vector.

  * •

$\mathbf{A}$ is a $N\times KMC$ block matrix comprising submatrices of size
$E_{m}\times K$ whose entries take the value $1$ if the sample $e$ of the
segment $m$ corresponds to the k-space location indexed by $k$ and $0$
otherwise.

  * •

$\boldsymbol{\mathcal{F}}$ is a $KMC\times VMC$ block diagonal matrix
comprising submatrices of size $K\times V$ representing 3D DFT’s with applied
k-space sampling.

  * •

$\mathbf{S}$ is a $VMC\times VM$ block matrix comprising diagonal submatrices
of size $V\times V$ whose diagonal elements correspond to the spatial
sensitivity of the coil $c$.

  * •

$\mathbf{T}$ is a $VM\times V$ block matrix comprising unitary (Unser et al.,
, 1995) submatrices of size $V\times V$ corresponding to the 3D rigid
transformation modeling the motion state $m$ by three translations and three
Euler rotation angles codified in the parameter vector
$\boldsymbol{\theta}_{m}$.

  * •

$\mathbf{x}$ is a $V\times 1$ vector.

Eq. (1) is a separable nonlinear least squares problem (Gan et al., , 2018;
Herring et al., , 2018). We confront it by iteratively addressing the
subproblems:

$\begin{split}\hat{\mathbf{x}}^{(i+1)}=&\displaystyle\operatornamewithlimits{argmin}_{\mathbf{x}}\|\mathbf{A}\boldsymbol{\mathcal{F}}\mathbf{S}\mathbf{T}_{\hat{\boldsymbol{\theta}}^{(i)}}\mathbf{x}-\mathbf{y}\|_{2}^{2}\\\
\hat{\boldsymbol{\theta}}^{(i+1)}=&\displaystyle\operatornamewithlimits{argmin}_{\boldsymbol{\theta}}\|\mathbf{A}\boldsymbol{\mathcal{F}}\mathbf{S}\mathbf{T}_{\boldsymbol{\theta}}\hat{\mathbf{x}}^{(i+1)}-\mathbf{y}\|_{2}^{2}.\end{split}$
(2)

The first subproblem, reconstructing the image $\mathbf{x}$ in the presence of
rigid motion (Batchelor et al., , 2005), can be solved by conjugate gradient
(CG) (Pruessmann et al., , 2001). As for the second, the solution must null
the gradient of the objective function against the motion parameters (Cordero-
Grande et al., , 2016), which is tackled by a Levenberg-Marquardt (LM)
algorithm using a simplified Jacobian (Ruano et al., , 1991). A natural
initialization is a zero-motion condition
$\hat{\boldsymbol{\theta}}^{(0)}=\mathbf{0}$, so the first step corresponds to
a standard sensitivity encoding (SENSE) reconstruction. Further in this paper,
we describe how to temporally arrange the k-space samples into segments to
improve the aligned reconstruction convergence.

###  2.2 DISORDER sampling

We focus on Cartesian 3D k-space grids with uniform sampling as sketched in
Fig. 1. Fig. 1a shows $K_{1}=4$ collected samples after the first readout or
profile in the $k_{1}$ direction. Fig. 1b shows the first segment, in this
example corresponding to the full acquisition of the $k_{1}$-$k_{2}$ plane.
Fig. 1c shows that segments can be used to define an ordered partition of the
$k_{1}$-$k_{2}$-$k_{3}$ grid. Due to short duration, we assume negligible
motion during the readout and focus on the phase encode (PE) plane in Fig. 1d.
We define $E^{\text{PE}}_{m}$ as the number of profiles per segment, so
$E^{\text{PE}}_{m}=E_{m}/K_{1}$, and hereinafter we adopt the replacement
$E_{m}\leftarrow E^{\text{PE}}_{m}$.

\begin{overpic}[width=108.18665pt,draft=false]{Fig01/fig01}\put(0.0,0.0){{a)}}\end{overpic}

\begin{overpic}[width=108.18665pt,draft=false]{Fig01/fig02}\put(0.0,0.0){{b)}}\end{overpic}

\begin{overpic}[width=108.18665pt,draft=false]{Fig01/fig03}\put(0.0,0.0){{c)}}\end{overpic}

\begin{overpic}[width=108.18665pt,draft=false]{Fig01/fig04}\put(0.0,0.0){{d)}}\end{overpic}

Figure 1: Sketch of volumetric sampling for an exemplary $K=4\times 4\times 4$
space (in radians). a) Measured samples after acquiring the first readout. b)
Measured samples after acquiring the first segment. c) Measured samples after
whole sequence acquisition with color coding used to differentiate each
segment. d) View of the $k_{2}$-$k_{3}$ PE plane.

By modifying the PE gradients before each readout it is possible to design
different encoding or view orders. These can be defined as a temporally
ordered set of profiles
$p\in\mathcal{P}=\\{(k^{1,1}_{2},k^{1,1}_{3}),\ldots,(k^{1,E_{1}}_{2},k^{1,E_{1}}_{3}),\ldots,(k^{M,E_{M}}_{2},k^{M,E_{M}}_{3})\\}$
with cardinality $\displaystyle P=|\mathcal{P}|=\sum_{m}E_{m}$. Fig. 2a shows
the first segment of a commonly used Sequential ordering scheme. In this case,
due to the partition definition, a segment includes two consecutive $k_{3}$
planes. Fig. 2b, introduces the Checkered traversal. First, a rectangular
tiling of the PE plane is built using tiles of size $U_{2}\times U_{3}$ such
that $U_{2}U_{3}=M$. Second, a spectral lexicographic order for the profiles
within a tile $\mathcal{K}_{U}$ is defined by
$\mathcal{K}_{U}\to\mathcal{M}_{U}=\\{1,\ldots,M\\}$, which can be extended to
different tiles by translation. Third, interleaved segments are defined such
that the same profile $m_{u}\in\mathcal{M}_{U}$ is used $\forall
e\in\mathcal{E}=\\{1,\ldots,E\\}$, with $\mathcal{E}$ a temporally ordered set
of tiles. Finally, the profile sequence to traverse each tile is defined by
mapping from the set of profiles within a tile to a temporally ordered set of
segments $\mathcal{M}_{U}\to\mathcal{M}_{T}=\\{1,\ldots,M\\}$ using an
electrostatic repulsion criterion with periodic boundary conditions. Hence, a
_distributed_ temporal coverage is guaranteed both for the whole spectrum and
within each tile. This strategy aids the aligned reconstruction conditioning
by reducing the chances for large uncovered spectral areas due to head
rotation. Fig. 2c presents the Random-checkered modification. Tiles are built
as in the Checkered approach but segments are constructed by a random
permutation of the tile elements drawn independently for each tile, so we have
$m_{u_{e}}$. This guarantees a distributed coverage in probability and
introduces some _incoherence_ among profiles within and across segments. Fig.
2d shows a Random view order where $\mathcal{P}$ is a random permutation of
the profiles conforming the sequence. Note that, considering a free definition
of segments, Sequential and Random schemes are particular cases of the Random-
checkered traversal respectively with tiling sizes of $1\times 1$ and
$K_{2}\times K_{3}$.

\begin{overpic}[width=108.18665pt,draft=false]{Fig02/fig01}\put(8.0,0.0){{a)}}\end{overpic}

\begin{overpic}[width=108.18665pt,draft=false]{Fig02/fig02}\put(8.0,0.0){{b)}}\end{overpic}

\begin{overpic}[width=108.18665pt,draft=false]{Fig02/fig03}\put(8.0,0.0){{c)}}\end{overpic}

\begin{overpic}[width=108.18665pt,draft=false]{Fig02/fig04}\put(8.0,0.0){{d)}}\end{overpic}

Figure 2: Example of segments for different encoding orders to cover a
$P=K_{2}K_{3}$ PE plane (in radians) with $K_{2}=K_{3}=32$ using an
acquisition partitioned into $M=16$ segments with equal number of profiles per
segment $E_{s}=E=P/M=64$. Full set of profiles to cover the PE plane as red
circles, profiles within a segment filled in blue, underlying tilings
($U_{2}=U_{3}=4$) in purple, samples at half the spatial resolution enclosed
in cyan, and areas covered by $4$ intra-segment temporal subdivisions of
considered segment enclosed in yellow. a) Sequential; b) Checkered; c) Random-
checkered; and d) Random traversals.

View orders should preserve the contrast and the trajectory consistency. We
establish the following differentiation:

  * •

Non steady-state sequences (MP-RAGE, FSE, FLAIR). They acquire a fraction of
the k-space or _shot_ after each radiofrequency (RF) preparation. Thus, they
induce a natural sampling partition where segments are in correspondence to
shots. In addition, magnetic properties are not invariant for the different
shot samples. Typically, middle samples within each shot cover the central
area of the spectrum (Busse et al., , 2008), which our orders can fulfill by
jumping from tile to tile in a Zig-zag manner. An example is shown in Fig. 3
using an elliptical sampling area. The Checkered traversal produces regular
segment patterns (first column), with neighboring colors maximally separated
within the tile, whereas the Random-checkered traversal produces non-regular
patterns. Tiling orders generate smooth color transitions across the spectrum
for all presented traversals (second column), which translates into smooth
magnetic properties of the profiles.

  * •

Steady-state sequences (SPGR, bSSFP). They produce a temporally stable
magnetization after reaching the steady-state, usually facilitated by some
preparatory dummy profiles, so the contrast becomes independent of the
encoding order. For estimates attempted at the segment level, temporal
resolvability of motion will increase with bigger tiling sizes $U_{2}U_{3}$.
However, large jumps in the spectrum may induce inconsistencies due to eddy
currents, especially for low repetition times, so a trade-off may be required.
If the profiles are covered by the application of $M$ spectral _sweeps_ ,
analogously to Tsao et al.,  (2005), eddy currents can be minimized by an
Alternating zig-zag tiling order where the traversal polarity is reversed for
consecutive sweeps. This is illustrated in Fig. 4. The segment structure
(first column) matches that of Fig. 3 but for some minor differences in the
Sequential case due to smooth magnetization requirements in shot-based
sequences (Busse et al., , 2008). The Sequential scheme guarantees a smooth
passage through k-space (third and fourth columns). Although quicker k-space
sweeps of our traversals imply larger $\differential k_{2}$ and $\differential
k_{3}$ steps, these remain substantially lower than for the Random order,
which should limit the impact of eddy currents. Finally, the Alternating zig-
zag suppresses the undesirable spikes in the fourth column of Fig. 3.

\begin{overpic}[width=106.67146pt,draft=false]{Fig03/fig01-01}\put(0.0,0.0){{a)}}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig03/fig01-02}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig03/fig01-03}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig03/fig01-04}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig03/fig02-01}\put(0.0,0.0){{b)}}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig03/fig02-02}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig03/fig02-03}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig03/fig02-04}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig03/fig03-01}\put(0.0,0.0){{c)}}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig03/fig03-02}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig03/fig03-03}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig03/fig03-04}\end{overpic}

Figure 3: Segment $m$ and tiling $e$ temporal orders together with trajectory
derivatives $\text{d}k_{2}$ and $\text{d}k_{3}$ (left to right) for Zig-zag
tiling orders used in non steady-state sequences. The example corresponds to a
case with $P=3630$ profiles, $M=30$ segments and tiling pattern $U_{2}\times
U_{3}=6\times 5$. a) Sequential; b) Checkered; and c) Random-checkered
segments.

\begin{overpic}[width=106.67146pt,draft=false]{Fig04/fig01-01}\put(0.0,0.0){{a)}}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig04/fig01-02}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig04/fig01-03}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig04/fig01-04}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig04/fig02-01}\put(0.0,0.0){{b)}}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig04/fig02-02}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig04/fig02-03}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig04/fig02-04}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig04/fig03-01}\put(0.0,0.0){{c)}}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig04/fig03-02}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig04/fig03-03}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig04/fig03-04}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig04/fig04-01}\put(0.0,0.0){{d)}}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig04/fig04-02}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig04/fig04-03}\end{overpic}

\begin{overpic}[width=106.67146pt,draft=false]{Fig04/fig04-04}\end{overpic}

Figure 4: Segment $m$ and tiling $e$ temporal orders together with trajectory
derivatives $\text{d}k_{2}$ and $\text{d}k_{3}$ (left to right) for
Alternating zig-zag tiling orders used in steady-state sequences. The example
corresponds to a case with $P=3630$ profiles, $M=30$ segments and tiling
pattern $U_{2}\times U_{3}=6\times 5$. a) Sequential; b) Checkered; c) Random-
checkered; and d) Random segments.

###  2.3 Aligned reconstruction refinements

We propose a series of refinements for improved and more efficient aligned
reconstructions:

  * •

Spatial multiresolution. The spatial and spectral grids for both subproblems
in Eq. (2) can be refined according to a given multiresolution pyramid as
commonly used in image registration (Unser et al., , 1993). In contrast to
sequential sampling, the proposed orders allow for motion estimates from
samples at coarse scales (area enclosed in cyan in Fig. 2) to be completely
exploited when reconstructing at fine scales. This is useful for quick aligned
reconstructions as adequate motion estimates are often possible at coarse
scales.

  * •

Temporal multiresolution. Motion estimation can also be attempted at intra-
shot (or intra-sweep) levels (for instance using the samples enclosed within
the yellow areas for the $4$ intra-segment subdivisions in Fig. 2). However,
estimates using low spatial harmonics localize the structures at coarse scales
only and, conversely, motion estimation using high spatial harmonics alone is
limited by lower SNR and prone to local optima. These limitations can be
alleviated using hierarchical estimation refinements by temporally subdividing
the samples considered by the motion states within a shot in a coarse-to-fine
manner.

  * •

Coil compression. The two subproblems in Eq. (2) can operate on a reduced
number of virtual channels (Buehrer et al., , 2007).

  * •

Motion compression. The reconstruction subproblem complexity grows with the
number of motion states, which can be reduced by motion compression or
binning. Estimated motion parameter traces are approximated using piecewise
constant functions by Haar wavelet decomposition truncation with threshold
$\boldsymbol{\tau}$. Thereby, original motion states are packed into effective
states by grouping together those contiguous states with similar motion
parameters into an effective motion parameter vector
$\tilde{\boldsymbol{\theta}}$. Thus, the reconstruction complexity is driven
by the underlying motion complexity.

  * •

Robustness. Accurate intra-shot corrections may be infeasible, for instance
due to temporary inconsistencies in the magnetization. Denoting the real
motion parameters by $\boldsymbol{\theta}^{\ast}$ we can ideally characterize
the loss $r_{\mathbf{x},\boldsymbol{\theta}^{\ast}}$ using the sampling noise
properties. Sampling noise follows a circularly symmetric complex Gaussian
additive stationary distribution and, after whitening, it is independent
across channels, so the losses per profile
$r[m,e_{m}]=\sum_{c,k_{1}}r[m,e_{m},c,k_{1}]$ should ideally follow a
$\chi^{2}$ distribution. To account for the sensitivity of the residuals to
the underlying signal, we use trimmed statistics on a logarithmic scale
$r_{b}[m]=P^{\mathbb{E}_{m}}_{c(b)}\log(r[m,e_{m}])$ with $b\in{1,\ldots,B}$
indexing the $100c(b)\,\%$ centile $P_{c}$ of the loss distribution through
k-space $\mathbb{E}_{m}$. As we are concerned with anomalously high residuals,
robust estimates of the scale and mean of the statistic distribution across
segments $\mathbb{M}$ are obtained respectively by
$\sigma_{b}=\sqrt{2}(P^{\mathbb{M}}_{c_{\text{U}}}r_{b}[m]-P^{\mathbb{M}}_{c_{\text{L}}}r_{b}[m])/(\operatornamewithlimits{erfc}^{-1}(2c_{\text{U}})-\operatornamewithlimits{erfc}^{-1}(2c_{\text{L}}))$
and
$\mu_{b}=P^{\mathbb{M}}_{(c_{\text{U}}+c_{\text{L}})/2}r_{b}[m]+\sqrt{2}\sigma_{b}\operatornamewithlimits{erfc}^{-1}(c_{\text{U}}+c_{\text{L}})$
choosing $c_{\text{U}}=0.25$ and $c_{\text{L}}=0.125$. Using these estimates,
the statistics are normalized and averaged into
$\overline{r}[m]=\sum_{b}(r_{b}[m]-\mu_{b})/B\sigma_{b}$, and segments are
weighted in the reconstruction by a matrix $\mathbf{W}$ with entries
$w[m]=\min(M\operatornamewithlimits{erfc}(\overline{r}[m]/\sqrt{2})/(2\tau_{w}),1)$,
with $\tau_{w}$ an acceptance threshold corrected for multiple comparisons.

  * •

Regularization. If outlier segment rejection is activated or the
reconstruction is applied to accelerated scans, some form of regularization
may be advisable. This is considered by reformulating the reconstruction as:

$\hat{\mathbf{x}}^{(i+1)}=\displaystyle\operatornamewithlimits{argmin}_{\mathbf{x}}\|\mathbf{W}^{1/2}(\mathbf{A}\boldsymbol{\mathcal{F}}\mathbf{S}\mathbf{T}_{\tilde{\boldsymbol{\theta}}^{(i)}}\mathbf{x}-\mathbf{y})\|_{2}^{2}+2\lambda\|\boldsymbol{\mathcal{S}}\mathbf{x}\|_{1},$
(3)

where $\lambda$ controls the degree of regularization and
$\boldsymbol{\mathcal{S}}$ corresponds to a shearlet decomposition, which
provides nearly optimal approximation rates for piecewise smooth functions
with discontinuities on a piecewise smooth surface (Kutyniok et al., , 2012).
We resort to an iteratively reweighted least squares (IRWLS) solver, able to
produce high quality solutions in a few iterations (Voronin and Daubechies, ,
2017), with $\lambda$ adaptively updated according to Li and Zhong,  (2018)
using a normalized Rayleigh-quotient trace estimator (Avron and Toledo, ,
2011).

##  3 Methods

###  3.1 Synthetic experiments

Our contributions are validated using a synthetic dataset built from a $T_{2}$
neonatal brain axial ground truth (GT) image without perceptible motion
artifacts. This corresponds to a multi-slice TSE sequence acquired on a
$3\,\mbox{T}$ Philips Achieva TX (same scanner as for in-vivo tests) using a
$C=32$-element neonatal head coil array, $0.8\times 0.8\,\mbox{mm}$ in-plane
resolution, $1.6\,\mbox{mm}$ slice thickness, echo time
$T_{\text{E}}=145\,\mbox{ms}$, repetition time $T_{\text{R}}=12\,\mbox{s}$,
and flip angle $\alpha=90^{\circ}$. Coil sensitivities were estimated from a
separate reference scan (Allison et al., , 2013). We use a 2D dataset and no
regularization or outlier rejection for a concise presentation of results. We
assume that the simulated 2D k-space corresponds to the $k_{2}$-$k_{3}$ PE
plane of 3D scans and expect the driving conclusions to be extensible to 3D
because estimates should be easier along the missing fully-sampled readout
direction.

Simulations were conducted to compare the conventional sequential order to the
various proposed schemes as well as to characterize their performance. The
forward model in the presence of rigid motion is applied to the GT to generate
synthetically motion corrupted data. Synthesized measures are corrupted with
noise levels corresponding to a mean SNR of $30\,$dB for reconstructions in
the absence of motion or acceleration. Different degrees of motion are
generated by drawing independent motion states uniformly at random on an
interval of rotations $[-\theta/2,\theta/2]$ around the field of view (FOV)
center. Satisfactory convergence in the presence of noise can be ascertained
on the assumption of an identifiable global optimum basin by
$r_{\hat{\mathbf{x}},\hat{\boldsymbol{\theta}}}\leq
r_{\hat{\mathbf{x}},\boldsymbol{\theta}^{\ast}}$. In this case, the error in
the motion parameters $\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{\ast}$
is attributed to the uncertainty from the measurement noise and not to partial
convergence. Note that we can generally achieve a lower loss for the joint
problem ($r_{\hat{\mathbf{x}},\hat{\boldsymbol{\theta}}}$) than with the
knowledge of the motion parameters
($r_{\hat{\mathbf{x}},\boldsymbol{\theta}^{\ast}}$) due to the larger
complexity of the former. Reconstructions are terminated when
$r_{\hat{\mathbf{x}},\hat{\boldsymbol{\theta}}}\leq
r_{\hat{\mathbf{x}},\boldsymbol{\theta}^{\ast}}$ and the abscissa scale of the
convergence plots was chosen so that iterations have a direct translation into
computational costs.

###  3.2 In-vivo experiments

In-vivo experiments include the main families of volumetric sequences for
brain MRI (see Table 1). We have performed a controlled motion experiment on a
consented adult volunteer and applied the method to replace sedation on
pediatric subjects scanned after written informed parental consent for an
epilepsy study. Imaging is performed using a $32$-channel adult head coil.
Data is acquired in the inferior-superior (IS) $k_{1}$, anterior-posterior
(AP) $k_{2}$ and left-right (LR) $k_{3}$ orientation using our scanner
implementation of the Random-checkered traversal. This way, potentially
strongest rotations on the sagittal plane are captured by the $k_{1}$-$k_{2}$
coordinates, which may increase the resolvability of intra-shot motion. In
addition, IS readouts allow to easily downweight additional motion sources
within the FOV when estimating for motion by restricting the loss function to
the superior part of the FOV ($2/3$ in our implementation). Finally, this
orientation facilitates non-selective RF excitation pulses for shorter
$T_{\text{R}}$.

Subjects | Modality | $T_{\text{R}}$ | $T_{\text{E}}$ | $T_{\text{I}}$ | $\alpha$ | $P$ | $U_{2}\times U_{3}$ | $\Delta_{\mathbf{y}}$ | $R$ | $\text{FOV}_{\mathbf{x}}$ | $t$
---|---|---|---|---|---|---|---|---|---|---|---
Adult | MP-RAGE | $8.3\,\mbox{ms}$ | $3.7\,\mbox{ms}$ | $1.00\,\mbox{s}$ | $8.0^{\circ}$ | $3630$ | $6\times 5$ | $1.5\,\mbox{mm}$ | $2\times 2$ | $240\times 240\times 220\,\mbox{mm}$ | $1\,\mbox{min}\,29\,\mbox{s}$
Pediatric | MP-RAGE | $7.0\,\mbox{ms}$ | $2.3\,\mbox{ms}$ | $0.90\,\mbox{s}$ | $8.0^{\circ}$ | $17920$ | $8\times 10$ | $1.0\,\mbox{mm}$ | $1.4\times 1.4$ | $240\times 240\times 188\,\mbox{mm}$ | $3\,\mbox{min}\,07\,\mbox{s}$
TSE | $2.50\,\mbox{s}$ | $0.31\,\mbox{s}$ | – | $90.0^{\circ}$ | $17955$ | $9\times 15$ | $5\,\mbox{min}\,42\,\mbox{s}$
FLAIR | $5.00\,\mbox{s}$ | $0.41\,\mbox{s}$ | $1.80\,\mbox{s}$ | $18200$ | $10\times 10$ | $8\,\mbox{min}\,30\,\mbox{s}$
SPGR | $8.2\,\mbox{ms}$ | $3.0\,\mbox{ms}$ | – | $12.3^{\circ}$ | $21654$ | $1.3\times 1.3$ | $3\,\mbox{min}\,06\,\mbox{s}$
bSSFP | $6.0\,\mbox{ms}$ | $3.0\,\mbox{ms}$ | $46.0^{\circ}$ | $2\,\mbox{min}\,16\,\mbox{s}$

Table 1: Sequence parameters for the different modalities considered in the
experiments. $T_{\text{I}}$ is the inversion time, $\Delta_{\mathbf{y}}$ the
acquired resolution, $R$ the uniform acceleration factor,
$\text{FOV}_{\mathbf{x}}$ the reconstructed FOV, and $t$ the scan duration.

In the controlled motion experiment, the volunteer was first scanned without
deliberate motion, and then asked to perform extreme and continuous motion for
the entire scan, which was repeated three times. Reconstructed volumes are
jointly registered together for error comparisons. The pediatric cohort
includes $26$ subjects ranging from $3$ to $19$ years old (mean$\pm$std of
$12\pm 5$ years), typically acquiring one MP-RAGE, TSE and FLAIR, two SPGRs
and three bSSFPs for an approximate total of $208$ tested volumes across all
participants. Strongest artifacts in our data are generally arising from
motion, so the reported case has been separately chosen for each modality as
the most artifacted after reconstruction without motion correction.

###  3.3 Implementation details

In the in-vivo experiments sensitivities are compressed into a number of
channels corresponding to a $10\%$ SNR loss. The number of resolution levels
is defined as
$L=\displaystyle\left\lfloor\log_{2}(4\,\mbox{mm}/\Delta_{\mathbf{y}}^{\text{min}})\right\rfloor+1$,
with $\lfloor\cdot\rfloor$ denoting the biggest integer lower or equal than
the argument and $\Delta_{\mathbf{y}}^{\text{min}}$ the minimum of the voxel
sizes along different directions. As we use $2\times$ subsampling ratios, we
operate at a minimum resolution of $4\,\mbox{mm}$. In the first iteration at
level $l$, a soft-masked (Fuderer et al., , 2004) full CG reconstruction is
run till loss reduction saturation. Then, the method quickly alternates
between reconstruction and motion correction using one CG and one LM iteration
with heuristically updated damping and line search. We activate a flag for
provisional convergence of the parameters of a given motion state when the
maximum update is smaller than a threshold
$\boldsymbol{\tau}_{\Delta\boldsymbol{\theta}}=\\{0.05,0.02^{\circ}\mbox{mm}^{-1}\\}\Delta_{\mathbf{y}^{l}}^{\text{max}}$,
with same values used for motion compression. This saves computations by
considering motion updates only on not-converged parameters. However, this
flag is reset to $0$ whenever $i=n(n-1)/2+1$ ($n\in\mathbb{N}_{>0}$) to
account for the impact of the updated reconstructions in the motion parameter
estimates. Joint convergence is achieved when provisional convergence is
achieved for all motion states. Then, the method runs a full CG reconstruction
with the consolidated motion parameters. If regularized outlier rejection
reconstructions are activated, artifacted segments are rejected at levels $l$
such that $\Delta_{\mathbf{y}^{l}}^{\text{max}}\leq 2\,\mbox{mm}$ by densely
sampling within
$[c(1),c(B)]=0.5+(0.35/\max(\Delta_{\mathbf{y}^{l}}^{\text{max}},1))[-1,1]$
and using $\tau_{w}=0.05$. If regularization is applied, shearlets are
designed based on Kutyniok et al.,  (2016) and a final reconstruction is
launched with $3$ CG iterations within $2$ updates of the IRWLS-induced cost
function and $\lambda$. Reconstructions are performed on a $8(16)\times$
Intel(R) Core(TM) i7-5960X $3.00\,\mbox{GHz}$ CPU, $64\,\mbox{GB}$ RAM,
GeForce GTX TITAN X GPU. For further implementation details, readers can refer
to the source code.

##  4 Results

###  4.1 Validation

In Fig. 5 we compare different simulated reconstruction scenarios showing the
losses when iterating the method,
$r_{\hat{\mathbf{x}}^{(i+1)},\hat{\boldsymbol{\theta}}^{(i)}}$ as solid
colored lines with joint iterations represented by markers. Losses in the
convergence plots are normalized to the minimum of the reference levels
$r_{\hat{\mathbf{x}},\boldsymbol{\theta}^{\ast}}$, which are shown as dashed
lines strongly overlapped for the different alternatives. Fig. 6 includes
reconstructions with and without motion correction for different
reconstruction scenarios and provides absolute value error maps and mean SNR
for the compared cases.

\begin{overpic}[width=429.28616pt,draft=false]{Fig05/fig01}\put(0.0,0.0){{a)}}\end{overpic}

\begin{overpic}[width=429.28616pt,draft=false]{Fig05/fig02}\put(0.0,0.0){{b)}}\end{overpic}

\begin{overpic}[width=429.28616pt,draft=false]{Fig05/fig03}\put(0.0,0.0){{c)}}\end{overpic}

Figure 5: Aligned reconstruction convergence against effective iterations $j$
defined as a single application of the encoding or decoding operator for a
single motion state and coil channel at full resolution. a) Different encoding
orders, number of segments $M\in\\{4,64\\}$ (rows), motion levels
$\theta\in\\{2^{\circ},5^{\circ},10^{\circ}\\}$ (columns),
$j_{\text{max}}=20000$. b) Different encoding orders and number of
multiresolution levels $L$, number of segments $M\in\\{16,256\\}$ (rows),
motion levels $\theta\in\\{5^{\circ},10^{\circ},20^{\circ}\\}$ (columns),
$j_{\text{max}}=200000$. c) Different encoding orders, $L=2$, acceleration
factors $R\in\\{1\times 1,2\times 2\\}$ with matched number of segments
$M\in\\{64,16\\}$ (rows), motion levels
$\theta\in\\{5^{\circ},10^{\circ},20^{\circ}\\}$ (columns),
$j_{\text{max}}=200000$.

\begin{overpic}[width=432.31653pt,draft=false]{Fig06/fig01}\end{overpic}

\begin{overpic}[width=433.62pt,draft=false]{Fig06/fig02}\put(0.0,-4.0){{a)}}\footnotesize{
\put(9.2593,-4.0){\makebox[0.0pt]{\Centerstack{No motion\\\
$\mbox{SNR}=29.99\,\mbox{dB}$}}}
\put(27.7778,-4.0){\makebox[0.0pt]{\Centerstack{Uncorrected\\\ {Sequential}\\\
$\mbox{SNR}=12.72\,\mbox{dB}$}}}
\put(46.2963,-4.0){\makebox[0.0pt]{\Centerstack{Corrected\\\ {Sequential}\\\
$\mbox{SNR}=13.37\,\mbox{dB}$}}}
\put(64.8148,-4.0){\makebox[0.0pt]{\Centerstack{Uncorrected\\\ {Random-
checkered}\\\ $\mbox{SNR}=14.51\,\mbox{dB}$}}}
\put(83.3333,-4.0){\makebox[0.0pt]{\Centerstack{Corrected\\\ {Random-
checkered}\\\ $\mbox{SNR}=28.40\,\mbox{dB}$}}} }\end{overpic}

\begin{overpic}[width=432.31653pt,draft=false]{Fig06/fig03}\end{overpic}

\begin{overpic}[width=433.62pt,draft=false]{Fig06/fig04}\put(0.0,-4.0){{b)}}\footnotesize{
\put(7.8182,-4.0){\makebox[0.0pt]{\Centerstack{$R=1\times 1$\\\ No motion\\\
$\mbox{SNR}=29.98\,\mbox{dB}$}}}
\put(23.4545,-4.0){\makebox[0.0pt]{\Centerstack{$R=1\times 1$\\\ Known
motion\\\ $\mbox{SNR}=28.42\,\mbox{dB}$}}}
\put(39.0909,-4.0){\makebox[0.0pt]{\Centerstack{$R=1\times 1$\\\ Estimated
motion\\\ $\mbox{SNR}=28.40\,\mbox{dB}$}}}
\put(54.7273,-4.0){\makebox[0.0pt]{\Centerstack{$R=2\times 2$\\\ No motion\\\
$\mbox{SNR}=20.78\,\mbox{dB}$}}}
\put(70.3636,-4.0){\makebox[0.0pt]{\Centerstack{$R=2\times 2$\\\ Known
motion\\\ $\mbox{SNR}=16.42\,\mbox{dB}$}}}
\put(86.0,-4.0){\makebox[0.0pt]{\Centerstack{$R=2\times 2$\\\ Estimated
motion\\\ $\mbox{SNR}=15.02\,\mbox{dB}$}}} }\end{overpic}

Figure 6: Reconstruction results a) without and with motion correction for
Sequential and Random-checkered traversals compared to the GT (top row) and
corresponding error maps (bottom row) for $M=64$ / $\theta=10^{\circ}$ and b)
for a Random-checkered traversal without motion, and with known and estimated
$\theta=10^{\circ}$ motion for $R=1\times 1$ and $R=2\times 2$ (top row) and
corresponding error maps (bottom row).

####  4.1.1 Encoding orders

Fig. 5a compares the Sequential, Checkered, Random-checkered and Random
traversals. Global convergence is achieved for all considered $M$ and $\theta$
when using any of the Checkered, Random-checkered or Random traversals. In
contrast, when using the Sequential order, the method converges to a local
optimum or fails to converge in the prescribed iterations except for
$\theta\in\\{2^{\circ},5^{\circ}\\}$ / $M=4$. The loss at the first iteration
$r_{\hat{\mathbf{x}}^{(1)},\hat{\boldsymbol{\theta}}^{(0)}}$ is always bigger
when using non-sequential traversals. This increased inconsistency in the
measurement domain relates to the aligned reconstruction sensitivity to motion
degradation.

Fig. 6a shows reconstructions and errors with and without motion correction
for the Sequential and Random-checkered traversals together with GT motion-
free reconstructions. Motion corrected reconstructions using the Random-
checkered data appear similar to the GT despite the strong blurring in
uncorrected reconstructions. This is confirmed by the lack of perceptible
structure in the residuals and a moderate noise amplification. In contrast,
corrections using the Sequential traversal provide only a modest visual
benefit.

####  4.1.2 Multiresolution

Fig. 5b compares the Checkered, Random-checkered and Random traversals when
using a single scale for joint motion estimation and reconstruction ($L=1$)
and when first approximating the motion solution at half the acquired
resolution to initialize the joint problem at full resolution ($L=2$). The
Sequential traversal was excluded because, as discussed when introducing the
multiresolution strategy, it has no opportunity to improve from the poor
relative performance showed in Fig. 5a by exploiting multiresolution. Plots
also include $r_{\hat{\mathbf{x}},\boldsymbol{\theta}^{\ast}}$ at the coarse
scale. Global convergence is achieved for all traversals at all considered
configurations except at $M=256$ / $\theta=20^{\circ}$. However, the
multiresolution strategy ($L=2$) achieves global convergence in less
iterations or provides a solution with lower residuals ($M=256$ /
$\theta=20^{\circ}$). For moderate levels of motion, convergence is generally
quick. For instance, it takes approximately $10$ joint iterations $i$ when
using the Random-checkered traversal in a case where random excursions of up
to $\theta=10^{\circ}$ are imposed in every one of the $M=256$ segments,
probably a more challenging scenario than expected in practice.

####  4.1.3 Acceleration

Fig. 5c tests the ability of the Checkered, Random-checkered and Random
traversals (using $L=2$ scales) to operate in uniformly accelerated regimes as
given by different acceleration factors $R$. We observe convergence to the
global solution in all tested scenarios apart from the Random traversal at
$R=2\times 2$ / $M=16$ / $\theta=20^{\circ}$. Considering all conducted
simulations, the Random-checkered traversal is generally providing the
quickest solutions.

Fig. 6b shows an example of reconstructions and errors in the absence of
motion, with known motion and with estimated motion at $R=1\times 1$ and
$R=2\times 2$. The SNR figures for $R=1\times 1$ and known motion show a
degradation of $1.56\,\text{dB}$ with respect to the reference due to noise
amplification from non-uniform effective k-space sampling after motion. No
further degradation is introduced from motion estimation errors, as
approximately the same SNR figures are obtained for known and estimated
motion. $R=2\times 2$ acceleration in the absence of motion introduces a
degradation of $9.20\,\text{dB}$ with respect to the reference, which stems
from the reduced number of samples and the g-factor (Pruessmann et al., ,
1999). The presence of motion adds further degradations quantified as
$4.36\,\text{dB}$, thus stronger than in the non-accelerated case. Therefore,
without regularization, the limiting reconstruction quality in the presence of
motion decreases with larger distances between neighboring k-space points.
Finally, accelerating the scan has also an impact in the uncertainty of motion
estimates, as we observe a degradation of $1.40\,\text{dB}$ from known to
estimated motion, although the errors show no perceptible structure.

###  4.2 Redundancy for motion tolerance

Fig. 7 compares reconstructions without correction, with inter-shot
corrections and when activating intra-shot corrections in the presence of
extreme motion during in-vivo data acquisition. Intra-shot corrections are
triggered by subsequent temporal binary subdivisions of the sampled
information within each shot until $16$ motion states are estimated per shot.
We show reconstructions without deliberate motion (GT reconstructions), and
reconstructions and absolute differences with respect to the GT using $Q=1$,
$Q=2$ and $Q=3$ repeats under extreme motion (Extreme motion reconstructions /
errors, $Q=\\{1,2,3\\}$). Results for $Q=1$ and $Q=2$ correspond to the first
repeats, with no remarkable differences observed when choosing any other
combination. Reconstructions are provided without regularization or outlier
rejection. Results without deliberate motion show that inter- and intra-shot
corrections do not reduce the reconstruction quality, which demonstrates a
safe application of generalized reconstructions in the absence of motion.
Degradation is noticeable for uncorrected reconstructions in the presence of
motion for all values of $Q$, although with less coherent ghosts as $Q$
increases due to incoherent blurring by Random-checkered motion averaging.
Inter-shot corrections increase the reconstruction quality in all cases, with
more finely resolved cortical structures as $Q$ increases but with noticeably
inferior quality than without deliberate motion. Residual degradation is only
partially accounted when using intra-shot corrections on a single repeat, but
can be more satisfactorily addressed with $Q=2$ and even more with $Q=3$.
Namely, the level of deblurring in the fourth and sixth columns of Fig. 7c
makes corresponding reconstructions visually comparable to those of the first
column despite the extreme and continuous motion (estimated excursions up to
$25^{\circ}$). Thus, we can reason that powerful tolerance is achieved for
$R=2\times 2$ and $Q=2$, so that $Q=1$ with acceleration
$R=\sqrt{2}\times\sqrt{2}$111Both alternatives involve the same scanning time
but the latter generates a lower g-factor. The former was used in this
experiment because it was more convenient to our scanner implementation of the
traversals. may be adequate for motion tolerance in practice, which has been
used to guide the acceleration in the pediatric cohort (see Table 1). However,
in contrast to computation times of $2\,\mbox{min}$ (non-deliberate motion),
and $11\,\mbox{min}$ (extreme motion, $Q=3$) for inter-shot corrections,
corresponding intra-shot corrections required $21\,\mbox{min}$ and
$20\,\mbox{h}\,36\,\mbox{min}$. Thus, despite being technically feasible,
intra-shot corrections may have limited applicability due to computational
costs. Computational cost increase with the complexity of motion is due to the
larger number of iterations for convergence and to the proposed motion
compression strategy, with $13/30$ binned inter-shot motion states without
deliberate motion and $30/30$ with extreme motion ($Q=1$), with proportional
savings in the reconstruction steps in the former.

\begin{overpic}[width=59.83977pt,draft=false]{Fig07/fig01-01}\put(1.0,2.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{a)}}}\end{overpic}

\begin{overpic}[width=59.83977pt,draft=false]{Fig07/fig02-01}\end{overpic}

\begin{overpic}[width=59.83977pt,draft=false]{Fig07/fig03-01}\end{overpic}

\begin{overpic}[width=59.83977pt,draft=false]{Fig07/fig04-01}\end{overpic}

\begin{overpic}[width=59.83977pt,draft=false]{Fig07/fig05-01}\end{overpic}

\begin{overpic}[width=59.83977pt,draft=false]{Fig07/fig06-01}\end{overpic}

\begin{overpic}[width=59.83977pt,draft=false]{Fig07/fig07-01}\end{overpic}

\begin{overpic}[width=59.83977pt,draft=false]{Fig07/fig01-02}\put(1.0,2.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{b)}}}\end{overpic}

\begin{overpic}[width=59.83977pt,draft=false]{Fig07/fig02-02}\end{overpic}

\begin{overpic}[width=59.83977pt,draft=false]{Fig07/fig03-02}\end{overpic}

\begin{overpic}[width=59.83977pt,draft=false]{Fig07/fig04-02}\end{overpic}

\begin{overpic}[width=59.83977pt,draft=false]{Fig07/fig05-02}\end{overpic}

\begin{overpic}[width=59.83977pt,draft=false]{Fig07/fig06-02}\end{overpic}

\begin{overpic}[width=59.83977pt,draft=false]{Fig07/fig07-02}\end{overpic}

\begin{overpic}[width=59.83977pt,draft=false]{Fig07/fig01-03}\put(1.0,2.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{c)}}}\footnotesize{\put(42.85,-26.0){\makebox[0.0pt]{\Centerstack{GT\\\
reconstructions}}}}\end{overpic}

\begin{overpic}[width=59.83977pt,draft=false]{Fig07/fig02-03}\footnotesize{\put(42.85,-26.0){\makebox[0.0pt]{\Centerstack{Extreme
motion\\\ reconstructions\\\ $Q=1$}}}}\end{overpic}

\begin{overpic}[width=59.83977pt,draft=false]{Fig07/fig03-03}\footnotesize{\put(42.85,-26.0){\makebox[0.0pt]{\Centerstack{Extreme
motion\\\ errors\\\ $Q=1$}}}}\end{overpic}

\begin{overpic}[width=59.83977pt,draft=false]{Fig07/fig04-03}\footnotesize{\put(42.85,-26.0){\makebox[0.0pt]{\Centerstack{Extreme
motion\\\ reconstructions\\\ $Q=2$}}}}\end{overpic}

\begin{overpic}[width=59.83977pt,draft=false]{Fig07/fig05-03}\footnotesize{\put(42.85,-26.0){\makebox[0.0pt]{\Centerstack{Extreme
motion\\\ errors\\\ $Q=2$}}}}\end{overpic}

\begin{overpic}[width=59.83977pt,draft=false]{Fig07/fig06-03}\footnotesize{\put(42.85,-26.0){\makebox[0.0pt]{\Centerstack{Extreme
motion\\\ reconstructions\\\ $Q=3$}}}}\end{overpic}

\begin{overpic}[width=59.83977pt,draft=false]{Fig07/fig07-03}\footnotesize{\put(42.85,-26.0){\makebox[0.0pt]{\Centerstack{Extreme
motion\\\ errors\\\ $Q=3$}}}}\end{overpic}

Figure 7: Reconstruction results for extreme motion in vivo. a) Uncorrected;
b) inter-shot corrections; and c) intra-shot corrections. From left to right,
results without deliberate motion and reconstructions and errors in the
presence of extreme motion for $Q=\\{1,2,3\\}$ repeats of a $R=2\times 2$
accelerated baseline scan.

###  4.3 Non-compliant subjects

Fig. 8 shows worst-case reconstructions without motion correction, with motion
correction alone and with motion correction and the regularized outlier
segment rejection. Results are shown for main structural brain MRI modalities,
MP-RAGE, TSE, FLAIR, SPGR and bSSFP. In all sequences we observe a substantial
improvement when activating motion-corrected reconstructions alone, with
better delineated cortical structures. However, subtle artifacts are still
present, either in the form of ghosts or of coloured noise. Fig. 8c shows that
quality can be further improved by rejecting the less consistent segments and
performing a regularized reconstruction. In some sequences discarding the
artifacted segments seems to reduce residual artifacts from uncorrected fast
motion (see for instance fine details in SPGR) while in others it seems to
mainly improve the magnetization consistency (see TSE contrast). Across the
cohort we have observed that motion artifact levels always decrease when
compensating for motion, with no remarkable differences when activating the
corrections in the absence of artifacts. This is along the lines of the
quantitative population metrics obtained for the less favourable sequential
sampling (Cordero-Grande et al., , 2016) or for multi-slice scans (Cordero-
Grande et al., , 2018). Worst-case results of Fig. 8 have been judged
satisfactory by the practitioners and researchers involved in the project.
Therefore, the proposed methodology is delivering reliable examinations for
unsedated pediatric subjects challenging to comply to the MRI motion
requirements. In this experiment motion estimates were performed at half the
acquired resolution with joint motion estimation and reconstruction always
taking less time than final reconstructions at full resolution. Computation
times range between $5\,\mbox{min}$ in least artifacted and $40\,\mbox{min}$
in most artifacted volumes in our cohort. Estimated motion traces and outlier
segments are reported in Supporting Information Fig. S1.

\begin{overpic}[width=84.55907pt,draft=false]{Fig08/fig01-01}\put(1.0,2.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{a)}}}\end{overpic}

\begin{overpic}[width=84.55907pt,draft=false]{Fig08/fig02-01}\end{overpic}

\begin{overpic}[width=84.55907pt,draft=false]{Fig08/fig03-01}\end{overpic}

\begin{overpic}[width=84.55907pt,draft=false]{Fig08/fig04-01}\end{overpic}

\begin{overpic}[width=84.55907pt,draft=false]{Fig08/fig05-01}\end{overpic}

\begin{overpic}[width=84.55907pt,draft=false]{Fig08/fig01-02}\put(1.0,2.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{b)}}}\end{overpic}

\begin{overpic}[width=84.55907pt,draft=false]{Fig08/fig02-02}\end{overpic}

\begin{overpic}[width=84.55907pt,draft=false]{Fig08/fig03-02}\end{overpic}

\begin{overpic}[width=84.55907pt,draft=false]{Fig08/fig04-02}\end{overpic}

\begin{overpic}[width=84.55907pt,draft=false]{Fig08/fig05-02}\end{overpic}

\begin{overpic}[width=84.55907pt,draft=false]{Fig08/fig01-03}\put(1.0,2.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{c)}}}\footnotesize{\put(37.5,-9.0){\makebox[0.0pt]{\Centerstack{MP-
RAGE}}}}\end{overpic}

\begin{overpic}[width=84.55907pt,draft=false]{Fig08/fig02-03}\footnotesize{\put(37.5,-9.0){\makebox[0.0pt]{\Centerstack{TSE}}}}\end{overpic}

\begin{overpic}[width=84.55907pt,draft=false]{Fig08/fig03-03}\footnotesize{\put(37.5,-9.0){\makebox[0.0pt]{\Centerstack{FLAIR}}}}\end{overpic}

\begin{overpic}[width=84.55907pt,draft=false]{Fig08/fig04-03}\footnotesize{\put(37.5,-9.0){\makebox[0.0pt]{\Centerstack{SPGR}}}}\end{overpic}

\begin{overpic}[width=84.55907pt,draft=false]{Fig08/fig05-03}\footnotesize{\put(37.5,-9.0){\makebox[0.0pt]{\Centerstack{bSSFP}}}}\end{overpic}

Figure 8: Reconstruction results for pediatric cases with largest intra-scan
degradations. a) Uncorrected; b) motion-corrected; and c) motion-corrected and
regularized outlier segment rejection reconstructions. From left to right,
results for the main families of sequences for structural brain imaging.

##  5 Discussion

We have presented DISORDER, a retrospective framework for motion tolerant
structural 3D k-space encoded brain imaging that combines optimized view
orders with an improved aligned reconstruction. The proposed distributed and
incoherent orders increase the motion sensitivity of the information sampled
within a given time window which, provided a certain degree of redundancy,
enables the resolvability of motion in the reconstruction. Conducted
simulations have shown that reordering the k-space traversals introduces a
significant boost in the ability to estimate the head pose and suppress motion
artifacts. Tolerance to motion has been demonstrated in-vivo on a controlled
experiment involving extreme and continuous motion throughout the examination
as well as for the main families of sequences used for structural brain
imaging by presenting the reconstruction results on the most challenging
datasets from a pediatric cohort of $26$ subjects.

Although DISORDER is robust enough in its current form so as to be of
practical interest for reliable structural brain MR examinations in non-
compliant cohorts, with plans in our center to use it to progressively replace
unnecessary sedation or anesthesia in pediatric and neonatal populations
(Barton et al., , 2011), it is obviously not free from limitations. First,
data consistency may be affected by additional degrading factors. These
include inaccuracies in sensitivities but also water-fat shifts, eddy currents
or flow artifacts. In practice, applying fat suppression when possible,
designing the tiles for adequate trade-offs between eddy currents and motion
resolvability in bSSFP sequences, and adequate planning and scanning
procedures are usually sufficient to address these issues. Differently,
correction of non-rigid motion components would require an extension of the
formulation. Although analogous methodologies (Loktyushin et al., , 2015) have
shown potential in this context, a robust and efficient extension to non-rigid
motion models will probably require a careful computational design. This may
be particularly the case for high resolution applications, where both rigid
and non-rigid motion become more important and challenging to correct (Budde
et al., , 2014). Moreover, coarse scale motion at ultra high field may require
additional corrections of high order effects. Finally, in this manuscript we
have restricted ourselves to uniform sampling, with further work required to
generalize the incoherent and distributed orders and characterize motion
correction and resolution retrieval when using variable densities.

In the in-vivo experiments of Fig. 8 we have shown that inter-shot corrections
can be sufficient in practical brain imaging scenarios requiring motion
tolerance. Our underlying assumption is that the subject remains approximately
still for a significant portion of the acquisition. In this case, inter-shot
corrections are enough to reconcile the brain pose amongst the stable periods
and data rejection can be applied to the transitions, again, provided that
sampling is redundant enough. However, intra-shot corrections may become more
important in challenging situations, as illustrated in Fig. 7. Despite its
computational limitations, our method is able to provide stable intra-shot
estimates in the absence of motion while offering certain motion correction
potential. Although a prior model for the temporal evolution of motion may aid
in certain applications, in general, limitations arising from available
computational resources and SNR per motion state are likely to complicate
intra-shot tractability.

The situation may perhaps be different if using supervised learning strategies
to inform the exploration of the motion parameter space. These may help to
improve the spatio-temporal resolvability of motion by aiding the intra-shot
corrections to find better motion solutions. Training may also help to enlarge
the motion capture range at a given level of redundancy or decrease the
required level of redundancy for a given motion capture range. Although direct
learning of motion-corrected reconstructions could also be attempted, it is
likely that, in many circumstances, better results will be obtained when
concatenating learned reconstructions with model-based strategies, as recently
suggested in Haskell et al.,  (2019). Further integration of both approaches
could be tackled, for instance, by incorporating the motion operator into the
model-based learning framework in Aggarwal et al.,  (2019), which may be
effective in dealing with the residual penalties from g-factor amplification
due to motion (see Fig. 6b). Thereby, future work will explore the
opportunities for extending the ranges of motion resilience by supervised
learning.

##  6 Conclusion

We have proposed a simple modification of standard 3D Cartesian sequences for
structural brain imaging, involving only a distributed and incoherent
reordering of the sampled profiles, for high quality imaging in the presence
of motion. Improved convergence has been demonstrated when using a separable
nonlinear least squares formulation for joint motion estimation and
reconstruction. Feasibility and conditions for inter- and intra-shot
corrections have been characterized by simulations and in-vivo reconstructions
under extreme motion. The DISORDER method has been successfully applied to
replace sedation in a pediatric population scanned using common clinical
examination protocols by combining inter-shot corrections with regularized
outlier segment rejection reconstructions. Future work will focus on applying
DISORDER to other cohorts and on strengthening its performance by integrating
motion learning strategies.

## Acknowledgments

This work received funding from the European Research Council under the
European Union’s Seventh Framework Programme (FP7/20072013/ERC, grant
agreement no. [319456], dHCP project). The research was supported by the
Wellcome/EPSRC Centre for Medical Engineering at King’s College London [WT
203148/Z/16/Z]; the Medical Research Council [MR/K006355/1]; and the National
Institute for Health Research (NIHR) Biomedical Research Centre based at Guy’s
and St Thomas’ NHS Foundation Trust and King’s College London. Jonathan
O’Muircheartaigh is supported by a Sir Henry Dale Fellowship jointly funded by
the Wellcome Trust and the Royal Society [206675/Z/17/Z]. The views expressed
are those of the authors and not necessarily those of the NHS, the NIHR or the
Department of Health. The authors acknowledge the Department of Perinatal
Imaging & Health at King’s College London.

## References

  * Adcock et al., (2014) Adcock, B., Hansen, A., Roman, B., and Teschke, G. (2014).  Generalized sampling: Stable reconstructions, inverse problems and compressed sensing over the continuum.  Adv. Imag. Electron Phys., 182(1):187–279.
  * Aggarwal et al., (2019) Aggarwal, H. K., Mani, M. P., and Jacob, M. (2019).  MoDL: Model based deep learning architecture for inverse problems.  IEEE Trans. Med. Imaging, 38(2):394–405.
  * Aksoy et al., (2012) Aksoy, M., Forman, C., Straka, M., Çukur, T., Hornegger, J., and Bammer, R. (2012).  Hybrid prospective and retrospective head motion correction to mitigate cross-calibration errors.  Magn. Reson. Med., 67(5):1237–1251.
  * Allison et al., (2013) Allison, M. J., Ramani, S., and Fessler, J. A. (2013).  Accelerated regularized estimation of MR coil sensitivities using augmented Lagrangian methods.  IEEE Trans. Med. Imaging, 32(3):556–564.
  * Anderson III et al., (2013) Anderson III, A. G., Velikina, J., Block, W., Wieben, O., and Samsonov, A. (2013).  Adaptive retrospective correction of motion artifacts in cranial MRI with multicoil three-dimensional radial acquisitions.  Magn. Reson. Med., 69(4):1094–1103.
  * Andre et al., (2015) Andre, J. B., Bresnahan, B. W., Mossa-Basha, M., Hoff, M. N., Smith, C. P., Anzai, Y., and Cohen, W. A. (2015).  Towards quantifying the prevalence, severity, and cost associated with patient motion during clinical MR examinations.  J. Am. Coll. Radiol., 12(7):689–695.
  * Avron and Toledo, (2011) Avron, H. and Toledo, S. (2011).  Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix.  J. Assoc. Comput. Mach., 58(2):8:1–8:17.
  * Bammer et al., (2007) Bammer, R., Aksoy, M., and Liu, C. (2007).  Augmented generalized SENSE reconstruction to correct for rigid body motion.  Magn. Reson. Med., 57(1):90–102.
  * Barton et al., (2011) Barton, K., Nickerson, J. P., Higgins, T., and Williams, R. K. (2011).  Pediatric anesthesia and neurotoxicity: what the radiologist needs to know.  Pediatr. Radiol., 48(1):31–36.
  * Batchelor et al., (2005) Batchelor, P. G., Atkinson, D., Irarrazaval, P., Hill, D. L. G., Hajnal, J., and Larkman, D. (2005).  Matrix description of general motion correction applied to multishot images.  Magn. Reson. Med., 54(5):1273–1280.
  * Budde et al., (2014) Budde, J., Shajan, G., Scheffler, K., and Pohmann, R. (2014).  Ultra-high resolution imaging of the human brain using acquisition-weighted imaging at 9.4T.  NeuroImage, 86(1):592–598.
  * Buehrer et al., (2007) Buehrer, M., Pruessmann, K. P., Boesiger, P., and Kozerke, S. (2007).  Array compression for MRI with large coil arrays.  Magn. Reson. Med., 57(6):1131–1139.
  * Busse et al., (2008) Busse, R. F., Brau, A. C., Vu, A., Michelich, C. R., Bayram, E., Kijowski, R., Reeder, S. B., and Rowley, H. A. (2008).  Effects of refocusing flip angle modulation and view ordering in 3D fast spin echo.  Magn. Reson. Med., 60(3):640–649.
  * Bydder et al., (2002) Bydder, M., Larkman, D. J., and Hajnal, J. V. (2002).  Detection and elimination of motion artifacts by regeneration of $k$-space.  Magn. Reson. Med., 47(4):677–686.
  * Cordero-Grande et al., (2018) Cordero-Grande, L., Hughes, E. J., Hutter, J., Price, A. N., and Hajnal, J. V. (2018).  Three-dimensional motion corrected sensitivity encoding reconstruction for multi-shot multi-slice MRI: Application to neonatal brain imaging.  Magn. Reson. Med., 79(3):1365–1376.
  * Cordero-Grande et al., (2016) Cordero-Grande, L., Teixeira, R. P. A. G., Hughes, E. J., Hutter, J., Price, A. N., and Hajnal, J. V. (2016).  Sensitivity encoding for aligned multishot magnetic resonance reconstruction.  IEEE Trans. Comput. Imaging, 2(3):266–280.
  * Fuderer et al., (2004) Fuderer, M., van den Brink, J., and Jurrissen, M. (2004).  SENSE reconstruction using feed forward regularization.  In 12th Proc. Intl. Soc. Mag. Reson. Med., page 2130, Kyoto, Japan.
  * Gallichan and Marques, (2017) Gallichan, D. and Marques, J. P. (2017).  Optimizing the acceleration and resolution of three-dimensional fat image navigators for high-resolution motion correction at 7T.  Magn. Reson. Imaging, 77(2):547–558.
  * Gan et al., (2018) Gan, M., Chen, C. L. P., Chen, G.-Y., and Chen, L. (2018).  On some separated algorithms for separable nonlinear least squares problems.  IEEE Trans. Cybern., 48(10):2866–2874.
  * Godenschweger et al., (2016) Godenschweger, F., Kägebein, U., Stucht, D., Yarach, U., Sciarra, A., Yakupov, R., Lüsebrink, F., Schulze, P., and Speck, O. (2016).  Motion correction in MRI of the brain.  Phys. Med. Biol., 61(5):R32–R56.
  * Haskell et al., (2019) Haskell, M. W., Cauley, S. F., Bilgic, B., Hossbach, J., Splitthoff, D. N., Pfeuffer, J., Setsompop, K., and Wald, L. L. (2019).  Network accelerated motion estimation and reduction (NAMER): Convolutional neural network guided retrospective motion correction using a separable motion model.  Magn. Reson. Med., 82(4):1452–1461.
  * Herring et al., (2018) Herring, J. L., Nagy, J. G., and Ruthotto, L. (2018).  LAP: A linearize and project method for solving inverse problems with coupled variables.  Sampling Theory Signal Image Process., 17(2):127–151.
  * Johnson et al., (2016) Johnson, P. M., Liu, J., Wade, T., Tavallaei, M. A., and Drangova, M. (2016).  Retrospective 3D motion correction using spherical navigator echoes.  Magn. Reson. Imaging, 34(9):1274–1282.
  * Kutyniok et al., (2012) Kutyniok, G., Lemvig, J., and Lim, W.-Q. (2012).  Optimally sparse approximations of 3D functions by compactly supported shearlet frames.  SIAM J. Math. Anal., 44(4):2962–3017.
  * Kutyniok et al., (2016) Kutyniok, G., Lim, W.-Q., and Reisenhofer, R. (2016).  ShearLab 3D: Faithful digital shearlet transforms based on compactly supported shearlets.  ACM Trans. Math. Softw., 42(1):5:1–5:42.
  * Li and Zhong, (2018) Li, C.-J. and Zhong, Y.-J. (2018).  A pseudo-heuristic parameter selection rule for $l^{1}$-regularized minimization problems.  J. Comput. Appl. Math., 333:1–19.
  * Loktyushin et al., (2013) Loktyushin, A., Nickisch, H., Pohmann, R., and Schölkopf, B. (2013).  Blind retrospective motion correction of MR images.  Magn. Reson. Med., 70(6):1608–1618.
  * Loktyushin et al., (2015) Loktyushin, A., Nickisch, H., Pohmann, R., and Schölkopf, B. (2015).  Blind multirigid retrospective motion correction of MR images.  Magn. Reson. Med., 73(4):1457–1468.
  * Maclaren et al., (2013) Maclaren, J., Herbst, M., Speck, O., and Zaitsev, M. (2013).  Prospective motion correction in brain imaging: A review.  Magn. Reson. Med., 69(3):621–636.
  * Odille et al., (2008) Odille, F., Vuissoz, P.-A., Marie, P.-Y., and Felblinger, J. (2008).  Generalized reconstruction by inversion of coupled systems (GRICS) applied to free-breathing MRI.  Magn. Reson. Med., 60(1):146–157.
  * Pipe et al., (2014) Pipe, J. G., Gibbs, W. N., Li, Z., Karis, J. P., Schar, M., and Zwart, N. R. (2014).  Revised motion estimation algorithm for PROPELLER MRI.  Magn. Reson. Med., 72(2):430–437.
  * Pruessmann et al., (2001) Pruessmann, K. P., Weiger, M., Börnert, P., and Boesiger, P. (2001).  Advances in sensitivity encoding with arbitrary k-space trajectories.  Magn. Reson. Med., 46(4):638–651.
  * Pruessmann et al., (1999) Pruessmann, K. P., Weiger, M., Scheidegger, M. B., and Boesiger, P. (1999).  SENSE: sensitivity encoding for fast MRI.  Magn. Reson. Med., 42(5):952–962.
  * Ruano et al., (1991) Ruano, A. E. B., Jones, D. I., and Fleming, P. J. (1991).  A new formulation of the learning problem of a neural network controller.  In Proc. 30th IEEE Conf. Decision Control, pages 865–866, Brighton, England.
  * Samsonov et al., (2010) Samsonov, A. A., Velikina, J., Jung, Y., Kholmovski, E. G., Johnson, C. R., and Block, W. F. (2010).  POCS-enhanced correction of motion artifacts in parallel MRI.  Magn. Reson. Med., 63(4):1104–1110.
  * Schulz et al., (2012) Schulz, J., Siegert, T., Reimer, E., Labadie, C., Maclaren, J., Herbst, M., Zaitsev, M., and Turner, R. (2012).  An embedded optical tracking system for motion-corrected magnetic resonance imaging at 7T.  Magn. Reson. Mater. Phy., 25(6):443–453.
  * Tisdall et al., (2016) Tisdall, M. D., Reuter, M., Qureshi, A., Buckner, R. L., Fischl, B., and van der Kouwe, A. J. (2016).  Prospective motion correction with volumetric navigators (vNavs) reduces the bias and variance in brain morphometry induced by subject motion.  NeuroImage, 127:11–22.
  * Tsao et al., (2005) Tsao, J., Kozerke, S., Boesiger, P., and Pruessmann, K. P. (2005).  Optimizing spatiotemporal sampling for k-t BLAST and k-t SENSE: Application to high-resolution real-time cardiac steady-state free precession.  Magn. Reson. Med., 53(6):1372–1382.
  * Unser et al., (1993) Unser, M., Aldroubi, A., and Gerfen, C. R. (1993).  A multiresolution image registration procedure using spline pyramids.  In Proc. SPIE, Wavelet Applications in Signal and Image Processing, volume 2034, pages 160–170.
  * Unser et al., (1995) Unser, M., Thévenaz, P., and Yaroslavsky, L. (1995).  Convolution-based interpolation for fast, high-quality rotation of images.  IEEE Trans. Image Process., 4(10):1375–1381.
  * Voronin and Daubechies, (2017) Voronin, S. and Daubechies, I. (2017).  Functional Analysis, Harmonic Analysis, and Image Processing, volume 693 of Contemp. Math., chapter An iteratively reweighted least squares algorithm for sparse regularization, pages 391–411.  Am. Math. Soc.
  * Zaitsev et al., (2015) Zaitsev, M., Maclaren, J., and Herbst, M. (2015).  Motion artifacts in MRI: A complex problem with many partial solutions.  J. Magn. Reson. Imaging, 42(4):887–901.

## Supporting information

Fig. S1a collects the estimated motion traces for the cases in Fig. 8.
Although no temporal regularization is used, all traces show periods of
stability, which suggests accurate estimates at least in these periods. In
Fig. S1b the opacity of the traces is driven by the segment weights from the
proposed outlier detection method. We observe that outliers generally
correspond to main motion transients, in agreement to higher chances for
intra-sweep degradation in these periods.

\begin{overpic}[width=177.78578pt,draft=false]{Fig09/fig01-01}\put(-16.0,28.0){\makebox[0.0pt]{\Centerstack{MP-
RAGE}}}\end{overpic}

\begin{overpic}[width=177.78578pt,draft=false]{Fig09/fig01-02}\end{overpic}

\begin{overpic}[width=177.78578pt,draft=false]{Fig09/fig02-01}\put(-16.0,28.0){\makebox[0.0pt]{\Centerstack{TSE}}}\end{overpic}

\begin{overpic}[width=177.78578pt,draft=false]{Fig09/fig02-02}\end{overpic}

\begin{overpic}[width=177.78578pt,draft=false]{Fig09/fig03-01}\put(-16.0,28.0){\makebox[0.0pt]{\Centerstack{FLAIR}}}\end{overpic}

\begin{overpic}[width=177.78578pt,draft=false]{Fig09/fig03-02}\end{overpic}

\begin{overpic}[width=177.78578pt,draft=false]{Fig09/fig04-01}\put(-16.0,28.0){\makebox[0.0pt]{\Centerstack{SPGR}}}\end{overpic}

\begin{overpic}[width=177.78578pt,draft=false]{Fig09/fig04-02}\end{overpic}

\begin{overpic}[width=177.78578pt,draft=false]{Fig09/fig05-01}\put(-16.0,28.0){\makebox[0.0pt]{\Centerstack{BSSFP}}}\put(1.0,-6.0){{a)}}\end{overpic}

\begin{overpic}[width=177.78578pt,draft=false]{Fig09/fig05-02}\put(1.0,-6.0){{b)}}\end{overpic}

Figure S1. First $120\,\mbox{s}$ of estimated motion traces for pediatric
cases with largest intra-scan degradations for each sequence. Left: original
motion traces. Right: motion traces with segment opacity given by
corresponding reliability $w$. Solid lines indicate data collection periods
for each segment with dotted lines used to connect these. Tra and Rot in the
legends refer respectively to Translation and Rotation.