Three‐dimensional motion corrected sensitivity encoding reconstruction for multi‐shot multi‐slice MRI: Application to neonatal brain imaging

Purpose To introduce a methodology for the reconstruction of multi‐shot, multi‐slice magnetic resonance imaging able to cope with both within‐plane and through‐plane rigid motion and to describe its application in structural brain imaging. Theory and Methods The method alternates between motion estimation and reconstruction using a common objective function for both. Estimates of three‐dimensional motion states for each shot and slice are gradually refined by improving on the fit of current reconstructions to the partial k‐space information from multiple coils. Overlapped slices and super‐resolution allow recovery of through‐plane motion and outlier rejection discards artifacted shots. The method is applied to T 2 and T 1 brain scans acquired in different views. Results The procedure has greatly diminished artifacts in a database of 1883 neonatal image volumes, as assessed by image quality metrics and visual inspection. Examples showing the ability to correct for motion and robustness against damaged shots are provided. Combination of motion corrected reconstructions for different views has shown further artifact suppression and resolution recovery. Conclusion The proposed method addresses the problem of rigid motion in multi‐shot multi‐slice anatomical brain scans. Tests on a large collection of potentially corrupted datasets have shown a remarkable image quality improvement. Magn Reson Med 79:1365–1376, 2018. © 2017 The Authors Magnetic Resonance in Medicine published by Wiley Periodicals, Inc. on behalf of International Society for Magnetic Resonance in Medicine. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.


INTRODUCTION
Common magnetic resonance imaging (MRI) acquisition protocols take seconds to minutes to complete so that the image quality is vulnerable to subject motion during the scan. Strategies for robustness against motion can be adopted at different levels. First, when possible, motion occurrence can be prevented or limited at the time of scanning. Second, sequence design can be made as tolerant as possible against motion disturbance, which should be made without compromising scan efficiency. Third, remaining artifacts can be reduced and sampled information can be optimally integrated when reconstructing the data.
As for brain imaging, although non-rigid motion (1) can occur due to pulsatile motion (2) and localized movements can arise such as from eyeball motion (3), the current clinical paradigm relies on the subjects holding their head still enough for the duration of the acquisition to avoid gross motion artefacts. However, in many brain studies, such as for those subjects who have difficulty remaining still enough (4) and in ultra-high resolution imaging applications (5), large or small motion inconsistencies become a limiting factor, so that appropriate acquisition or reconstruction strategies have to be put in place.
Usual sequences in structural brain studies are based on multi-shot methods, where a fraction of the k-space is acquired after a single radio frequency excitation or shot. This type of sampling, particularly when combined with slice selective excitation, provides flexibility to achieve the desired contrast while simultaneously balancing the signal to noise ratio, image resolution, and scanning time requirements. A common operating regime is to use an interleaved scheme in which several slices are sequentially excited within a given repeat time (T R ), that way permitting efficient sampling for large T R 's (6). In this setting, changes in the head position among different shots and slices, which may be acquired very distant in time, will provoke degradation of individual slices and inconsistencies in volumetric information.
Head motion estimation and/or retrospective motioncompensated reconstruction in multi-shot sequences has been extensively studied in the past (7)(8)(9)(10)(11)(12)(13)(14)(15). However, to the best of our knowledge, none of the proposed reconstruction methods can be used to correct for throughplane motion in multi-slice acquisitions (1,3). On the other hand, image-based alignment methods have been introduced to assemble the information coming from snapshot acquisitions, where individual slices are acquired fast enough to approximately freeze movement, into self-consistent three-dimensional (3D) representations of the imaged structures. This has been performed either by matching the structures along slice intersections (16,17) or by slice to volume registration (18)(19)(20).
In Ref. 15, we proposed a framework for rigid body motion estimation and motion-compensated reconstruction in volumetric multi-shot sequences using parallel imaging that is grounded on the sensitivity encoding (SENSE) reconstruction paradigm (21). The method is fully data driven, using a common functional to estimate the motion and the image. Moreover, it is suitable for a wide range of settings as it does not require external sensors, navigators, or particular samplings. In that work, we provided an empirical characterization of the conditions for which, in the absence of noise, fully rigid motion corrected reconstructions are still possible in terms of the amount of motion, number of shots, encoding trajectories, and availability of prior information. Here, the framework is extended to multi-slice sequences, thus fusing the multishot and aligned snapshot families of motion compensation methods. Moreover, differently from most imagebased alignment methods, our procedure aims to correct for within-plane and moderate through-plane motion requiring only a single slice orientation. To support this aim, on the acquisition side, slices are sampled in an overlapped manner (where the slice separation is maintained while the slice thickness is increased for same total acquisition time). On the reconstruction side, first, a motion corrected super-resolution technique is used to recover from through-plane motion while considering slice blurring. Second, an outlier rejection strategy is introduced to discard those shots for which the basic assumptions of our framework are not satisfied. Finally, the information from different receiver coils in a head array is used to detect magnetic inconsistencies, to estimate rigid motion, and to interpolate the corrupted spectral information of outlier shots. Once slice and volumetric consistency is improved this way, if two or more orthogonal stacks are available, the result could be connected with a slice to volume reconstruction procedure to obtain 3D nearly isotropic representations of multi-shot multi-slice multi-view datasets. The proposed methodology is applied to the problem of motion correction in T 2 -and inversion recovery T 1 -weighted fast spin echo scans. The approach is tested in a collection of 517 examinations of neonates studied in natural sleep. Motion correction performance is quantified using two complementary image quality metrics and the ability to correct for inter-shot and inter-slice motion, both in-plane and through-plane, as well as for artifacted shots is visually illustrated. The source code of a MATLAB implementation of the reconstruction method and the required data to reproduce the results in Figures 3 and 4 of "Visual Validation" section is made available at https://github.com/ mriphysics/multiSliceAlignedSENSE/releases/tag/1.0.1.

Rigid Motion Corrected Multi-Shot Reconstruction
As introduced in Ref. 15, assuming noise pre-whitening, the reconstruction with rigid motion correction for parallel multi-shot volumetric imaging can be formulated in matrix form as: ðx;TÞ ¼ argmin where y denotes the measured data, x the image to be reconstructed, T a rigid motion transformation, S a coil sensitivity operator, F a discrete Fourier transform (DFT), and A a sampling mask. This formulation is used to reconstruct a 3D image of size N ¼ N 1 N 2 N 3 with N i the number of voxels along dimension i using a coil array of C elements from M ¼ ESC samples of a discretized k-space grid of size y is a vector of size M Â 1.
A is comprised of SC Â SC submatrices of size E Â K that take the value 1 if the sample e of the shot s corresponds to the sampled k-space location indexed by k and 0 otherwise. F is comprised of SC Â SC submatrices of size K Â N corresponding to the 3D DFT with applied k-space oversampling or downsampling. S is comprised of SC Â S diagonal submatrices of size N Â N whose diagonal elements correspond to the spatial profile of a given coil c.
T is comprised of S Â 1 unitary submatrices of size N Â N corresponding to the 3D rigid transformation the underlying structure has been subject to when acquiring the shot s. A convolution-based interpolation technique (22) is used to build this matrices, so that the acquired resolution is fully preserved.
x is a vector of size N Â 1.

Rigid Motion Corrected Multi-Shot Multi-Slice Reconstruction
The extension of the model in Equation [1] to multi-shot multi-slice motion correction in the context of both conventional slice selective excitations as well as simultaneous multi-slice or multiband MRI (23) can be generically formulated as: ðx;TÞ ¼ argmin where some new terms have been introduced: H, a slice profile filter, W, a stabilizer of the reconstruction, and k, a parameter that weights the degree of stabilization required. The sampling matrix A has been split into a spectral part, A ðsÞ , and a spatial one, A ðrÞ . In this case, assuming that the slices are arranged along the third dimension, i ¼ 3, we want to reconstruct a 3D image from M ¼ ESRC samples of a discretized hybrid k-space where now S is the number of shots per simultaneously excited slices, R ¼ P=Q is the number of excitations to cover the field of view (FOV) in the slice direction, with P the total number of excited slices and Q the number of simultaneously excited slices or multiband factor, and K 3 corresponds to the blipping pattern cycling period (23). The structure of those terms involved in Equation [2] that differ from their counterparts in Equation [1] is as follows: A ðsÞ is comprised of SRC Â SRC submatrices of size E Â K that take the value 1 if the sample e of the shot s corresponds to the sampled k-space location indexed by k and 0 otherwise. F is comprised of SRC Â SRC submatrices of size K ÂN 1 N 2 Q corresponding to the 3D DFT with applied k-space resampling and multiband blipping pattern.
A ðrÞ is comprised of N 1 N 2 SRC Â N 1 N 2 SRC submatrices of size Q Â P that take the value 1 if the slice indexed by q in the excitation r corresponds to the slice p and 0 otherwise. H is comprised of N 1 N 2 SRC Â N 1 N 2 SRC submatrices H of size P Â N 3 that account for the slice profile filter and are further discussed below. W is comprised of N 1 N 2 Â N 1 N 2 submatricesW of size N 3 Â N 3 that account for the reconstruction stabilization to be applied in the slice direction and are further discussed below.

Slice Profile Filter
A signal processing model for multi-slice acquisitions was proposed in Ref. 24. It assumes that the underlying continuous object has first been convolved by a slice profile and then it has been sampled at some slice locations. The effect of the slice profile is usually that of a low pass filter on the slice direction whereas the sampling period (i.e., slice distance) is inversely related to the spectral overlap. In our reconstruction technique, we will assume that the reconstructed resolution is dense enough so as to neglect the spectral overlapping of those harmonics above the Nyquist limit at that resolution. Then, the submatricesH represent a convolution filter and a spectral overlapping from the reconstructed to the acquired grid, so they can be diagonal- , withH D a diagonal matrix of size N 3 Â N 3 whose diagonal elements represent the slice profile filter, F N3!N3 the 1D DFT at the reconstructed resolution and F H N3!P the inverse 1D DFT with applied resampling from reconstructed to acquired slice resolution.*

Reconstruction Stabilizer
In Equation [2], we have chosen a stabilizer based on the ' 2norm. This is justified by the large computational demands involved in our procedure, so that a closed-form solution for x becomes highly desirable (see "Problem Solving" section below). The role of the stabilizer is twofold. On the one hand, it aids the treatment of the ill-conditioned deconvolution of the high spatial frequencies that were attenuated by the slice profile filter, that is, it underpins a certain degree of super-resolution. On the other hand, in the presence of through-plane motion, the sampling density in the slice direction can no longer be assumed homogeneous so that supporting a certain degree of slice interpolation appears necessary. This motivated the selection of a second-order finite difference regularizer for this term, as its larger support as compared to the first-order counterpart was observed to help in slice interpolation.
We should remark that the problem formulation in Equation [2] remains intractable due to the huge size of the matrices involved. Thus, certain term rearrangements and approximations are required to reduce the computational burden. These will be described in the "Algorithmic Implementation" section.

Problem Solving
The sampling structure of our problem can be described with a compound matrix E: The joint problem in Equation [2] may be addressed in an alternating fashion by iteratively solving the following sub-problems: The first sub-problem admits a closed form solution, which can be computed using the conjugate gradient algorithm.
The solution of the second sub-problem must satisfy: r ðzÞ jjETðzÞx À yjj 2 2 ¼ 0; [6] with z representing the set of translation and rotation parameters of the transformations involved. The solution can be obtained separately for each excitation and shot, r ðzrsÞ jjE rsTrs ðz rs Þx À y rs jj 2 2 ¼ 0; [7] where z rs is a 6-component vector describing the rigid transformation for excitation r and shot s. We have resorted to the Newton's method to solve this system of equations. The expressions for the gradient and Hessian can be consulted in Ref. 15. An outlier rejection criterion is introduced to improve the robustness of the reconstruction against the violation of the model assumptions (either provoked by within-shot motion, spin history, physiological motion, or other effects). Thus, a feature for outlier rejection is built using where PðrÞ denotes the set of slices P corresponding to excitation r and P fÀ1;þ1g ðrÞ ¼ P À1 ðrÞ [ P þ1 ðrÞ, with P j ðrÞ *We drop bold notation from these DFT's to avoid confusion with the generic DFT matrices in Equations [1] and [2].
the set of slices corresponding to excitation r with their indexes displaced by j. Therefore, this outlier rejection feature is based on the similar nature of the information contained in a given shot for adjacent (j ¼ fÀ1; þ1g) slices. On this basis, a given shot is rejected if t rs > t. In our implementation, outlier detection is performed at each iteration just after the motion estimation step. Then, shots classified as outliers are removed in the reconstruction step at next iteration (by making A ðsÞ rs ¼ 0). The information from these artifacted shots is retrieved by interpolation of neighboring spectral information using the spatial encoding of the coil array and the spatial information from adjacent slices. However, motion estimation is still performed for outliers in subsequent iterations and if it turns that the new motion estimation is able to explain the measured information and, consequently, the error drops below the outlier detection threshold, the data is reincorporated into the information used in the reconstruction. The joint procedure is sketched in Figure 1.

Acquired Database and Sequence Design
The algorithm was tested on data from examinations on 517 babies (gestational ages at scan ranging from 32 to 45 weeks) who were all imaged in natural sleep using a 3 T Philips Achieva TX with a dedicated C ¼ 32-channel neonatal head coil (RAPID Biomedical) and subject handling system (25). T 2 -and inversion recovery T 1weighted multi-shot multi-slice fast spin echo sequences were acquired for each subject in both the axial and sagittal views. This structural information was gathered as part of a broader study where volumetric structural, resting-state functional and diffusion weighted MRI scans were also performed to study brain development within the developing Human Connectome Project (26).
Written informed consent for each participant was obtained from someone with parental responsibility prior to them being scanned. Main sequence parameters and number of completed acquisitions for each modality and orientation are reported in Table 1.
The shot and slice sampling structure of our sequences is sketched in Figure 2 for the axially acquired data (similar patterns are observed for the sagittal data, with differences only in the number of slices). Figure 2a shows the different interleaving patterns used for T 2 and T 1 , where a hierarchical interleave is used both in the phase encode and slice directions. The differences in the interleaving strategies for the T 2 and T 1 cases can be explained by their different targeted contrasts which are achieved by manipulating their T R and T E . To help interpretation, the echo and shot structure for k-space sampling are included in Figure 2b. In the T 2 case, the center of the spectrum is reached at the end of the echo train, whereas in the T 1 case it is acquired at the beginning.

Algorithmic Implementation
Certain term rearrangements, simplifications, and acceleration strategies are necessary to alleviate the computational complexity of the formulation proposed in Equation [2]: 1. For common Cartesian acquisitions, such as the fast spin echo sequences described before, the spectral measurements are arranged along lines, so the DFT along the readout direction can be separately precomputed. Assuming that the readout corresponds to the first dimension, i ¼ 1, we have: ðx;TÞ ¼ argmin FIG. 1. Block diagram of the alternating minimization approach with dynamic outlier rejection. where F ð2;3Þ is comprised of submatrices of size K 2 K 3 ÂN 2 Q corresponding to the 2D DFT with applied kspace resampling and multiband blipping pattern. 2. The datasets used to validate our method have been acquired using only in-plane acceleration, without simultaneous multi-slice excitations. In this case, the multiband factor is Q ¼ 1, the number of excitations to cover the FOV corresponds to the number of excited slices, R ¼ P, and the blipping pattern cycling period is also [10] where F ð2Þ is comprised of submatrices of size K 2 ÂN 2 corresponding to the 1D DFT with applied phase encode resampling. 3. A common setting in MRI reconstruction is that the number of reconstructed slices is matched with the number of acquired slices, so R ¼ P ¼ N 3 . For high resolution applications, another reasonable assumption is that the coil sensitivities remain approximately constant along the slice support, as they usually present a slow spatial variation. In consequence, they can commute with the slice profile and excitation matrices and be applied on a slice by slice basis: ðx;TÞ ¼ argmin where now H is comprised of N 1 N 2 SR Â N 1 N 2 SR submatrices of size N 3 Â N 3 ; A ðrÞ is comprised of N 1 N 2 SR ÂN 1 N 2 SR submatrices of size 1 Â N 3 that take the value 1 if r ¼ n 3 and are 0 otherwise. 4. As the slice profile filter H can be assumed to have approximately compact spatial support, the rigid transformation T can be applied to a slab in the surroundings of the excited slice that it affects. The number of slices in this slab, which we define to be 2V þ 1, should be given by the maximum amount of motion to be expected in the through-plane direction plus the spatial support of the slice profile. We can express the slab extraction operation by a matrix U and write: ðx;TÞ ¼ argmin where U is comprised of N 1 N 2 Â N 1 N 2 submatrices of size ð2V þ 1ÞR Â N 3 that take the value 1 if jn 3 À rj V and 0 otherwise, T U is comprised of SR Â R unitary submatrices of size N 1 N 2 ð2V þ 1Þ Â N 1 N 2 ð2V þ 1Þ; H U is comprised of submatrices of size ð2V þ 1Þ Â ð2V þ 1Þ, and A ðrÞ U is comprised of submatrices of size 1 Â ð2V þ 1Þ that take the value 1 if v ¼ V þ 1 (i.e., at the center of the slab) and 0 otherwise. We have set V ¼ 3 in our tests. 5. Coil information is compressed (27) to 95% of its energy to accelerate computation, which has been observed to have a minor impact in signal to noise ratio and ability to resolve motion. Thus, if we respectively denote by S and y the compressed coil sensitivities and measurements, we have: 6. A multi-resolution strategy is adopted that progressively incorporates high-frequency components into the formulation. Two resolution levels are used with an in-plane subsampling ratio of 2. The alternation between motion estimation and reconstruction is not applied at the finest scale to accelerate computations, which is observed to have a minor impact in the quality of obtained reconstructions. Colorbar shows the order in which the scanned information is acquired, so the images reflect the order in which the phase encode's (horizontal axis) and slices (vertical axis) are acquired during scan time. This acquisition structure is used in our method to define a set of motion states for which corresponding acquired information is assumed to be subject to negligible motion inconsistencies. b: Sampling structure in the phase encode direction. Colorbar shows the phase encode ordering, so the images reflect the order in which the different shots (horizontal axis) and echoes (vertical axis) covered the acquired k-space. This spectral acquisition structure is used by our method to infer motion estimates from the partial k-space information corresponding to each shot.

A graphical processing unit (GPU) version of the algorithm is used in practice.
Coil sensitivity maps are estimated from a separate reference scan (28). A spatial mask is obtained from the reference scan and used to constrain the reconstruction to be zero in the background. In the presence of strong motion, planning and masking have to be made conservative enough so as to prevent the imaged structures from moving outside the prescribed region of interest. Finally, outlier rejection threshold (s) and weight of reconstruction stabilizer (k) parameters have been tuned by visual inspection of exemplary volumes from within the pilot data acquired in our study and kept constant for all subjects. The former has been set to t ¼ 1:2 for all the modalities and orientations explored. The latter has been selected taking account of signal to noise ratio/resolution trade-off for each modality but kept the same for different orientations. Using this scheme, motion corrected reconstructions are obtained with a GEFORCE GTX TITAN X GPU in 30 min to 3 h per volume, with computation times mainly dependent on the level of motion.

Image Quality Assessment
Two metrics previously suggested as good indicators of image quality degradation in the presence of motion (15,29,30) have been used to quantitatively assess the relative image quality of uncorrected and corrected reconstructions in our cohort: The ' 1 -norm of a wavelet decomposition of the images (29), V a ¼ jjW 3 DbÀa xjj 1 , with W 3 DbÀa denoting the a-vanishing moments Daubechies (Db) wavelet decomposition at level 3 and a ¼ f1; . . . ; 4g. The gradient entropy (GE), V 5 ¼ HðjrxjÞ, of the images, which was shown to have the strongest correlation with observer quality scores in Ref. 30.
As we are interested in the motion correction ability, to isolate the effects coming from the thick slice profiles, we have compared the results of applying a conjugate gradient-SENSE reconstruction that incorporates our proposed deconvolution of the slice profiles without any motion or outlier modelling, non motion corrected with slice profiles (NMC-SP), and the fully motion corrected (MC) results. A paired right-tailed sign test against the null hypothesis that the median of the difference of these metrics V a ; a ¼ f1; . . . ; 5g, with and without motion correction is lower than zero or zero is performed to assess whether the motion corrected reconstruction effectively improves the sparsity of the wavelet coefficients and decreases the entropy of the reconstructed image gradient. In addition, we report the distribution of the relative change of the different metrics from corrected to uncorrected reconstructions, given by well as the improvement ratio r, the fraction of cases in which the corresponding metric decreased after correction, which would suggest a data quality improvement.

Visual Validation
To show the benefits that follow from including the main elements of the proposed framework, two levels of motion corruption and two different native acquisition planes are illustrated in Figures 3 and 4, which show respectively a mildly damaged T 2 -weighted sagittal acquisition and a grossly damaged T 1 -weighted axial acquisition. In each case, the separate rows depict transverse, sagittal, and coronal sections through the stacks of slices, as well as magnified areas to highlight fine detailed image differences, and the columns present different configurations of the motion corrected reconstruction with or without its main constituents: slice profile model, within-plane motion, through-plane motion and outlier rejection. In the mildly corrupted T 2 sagittal acquisition of Figure  3, the introduction of the slice profile in the reconstruction model (Fig. 3b) tends to suppress the saw-tooth boundary artifacts observed in non-native views (i.e., axial and coronal views) for standard SENSE reconstruction (Fig. 3a). However, localized artifacts are still observed in the sagittal view. The strength of these artifacts gets reduced when either through-plane (Fig. 3c) or within-plane (Fig. 3d) motion correction is introduced, as illustrated by the blue ellipse area. Finally, we observe how when both components of motion are incorporated into the model, the artifacts get further reduced and the through-plane consistency of the data gets improved (Fig. 3e). In this example, we have not included the results without outlier rejection as the inclusion of this feature did not influence the results.
In the highly corrupted T 1 axial acquisition of Figure 4, the reconstruction alternatives that do not include outlier rejection (Fig. 4a,b) are unable to recover consistent information, as the magnetization properties of certain shots can be severely altered due to motion. Without outlier rejection, application of motion correction may distribute magnetization artifacts in the slice direction, so the quality of motion corrected data (Fig. 4b) is not necessarily superior to that of non-corrected data (Fig. 4a). Although not included here, a similar effect is observed when comparing standard SENSE reconstructions and reconstructions using the slice profile model, where corrupted slices will damage neighboring slices. When outlier rejection is introduced (Fig. 4c,d,e), shots with irreconcilable data are effectively discarded and the reconstructed information appears more consistent. However, interslice misalignments are observed in the blue circle area when either within-plane (Fig. 4c) or through-plane (Fig. 4d) motion is not considered in the motion model, whereas the cortical boundaries appear better aligned when both components of motion are incorporated (Fig. 4e).
Further examples of the ability to correct a variety of artifacts in the acquired slices are shown in Figure 5

Quantitative Validation
The data quality of motion corrected and uncorrected reconstructions using the metrics described in "Image Quality Assessment" section is compared in Figures 7 (T 2 case) and 8 (T 1 case). Box plots of the distribution of the relative metric change from corrected to uncorrected reconstructions r, significance levels P, and improvement ratios r for each metric are included. Results for the whole dataset of 1883 brain volumes from 517 examinations are shown for different views and contrasts, axial T 2 (Fig. 7a), sagittal T 2 (Fig. 7b), axial T 1 (Fig. 8a), and sagittal T 1 (Fig. 8b). They correspond to different shot and slice sampling structures as shown in Figure 2. T 2 sequences are acquired at the beginning of the study whereas T 1 sequences are acquired at the end. Total scan duration is also different for T 2 and T 1 , as shown in Table 1. Moreover, different relative strengths of within-plane and through-plane motion may be expected for different orientations. In addition, a representative sample of reconstructions for different contrasts has been included as Supporting Information and made available in the online version of the article. Þ, selected as an indicator of the degree of motion correction for a given modality. On the one hand, this material is intended to show the agreement between the image quality metric r 4 and actual motion correction, with lower percentiles correlated with small level of correction (either due to inability to correct or due to small corruption). On the other hand, it is conceived to show unbiased information about the overall performance of the method in our database; aside from blindly choosing cases using the percentile rule, shown slices have also been picked out blindly at FOV/2 in the axial and coronal views and at FOV/3 in the sagittal view. Results for T 2 and T 1 modalities are respectively included in Supporting Figures S1-S10 and S11-S20.
Differences between corrected and uncorrected reconstructions in Figures 7 and 8 are highly significant (P ( 0:05) for all metrics, contrasts, and orientations. Thus, they show the ability of the motion correction method to effectively improve the compressibility and minimize the entropy of the gradient of reconstructed images, which, we interpret, is derived from its ability to reduce motion artifacts (see Supporting Information for visual verification). Larger improvement ratios (r) are generally obtained for T 1 than for T 2 scans. Larger improvement ratios can be due to larger overall motion artifacts in the sample or better performance of the method. We think that, in this case, we are in the first situation, as uncorrected T 1 -weighted images usually show more prominent motion artifacts than T 2 -weighted images (see Supporting Information). This may be due to the former being acquired at the end of the one hour examination protocol, with attendant risk that the unsedated neonates may stir, their longer scan duration as compared to the T 2 scans, and/or the specific inner time structure of this sequence, which involves the recovery from inversion. Improvement ratios are substantial for the wavelet measure and more moderate for the GE. We think this may be related with the wavelet metric showing a reduction when ghosts are corrected, which we have observed to be usual in our database, while the GE metric may be mainly an indicator of blurring, a less common artifact in the tested sequences.

Assembling Different Orientations
Despite being able to significantly diminish the spurious information introduced by motion inconsistencies, our method operates separately on each view, so the resulting volumetric image resolution may be biased by the acquired orientation, particularly in cases of large through-plane motion. To get an isotropic representation of the imaged volumes, we have deployed the slice to volume reconstruction method in Ref. 20 to assemble the information coming from different views. This method is able to correct for motion inconsistencies and intensity biases among slices and provides a super-resolved representation of the information coming from thick slices acquired in different orientations, so it matches very well with the residual corrections that may be required after our method. An example of the proposed pipeline for robust neonatal structural brain imaging is shown in Figure 9. More examples of view assembling are included as Supporting Information. In our project, resulting T 2 volumes are used for brain segmentation and T 1 volumes are later mapped to the T 2 space to study brain maturation.

DISCUSSION
Proposed method for motion corrected reconstructions in multi-shot multi-slice imaging is based on combining a strategy to estimate the rigid motion states of the imaged structure on a per-shot basis on top of a conjugate gradient-SENSE reconstruction that accounts for slice profiles to integrate the volumetric information of the acquisition. This article is mainly focused on the ability FIG. 8. Box plots of relative metric change r, P-values of a paired right-tailed sign test, and percentage of cases r in which the metric decreased for the ' 1 norm of Db wavelet decompositions and for the GE in motion corrected versus uncorrected reconstructions. Negative values in the paired box plots indicate a decrease in the corresponding metrics when applying motion correction, which has been documented as associated with an improvement in image quality (29,30). a: Axial T 1 . b: Sagittal T 1 .
FIG. 9. Slice to volume reconstruction-based assembling of motion corrected information of different views for suppression of residual artifacts and isotropic resolution: T 1 example. Volumetric data consistency is substantially improved for each of the views after applying our method, but residual motion may still be present due to remaining inconsistencies between slices, and non-native views may appear blurred as compared to native views. After slice to volume reconstruction correction, information is made consistent between views and a nearly isotropic representation of the imaged volume is obtained. In the bottom right corner, we show a magnified example comparing the results for the skull in the sagittal slice of the axial acquisition (left image, enclosed in blue) versus the corresponding results after view assembling (right image, enclosed in green), with residual motion inconsistencies strongly suppressed in the latter.
to estimate the within-plane and through-plane components of motion while discarding those shots that have intrinsic magnetization differences to get a consistent representation of the imaged volume. This has been tested for different structural modalities and orientations in a large database of neonatal studies obtained from the subjects in natural sleep, when there is a risk of sporadic movement. We have provided evidence of significant improvement both from visual assessment in selected examples and quantitative measures obtained from all available images. The formulation can be applied when using in-plane acceleration and simultaneous multi-slice excitation. However, only the former was tested and it is possible that some aspects of the implementation presented in "Algorithmic Implementation" section may need review for efficient reconstruction in the latter case.
Although we have demonstrated certain capability to correct through-plane motion, full correction at the prescribed resolution may not be possible because subject movement can result in insufficiently dense samples in some regions of the brain. Two complementary main strategies are envisaged to deal with nonuniform sampling. The first is the one adopted in this article: by acquiring additional sets of orthogonal of slices, inherent anisotropy can be corrected either as a post-processing step or by direct integration of multi-view information in the formulation. Optimal within-view acceleration and number of acquired orientations should be properly balanced for efficiency and robustness. The second would be the adoption of prospective motion correction strategies (31) with additional added value in preventing spinhistory effects. In this case, retrospective motion correction techniques as presented here may still be useful in providing correction of cross-calibration errors (32) or in treating artifacted shots.
Our method builds on a rigid model assumption that is only approximately valid for the brain. Thus, full correction cannot be expected for non-rigid sources of motion (33) or when the motion description is not unique due to other moving structures within the FOV (2). Issues due to non-rigid motion in the brain stem or in the neck can be alleviated by encoding using a foothead readout and limiting the motion estimates to the brain regions where the rigid body model fits well. This was actually applied in Ref. 15 but we have not found it to be necessary in the testing database used in this article, as the convergence of the algorithm is governed by the larger contrast of the cortical structures and their proximity to the coils. More important in practice was the application of a fat suppression technique to avoid motion inconsistencies due to water-fat shift. More complex motion models could be incorporated to our matrix formulation similarly to Ref. 14, however, the ability to effectively constrain and solve the problem will be largely application dependent, particularly when dealing with through-plane motion.
The hard criterion adopted in this work to differentiate between consistent and inconsistent data has been inspired by the type of motion most usually encountered in our application regime, where neonates usually move in sparse bursts during the acquisition, so that data is clearly degraded during those periods. However, this approach may present problems when generalized to other scenarios such as those with more homogeneous motion degradation during the acquisition (for instance for correcting small motions in high-resolution imaging). Using a weighted least squares technique along the lines of Ref. 34, may allow more effective data integration, but some work is needed on how to best design data uncertainty estimation and control the potentially unstable convergence of the joint motion estimation and reconstruction.
Other future developments could include more comprehensive treatment of (1) the inhomogeneous signal to noise ratio of parallel imaging as given by the array of coils and noise penalty implications of adopted acceleration methods, (2) the uncertainty in sensitivity map estimates, (3) the presence of other potential artifacts, and (4) the inhomogeneous resolution as given by throughplane motion and, perhaps, discarded shots. In addition, different priors can be used for motion estimation and reconstruction. Thus, a two-step approach could be envisaged. First, motion is estimated using a linear reconstruction technique as proposed in this article, where the prior is selected for the sake of stability and computability of the procedure. Second, a nonlinear reconstruction technique is used as a final step, where better approximants of the underlying image (35) could be used to integrate the sampled information in the presence of motion.

CONCLUSIONS
We have presented a method to obtain reconstructions for multi-shot multi-slice brain imaging including correction for 3D rigid body motion. The method is an extension of a previously proposed framework for motion estimation from incomplete parallel spectral information in purely volumetric MRI data (15). In this extension, we have provided a general formulation to jointly tackle the within-plane and through-plane motion problems as well as a strategy to iteratively solve for motion and get a motion-free reconstructed volume. The forward modeling framework adopted naturally incorporates parallel imaging acceleration methods (both in-plane or simultaneous multi-slice) and addresses through-slice movement by considering the slice excitation profiles. In addition, a procedure has been proposed to detect and discard artifacted shots. The method has been tested by reconstructing motion corrected images in a large cohort of neonates. Both visual validation and quantitative assessment have demonstrated significant benefits for different contrasts, different acquisition orientations and different shot and slice acquisition orderings. When combined with appropriate acquisition strategies including overlapped slices and multi-view information, our method has demonstrated to provide a robust motion-tolerant pipeline for some of the most prevalent structural brain MRI protocols. Thus, although our test data comes from a neonatal population, this proposal may also improve brain imaging performance in general, particularly for non-compliant populations.