Correcting inter-scan motion artefacts in quantitative R1 mapping at 7T

by   Yaël Balbastre, et al.

Purpose: Inter-scan motion is a substantial source of error in R_1 estimation, and can be expected to increase at 7T where B_1 fields are more inhomogeneous. The established correction scheme does not translate to 7T since it requires a body coil reference. Here we introduce two alternatives that outperform the established method. Since they compute relative sensitivities they do not require body coil images. Theory: The proposed methods use coil-combined magnitude images to obtain the relative coil sensitivities. The first method efficiently computes the relative sensitivities via a simple ratio; the second by fitting a more sophisticated generative model. Methods: R_1 maps were computed using the variable flip angle (VFA) approach. Multiple datasets were acquired at 3T and 7T, with and without motion between the acquisition of the VFA volumes. R_1 maps were constructed without correction, with the proposed corrections, and (at 3T) with the previously established correction scheme. Results: At 3T, the proposed methods outperform the baseline method. Inter-scan motion artefacts were also reduced at 7T. However, reproducibility only converged on that of the no motion condition if position-specific transmit field effects were also incorporated. Conclusion: The proposed methods simplify inter-scan motion correction of R_1 maps and are applicable at both 3T and 7T, where a body coil is typically not available. The open-source code for all methods is made publicly available.


page 4

page 6

page 7


MoCoNet: Motion Correction in 3D MPRAGE images using a Convolutional Neural Network approach

Purpose: The suppression of motion artefacts from MR images is a challen...

Clinically Translatable Direct Patlak Reconstruction from Dynamic PET with Motion Correction Using Convolutional Neural Network

Patlak model is widely used in 18F-FDG dynamic positron emission tomogra...

Signal to Noise and b-value Analysis for Optimal Intra-Voxel Incoherent Motion Imaging in the Brain

Intravoxel incoherent motion (IVIM) is a method that can provide quantit...

FlowNet-PET: Unsupervised Learning to Perform Respiratory Motion Correction in PET Imaging

To correct for breathing motion in PET imaging, an interpretable and uns...

1 Introduction

Quantitative MRI, and the push towards in vivo histology, aims to extract tissue-specific parameters from a series of weighted volumes [weiskopf2021quantitative]. For example, the longitudinal relaxation rate, , which is sensitive to important biological features, such as myelin and iron content, can be quantified with the variable flip angle (VFA) approach, e.g. deoni2003rapid, helms2008quantitative. A common assumption when computing quantitative metrics is that certain multiplicative factors, such as the signal intensity modulation imposed by the receiver coil’s net sensitivity profile, are constant across the weighted volumes. However, this is invalid if motion occurs between the volume acquisitions. In the case of neuroimaging, rigid body co-registration can be used to realign the brain but will not correct for the differential coil sensitivity modulation, which in maps computed with the VFA approach can lead to mean absolute error approaching 20% [papp2016correction].

A correction scheme has previously been proposed by papp2016correction and validated for mapping at 3T. The position-specific net receive sensitivity is estimated from two rapid low-resolution magnitude images, received on the body and array coils respectively prior to each VFA acquisition. The more homogeneous profile of the body coil is used as a reference to compute the net receiver sensitivity, which is then removed from the VFA acquisitions. This approach effectively assumes that the body coil’s modulation is consistent across volumes instead of that of the array coil. This in itself is a potential limitation, as is the general unavailability at body coils at higher field strengths.

Here we propose an alternative whereby we estimate the relative sensitivity between volumes. This approach does not fully remove the receiver’s sensitivity modulation but does remove the bias that differential modulation introduces in quantitative metrics. Only the calibration images obtained with the array coil are required, i.e. less data than the originally proposed method [papp2016correction]. To validate the approach, we focus on maps computed with the multi-parameter mapping (MPM) protocol [weiskopf2013quantitative]. We first compare performance with the established method of Papp et al. at 3T [papp2016correction] and then demonstrate a reduction of inter-scan motion artefacts at 7T under a range of different motion conditions. We further demonstrate that, unlike at 3T, the transmit field also exhibits substantial position-specific variability at 7T. As a result, the most precise estimates were obtained by accounting for both position-specific and effects.

2 Methods

2.1 Theory

mapping can be achieved by acquiring spoiled gradient echo volumes with at least two different flip angles [deoni2003rapid, helms2008quantitative]. At a given spatial location, , the image intensity, , for a given nominal flip angle is:


where is the receive sensitivity, is the proton density, is the transmit field, is the longitudinal relaxation rate, TR is the repetition time, and indexes the VFA acquisition. Co-registration allows for inter-scan motion by realigning anatomical structure across acquisitions. Under the small flip angle approximation [helms2008quantitative], with two nominal flip angles ( = {1,2}), can be computed as follows:


Typically, it is assumed that and the sensitivities simplify out. However, this assumption is invalid if inter-scan motion has occurred leading to substantial bias in [papp2016correction]. This can be avoided by accounting for the relative sensitivity across positions: . Substitution for in equation (2) gives:


The method of papp2016correction did not include the relative sensitivity but referenced to an additional calibration image acquired on the body coil, assuming that the body coil modulation was position-independent.

It is commonly assumed that the transmit field is sufficiently smooth as to be considered position-independent, i.e. , such that:


However, in this work we show that this assumption does not hold at 7T, and that incorporating position-specific transmit field estimates maximises the precision of .

Figure 1: 7T Example. The acquired calibration images, have different orientation due to participant movement between acquisitions. Co-registration can align the images spatially, but does not correct for their differential sensitivity field modulation, visible via their ratio, / = . Alternatively, the joint log-likelihood of a generative forward model that embeds the spatial transformation from a mean image, , to native space can be maximised to determine the mean image and modulating sensitivities, , that best explain the acquired images . The generative modelling approach produces a similar relative modulation (/) but allows the sensitivity modulations to be explicitly estimated, and for the corrected images to have the minimal modulation of the mean image.

Ratio approach The calibration data used to correct inter-scan motion artefacts comprised rapid low resolution, coil-combined magnitude images acquired immediately prior to each high resolution VFA acquisition. These images, , assumed to have been rigidly co-registered and resliced to a grid of voxels, can be written as the product of a common image and a net sensitivity field . The relative sensitivity, can be computed with respect to one of the calibration acquisitions, used as a reference:


Dividing each VFA acquisition by its relative sensitivity results in a common modulation, , which, although less homogeneous than the body coil used by Papp et al., more faithfully restores the validity of assuming common modulation when computing

. This ratio approach risks noise amplification, particularly in regions of low signal-to-noise ratio (SNR). This is combatted by isotropically smoothing

and before taking their ratio.

Generative approach A potentially more robust alternative is to cast the computation of the relative coil sensitivities, and a common image modulated by them, as an inference problem in a probabilistic generative model of that incorporates noise and can also embed knowledge about the spatial smoothness of the sensitivities. Full details are in Appendix A.

2.2 Experiment

Participants One participant (F, 31 years) was scanned at 3T (MAGNETOM Prisma, Siemens, Erlangen, Germany) using a body coil for transmission and either the body coil or a 32-channel head array coil for reception. Three additional participants (2F, 1M; 32 - 41 years) were scanned at 7T (MAGNETOM Terra, Siemens, Erlangen, Germany) using an 8-channel transmit, 32-channel receive head array coil (Nova Medical, Wilmington, MA, USA) in TrueForm mode. All data were acquired with approval from the UCL research ethics committee.

MPM Datasets MPM data were acquired using a multi-echo spoiled gradient echo sequence with flip angles of 6 (PD-weighted, ”PDw”) and 26 (T1-weighted, ”T1w”), a TR of 19.5 ms and an RF spoiling increment of 117

with a total dephasing gradient moment per TR of 6

. Eight echoes were acquired with TE ranging from 2.56 ms to 15.02 ms in steps of 1.78 ms using a bandwidth of 651 Hz/pixel. Data were acquired with 1mm isotropic resolution over a field of view of 160 mm right-left and 192 mm in the anterior-posterior and superior-inferior directions. Elliptical sampling and partial Fourier, with factor 6/8 in each phase-encoded direction, were used to accelerate the acquisition, leading to a scan time of 5 minutes per volume. A map was estimated by acquiring a series of spin and stimulated echoes using previously described 3T and 7T protocols [lutti2010optimization, lutti2012robust].

For inter-scan motion correction, additional single echo acquisitions were acquired prior to each VFA acquisition with a flip angle of 6, TE = 2.4 ms, TR = 6.5 ms, a bandwidth of 488 Hz/pixel and no acceleration schemes. At 3T, these data were acquired, receiving sequentially on the array and body coils, with 8 mm isotropic resolution leading to a scan time of 6 s per volume. To capture the greater spatial variation in the net sensitivity field at 7T, the resolution was increased to 4 mm isotropic leading to a scan times of 18 s per volume, but acquired only on the array coil due to the absence of a body coil.

Motion Conditions Two MPM datasets were acquired to define baseline reproducibility. Participants were then instructed to move to a new, arbitrary position within the confines of the coil. A localiser was acquired and the field of view repositioned as necessary to ensure appropriate brain coverage in the new position. A third MPM data set was then acquired.

Analysis maps were computed using the hMRI toolbox [tabelow2019hmri], which uses the small flip angle approximation [helms2008quantitative] and corrects for imperfect spoiling [preibisch2009influence], which always used the map acquired in the space of the PDw acquisition. Maps were computed with and without inter-scan motion using all possible PDw and T1w combinations. To ease comparisons, all maps were constructed in the space of the first PDw volume. Rigid transformations between all volumes and the first PDw volume were estimated using SPM12 (Wellcome Centre for Human Neuroimaging) having first corrected for intensity non-uniformity and skull-stripped the images. maps were computed with and without the proposed inter-scan motion correction schemes. At 3T, maps were also computed with the method of papp2016correction. An isotropic kernel of 12mm full-width-at-half-maximum – the default in the hMRI toolbox – was used to smooth the calibration images prior to computing the relative (proposed) and absolute (Papp et al.) sensitivities. The generative modelling approaches used , , and 15 iterations.

per contrast In each case, two different corrections for transmit field inhomogeneity were performed: the first assumed the transmit field was identical across head positions, and used that of the PD-weighted acquisition; the second used position-specific maps. In the latter case effective flip angle maps were computed for each volume and used in equation (3).

Error metric maps with no inter-scan motion, and no corrections, were averaged to produce a ‘ground-truth’ map, . This set of

maps was also segmented to create a mask selecting those voxels with a mean probability of being in WM, GM or CSF greater than 50%. For participant 3, the cerebellum was excluded, using the SUIT toolbox

[DIEDRICHSEN2006127] in SPM, as a result of mapping failure caused by excessively large off-resonance. For the voxels within the resulting participant-specific mask, the mean absolute error, MAE, for each map was computed with respect to the ‘ground-truth’ map, as:


These errors are reported as percentages.

3 Results

Exemplar images, relative sensitivities, and results from the generative modelling are shown in figure 1. The and error maps obtained at 3T and 7T are shown in figures 2 and 3

respectively. The means and standard deviations of the MAE are reported in Table

1. The differential impact of correcting for transmit and receive field effects is illustrated in figure 4. This shows and error maps without motion, and with motion having implemented (i) no correction, (ii) correction only for receive field effects, , (iii) only for transmit field effects, , or (iv) for both effects in combination.

Figure 2: Results at 3T. The first row shows example maps constructed with each method. The second row shows (normalised) error maps with respect to the ground-truth map . The third row shows histograms (filled area) of

within the GM (green) and WM (purple), and their log-Normal fit (solid line); these histograms display probability distributions and therefore integrate to 1.

3.1 3T Validation

When the T1w and PDw volumes were acquired in the same position, the MAE captured the test-retest variability, which was approximately 3% at 3T and 4-5% at 7T. In the absence of overt motion, correcting for the differential sensitivity modulation did not substantially change the MAE. In the presence of overt motion, the MAE rose to 10%. It was reduced to 4.7% by the method of Papp et al. and to less than 4.4% by the proposed correction schemes, with or without incorporating the body coil in the generative modelling approach. The histograms in figure 2 confirm that the method did not introduce any bias to the estimates.

3.2 Extension to 7T

At 7T, the range of motion varied across participants. The net (root-mean-square across axes) translation ranged from 1.6 to 7.9 mm, while the net rotation ranged from 2.7 to 11.2 degrees. The overall amplitude of motion dictated the increase in MAE, which reached a maximum of 13.4% under the tested conditions (c.f. motion summaries in figure 3 and MAE in table 1). The proposed correction scheme reduced the MAE (5-8%), though not to the level of no overt motion.

The variability of the transmit field, , across head positions was found to be much higher than at 3T. Incorporating position-specific maps reduced the MAE (5-12%) even without correcting for the differential receive sensitivity modulation.

The greatest reductions in MAE were obtained by correcting for position-specific transmit and receive fields, reaching 4-7%, converging on the level obtained in the absence of overt motion (i.e. 4-5%).

3.3 Comparison of Methods

Overall, the ratio and generative modelling approaches to estimating the relative sensitivities performed similarly. The generative modelling approach performed marginally better at 3T (0.1% lower MAE), while the simpler ratio approach performed better at 7T (up to 0.5% lower MAE).

rcc — cccccc & Dataset & Motion & No correction & Ratio & Generative & Generative (array + body) & Papp et al.
3T & #1 & No & 3.0 0.1 & 3.0 0.1 & 3.0 0.1 & 3.1 0.2 & 3.7 0.2
& & Yes & 10.1 0.8 & 4.4 0.4 & 4.3 0.3 & 4.3 0.1 & 4.7 0.3
& Dataset & Motion & No correction & Ratio & Generative & B per contrast & Ratio & B per contrast & Generative & B per contrast
& #1 & No & 4.0 0.3 & 4.1 0.4 & 4.1 0.3 & 4.1 0.2 & 4.1 0.4 & 4.1 0.3
& & Yes & 5.8 0.7 & 4.8 0.3 & 4.9 0.4 & 4.9 0.3 & 3.8 0.3 & 3.9 0.3
7T & #2 & No & 4.8 0.2 & 4.9 0.3 & 4.9 0.3 & 4.8 0.2 & 5.0 0.2 & 5.5 0.1
 & & Yes & 8.4 0.4 & 6.1 0.3 & 6.2 0.3 & 8.0 0.4 & 5.3 0.1 & 5.5 0.1

& #3 & No & 5.1 0.2 & 5.0 0.2 & 5.0 0.2 & 5.0 0.2 & 5.0 0.3 & 5.0 0.2
 & & Yes & 13.4 1.9 & 7.9 1.1 & 8.3 1.2 & 11.8 0.9 & 6.3 0.3 & 6.8 0.3

Table 1: MAE (mean s.d. across repeats, in %) with respect to the average reference .
Figure 3: Results at 7T. The first column shows the ground-truth map for the three participants sorted based on the magnitude of inter-scan motion. Uncorrected and corrected maps are shown in the middle with the corresponding (normalised) error maps with respect to the ground-truth on the right. Correction is only applied for net receive sensitivity modulation, and not for . Rows 1 to 3 of the figure correspond to datasets 1, 2 and 3 as reported in Table 1.
Figure 4: Combining and correction at 7T, for participant 3 (last row in figure 3). The first row shows an example map without and with inter-scan motion, before and after net receive sensitivity correction, and employing a separate map for each contrast in the case of inter-scan motion. The second row shows (normalised) error maps with respect to the ground-truth map . In this example, inter-scan motion correction was performed with the generative modelling approach.

4 Discussion

We have introduced methods for correcting inter-scan motion artefacts in quantitative MRI that do not rely on the availability or spatial homogeneity of a body coil. The approaches are based on estimating the relative sensitivity modulation across positions, and successfully reduced error in maps at both 3T and 7T.

At 3T, the proposed approaches outperformed a previously established correction method [papp2016correction]. This can be attributed to the fact that the method of Papp et al. assumes that the reference modulation of the body coil is independent of position, whereas the proposed methods do not. Instead they specifically account for the relative sensitivity across positions thereby restoring consistent modulations.

In the motion conditions tested here, both proposed approaches (ratio with Gaussian smoothing, or generative modelling) produced comparable improvements in reproducibility in the presence of inter-scan motion. Equally importantly, when there was no overt motion neither method decreased reproducibility, which was at a level in keeping with previous reports for similar resolution MPM data [leutritz2020multiparameter, weiskopf2013quantitative].

The simpler ratio method defines one calibration image as the reference (denominator in equation (5)), and may therefore be vulnerable to low SNR. The alternative generative modelling approach inherently adapts to variable SNR by estimating the position-specific net sensitivity modulation relative to a common image, which is their barycentre mean. This common image dictates the final modulation of all the corrected volumes. The generative model can also easily incorporate any additional data, e.g. body coil images as done at 3T, which flattens the final modulation. Furthermore, rigid registration could be interleaved with model fitting [ashburner2013symmetric] to reach a better global optimum. Finally, the generative model could naturally be integrated with any fitting approach that defines a joint probability over all acquired data, such as balbastre2021model in the context of MPM.

The impact of movement on the effective transmit field has previously been investigated in the context of specific absorption rate management [kopanoglu2020specific, le2017probabilistic, wolf2013sar, shajan201416]. An important additional finding of the present work is the impact this can have on estimates at 7T, which was negligible at 3T as demonstrated previously [papp2016correction]. Although errors were more substantially reduced by correcting for receive field effects, they could only be reduced to the approximate level of no overt motion by additionally accounting for the positional-dependence of the transmit field. However, acquiring a map at multiple positions comes at a cost of increased scan time.

5 Conclusions

Inter-scan motion causes serially acquired weighted volumes to be differentially modulated by position-specific coil sensitivities leading to substantial errors when they are combined to compute quantitative metrics. We have demonstrated the efficacy of two methods at reducing these artefacts in the context of mapping. The proposed methods do not require a body coil making them ideally suited for use at 7T, and can be extended to the computation of other quantitative metrics, such as magnetisation transfer saturation [helms2008high], that similarly assume constant modulation across multiple weighted acquisitions.


Appendix A Generative Model of Sensitivity-Modulated Images

The proposed generative model is similar to that of ashburner2013symmetric (section 2.3), without non-linear deformations. The magnitude images of the calibration dataset, , can be written as the voxel-wise product of a mean image and the net sensitivity field

plus additive noise, approximated as Gaussian with variance

, which is assumed to be uncorrelated across . The net sensitivity fields can be written as diagonal matrices , such that the corresponding conditional probability is:


If the images have been co-registered but not resliced, a mean space can be defined by computing the barycentre of all aligned orientation matrices [ashburner2013symmetric]. The linear operation of resampling an image from mean to acquired space can be encoded by the matrix such that the conditional probability becomes:


While this is the approach taken in practice, the following derivation is restricted to the resliced case for clarity.

The magnitude of each sensitivity field is unknown since the intensity scaling depends on many parameters. However, the sensitivities are known to vary smoothly in space, which is captured by a probability distribution that penalises the field’s bending energy [ashburner2007fast], which integrates the squared curvature of the sensitivity field, making it invariant to intensity scaling. The net sensitivities are encoded by their logs, (i.e., ), such that invariance under shifts in log-space implies invariance under scales in exponentiated space, which is also incorporated into the prior. In a discrete setting, computing the bending energy reduces to the quadratic term

; The prior distribution over sensitivities is therefore defined as a Normal distribution over their logs:


where is an image-specific regularisation factor.

The joint model likelihood, or its negative log, is obtained by combining the likelihood in (7) and the prior in (9): . This is minimised with respect to and . Neglecting terms that do not depend on these variables, yields the objective function:


Differentiating with respect to the mean image, , while keeping the sensitivities fixed gives a voxel-wise () closed-form update:


The log-sensitivities have no closed-form solution necessitating an iterative method. The objective function is not everywhere convex, but the likelihood term resembles that of a nonlinear least-squares problem, which can be solved using Gauss-Newton optimisation. Gauss-Newton is a modification of Newton-Raphson that uses Fisher’s method of scoring, which amounts to replacing the Hessian at any point with its value at the optimum. With and , the gradient and Hessian used in the Newton-Raphson iteration are:


This ensures that the Hessian used in the Newton-Raphson iteration is positive definite, but does not ensure that the iteration monotonically improves the objective function. We therefore replace the Gauss-Newton Hessian with the more robust preconditioner:


which has been shown to yield monotonic convergence [balbastre2021model]. The inversion in equation (15) is performed with a full multi-grid solver that leverages the sparsity and structure of the preconditioner [ashburner2007fast].

Finally, a global scaling field , applied to both the sensitivities () and mean image (), ensures that the product is unchanged. Keeping only terms of the objective function that depend on this scaling field gives:


By differentiating, the optimal scaling field is:


Therefore, at the optimum, the (weighted) mean log-sensitivity field must be zero, and the mean image is a barycentre of the calibration images. To accelerate convergence, this condition is enforced after each global iteration.

In general, the variability of the sensitivity modulation smoothness across images is unknown. Therefore all are given the same value . However, images acquired with the body coil will have a flatter sensitivity field modulation, which can be incorporated by setting .

Code availability: The code used to fit the generative model is available at A modified version of the hMRI toolbox that integrates this approach and enables correction on a per-contrast basis is available at The ratio approach can be performed natively with the hMRI toolbox: