Multiple motion encoding in Phase-Contrast MRI: A general theory and application to elastography imaging

While MRI allows to encode the motion of tissue in the magnetization’s phase, it remains yet a challenge to obtain high ﬁdelity motion images due to wraps in the phase for high encoding eﬃciencies. Therefore, we propose an optimal multiple motion encoding method (OMME) and exemplify it in Magnetic Resonance Elastography (MRE) data


Introduction
Phase-contrast Magnetic Resonance Imaging (PC-MRI) is a wellestablished method for measuring flow velocities ( Srichai et al., 2009;Markl et al., 2012;Dyverfeldt et al., 2015 ) or tissue displacements due to harmonic excitation as used in Magnetic Resonance Elastography (MRE) ( Muthupillai et al., 1995;Mariappan et al., 2010;Klatt et al., 2007;Singh et al., 2015;Pepin et al., 2015;Dong et al., 2018;Manduca et al., 2021 ).MRE is used for the non-invasive characterization of the mechani-wraps (abrupt jumps of 2 π k , with k ∈ Z ) occur if the encoded phase exceeds those limits.Consequently, for a given encoding efficiency, there is a fixed amplitude range or dynamic range , where motion can be acquired without phase wraps.In other words, if the encoding efficiency is too large and therefore the true motion amplitude exceeds the dynamic range, phase wraps occur.Unfortunately, selecting a large dynamic range leads to poor quality images since -for a given signal-to-noise-ratio (SNR) in the magnitude image -the "motion-to-noise-ratio" (MNR) is inversely proportional to the dynamic range.Moreover, acquiring and averaging several images with fixed but high dynamic range (i.e. with no wraps) decreases the standard deviation of the averaged image with order O(d G N −1 / 2 G ) , with N G the number of measurements for a fixed dynamic range d G .
Therefore, it is a common practice to use low dynamic ranges and then to remove the wraps afterwards.This allows a faster decrease of the standard deviation of the estimated motion encoded image than averaging.There are usually two type of approaches.
First, unwrapping algorithms have been developed by assuming that the motion field is smooth in space, see e.g.Barnhill et al. (2015) , Loecher et al. (2016) and references therein.Nevertheless, they cannot recover the true underlying motion and eventually fail when the aliased regions are highly heterogeneous, subject to noise or include nested wraps (i.e. when | k | > 1 ).In such cases, the unwrapped phase appears to be distorted and obstructs further data processing steps which leads to artifacts in the estimates of tissue properties ( Manduca et al., 2021 ).For instance, the simple 2 π -unwrapping Flynn (1997) algorithm is inherently two-dimensional and fails to unwrap noisy complex wraps when no well-defined enclosed region exists.The true motion cannot be recovered because arbitrary 2 π -offsets are added.Gradient based algorithms ( Sack et al., 2008 ) only yield the derivative of the phase and amplify noise.Laplacian based unwrapping algorithms ( Schofield and Zhu, 2003 ) remove the constant and linear terms from the data and induce spatial smoothing, altering the resulting phase where important details may be lost.
Second, voxelwise motion reconstructions using the so called dual motion encoding strategies have been proposed in PC-MRI which are based on unwrapping low dynamic-range data by exploiting high-dynamic range data ( Lee et al., 1995;Schnell et al., 2017;Carrillo et al., 2019;Yin et al., 2018 ).In other words, measurements with a reduced dynamic range (hence, improved MNR) are unwrapped using a measurement with a larger dynamic range.Those methods are performed at each voxel independently and therefore they do not assume or enforce smoothness of the motion-encoded phase field.This allows the reconstruction of the correct motion images, but, at the cost of additional measurements.However, dual motion encoding also fails in the presence of noise when the MEG amplitudes do not differ largely.
Hence, the aim of this work is threefold.Firstly, we analyze dual motion-encoding strategies showing that dual motion encoding methods are limited to low noise phase images.
Therefore, we then develop a mathematical framework for multiple motion encoding, henceforth termed Optimal Multiple Motion Encoding (OMME) as an extension of optimal dual motion encoding reported in Carrillo et al. (2019) .We show that OMME outperforms dual motion encoding in terms of unwrapping's robustness to noise.This improvement comes at the cost of additional scans, however being considerably more cost effective in terms of scan time than standard image averaging.
Finally, we propose a MRE scan protocol for OMME and test it on in vivo brain data.We first confirm that OMME provides either a more noise robust unwrapping with similar MNR or improved MNR with similar noise robustness compared to dual motion encoding strategies.Moreover overcoming the limitations of unwrapping algorithms not only increases MNR but also allows to recover more detailed stiffness maps than using standard unwrapping methods.

Theory
In this section, we first introduce the notation and the mathematical model of phase and motion, including theoretical statistical analysis of dual encoding methods.Then, we introduce OMME and derive its statistical properties.Please notice that the theory in this section is formulated for one voxel, one motion encoding direction and one time instance.This means, at the same time, that OMME does not rely on spatio(-temporal) regularity assumptions as conventional unwrapping methods.Therefore, in practice OMME is applied to each direction and spatio-temporal location independently.

Single motion encoding
For a given MEG "G ", the model of measured phase can be written in the form with the following notation: • u denotes the motion of the media (quantity of interest) in the direction of the MEG (dimensions are velocity or displacement, depending on the quantity encoded); • d G is the "dynamic range" of motion encoding, which depends on the MEG's amplitude, duration, shape and assumption on the motion (same dimensions as u ); • δ u G is the spatially varying background phase, which in general depends on gradient imperfections of spin-echo sequences, the MEG (e.g by eddy currents and Maxwell effects) and the motion since imaging gradients also encode motion (time dependent); • ε is a zero mean Gaussian random variable representing the measurement noise in the phase.Its variance depends on the SNR of the magnitude measurements.
The background phase can further be modelled as follows: with the following notation: • ϕ 0 is the time constant (static) background phase of the imaging gradients due to gradient imperfections and concomitant fields • δ G is a time constant MEG-dependent phase induced e.g by eddy currents and Maxwell effects • m (u ) is the motion-dependent phase encoded by the imaging gradients.
In order to separate the unwanted background phase from the desired dynamic (time dependent) contributions of the MEG, it is necessary to perform a series of acquisitions with and without applied vibrations and MEGs.In flow velocity-encoded MRI, the measurement of δ u G is the common practice (i.e.four points 4D Flow) ( Markl et al., 2012 ).In MRE, we will show in the methods section that δ u G can be measured considerably faster compared to ϕ G (u ) .
We define as single encoding the situation where u estimated from a measurement of the background phase δ u G and the full phase ϕ G (u ) , for a single MEG G .
The constant C δ represents the proportion in standard deviations between the motion-encoded phase and the background phase measurements, which may differ in practice.Consequently, u G ∼ , with u true the true motion.Therefore, for a fixed value of σ ϕ , d G should be chosen as small as possible in order to minimize Var (u G ) .However, phase can be measured only within the interval For a given dynamic range d G , a possible approach to reduce Var (u G ) is to average more measurements, say n G times.Then: From this relation it is clear that decreasing d G is more effective than increasing the number of measurements n , since Var (u G ) de-

Phase-contrast from two MEGs
We assume now that we measure phases with two differ- Without loss of generality, we assume 0 < β < 1 , obtaining dynamic ranges d 1 and d 2 = βd 1 , respectively.This results in four measured phases ϕ 1 , δ u 1 , ϕ 2 , δ u 2 .We assume that these values already include the additive noise as indicated above.We also assume for the subsequent computation of the variances that these phases have no wraps.
From the four phase measurements, four motion images can then be estimated: (5) The indexes pc and ps denote phase contrast and phase sum , respectively.Notice that , the variances of the different motion estimators satisfy:

Classical dual motion encoding unwrapping
As explained above, when the dynamic ranges are decreased, wraps appear.Dual encoding reconstructions aim to unwrap a motion u low obtained with low dynamic range d low using a motion u high obtained with a high dynamic range d high as follows Lee et al. (1995) : (11) with N.I. the nearest integer operator.This leads to Var (u uw ) = Var (u low ) when the unwrapping is successful.This method will be denoted in as standard dual encoding .We denote as effective dynamic range of the dual motion encoding unwrapping method d e f f = d high .
To choose u low and u high we apply the following reasoning.Firstly, we select u low = u 2 since it possess a higher dynamic range than u ps (e.g.d pc = 3 / 2 d ps for β = 1 / 2 ) only with a slightly higher variance (e.g.Var (u 2 ) / Var (u ps ) = 1 .125 for β = 1 / 2 ).Secondly, we select since u pc has the desired effective range Then, the lowest dynamic range is . This leads to both approaches be equivalent if β = 1 + (k − 1) / (k + 1) .In Yin et al. (2018) , it was for instance taken k = 0 .777 hence equivalent to a β = 0 .875 .Note that we excludes the case k = 1 from the analysis since it does not correspond formally to dual encoding but to single encoding because d ps → ∞ .

Optimal dual encoding unwrapping
In Carrillo et al. (2019) , a new method for unwrapping two motion-encoded images was introduced, Optimal Dual Venc (ODV).The method is based on the formulation of the phase contrast problem as the minimization of cost functional.For the single motion encoding case, the cost functional has the form: which comes from a least squares approximation for the angle by measuring the components of a vector.
It is easy to see that the period of J i (u ) is 2 d i , and therefore local minimum among with smallest value (in absolute terms) , corresponds to the single encoding phase-contrast motion.
For the dual encoding case, the problem shifts from finding the local minima of J i (u ) to find the global minima of It was proven in Carrillo et al. (2019) where the last equality follows from the proof in Appendix A.1 .
Therefore, in these scenarios both standard dual encoding and ODV methods have the same effective dynamic range and they can be fairly compared to each other.Another contribution of this work is the computation of the variance of the ODV estimate u * leads to The computation is detailed in Appendix A.2 for the more general case of multiple encoding.

Limitations of dual encoding
From Eq. ( 7) , notice that when β → 1 − (left limit), then d e f f → ∞ .Therefore, one may think that d 1 , d 2 could be chosen arbitrarily small to minimize Var (u 1 ) , Var (u 2 ) while keeping d e f f large.However, we will show here that noise affects the unwrapping performance of the methods.Therefore, dual motion encoding strategies have limitations which become more important the closer d 1 and d 2 are.Fig. 1 presents the previous findings in a graphical way.There, we show the standard deviations of the estimators (i.e. the square root of the variances) versus the effective dynamic ranges for various values of β.Each sub-figure was generated for a given value of β ∈ { 1 / 4 , 1 / 2 , 2 / 3 , 3 / 4 , 1 } and σ φ = { 0 .01 , 0 .05 } , by the following procedure: • Ground truth values are set as: u true = 1 and δ u 1 = δ u 2 = 0 .9 π .• For a fixed value of d e f f and β, d 1 are computed as: • Measurements of motion encoded phases ϕ 1 , ϕ 2 were generated using Eq. ( 1) , the ground truth values of the parameters defined above and adding Gaussian noise with standard deviation σ ϕ .Measurements were wrapped to the interval 2 are perturbed adding Gaussian noise with standard deviation σ ϕ .
• Similarly u e f f was generated being the single motion phase contrast estimate with d G = d e f f is computed for comparison.• The standard deviation of such estimates considering the 50 0 0 realizations is computed.
The quality of the results depends on both values of β and σ ϕ .For small values of σ ϕ , the empirical and theoretical standard deviations match, but the empirical increases -deviating from the theoretical one -when d e f f → | u true | , as expected, due to the wraps.This issue becomes more relevant when σ ϕ grows.In this low noise scenario, the maximum gain with respect to the case of repeating the same measurements (i.e.β = 1 ) is when 3 and d 2 = d e f f / 4 .However, the reconstruction with β = 3 / 4 becomes unstable when increasing σ ϕ .A similar effect can be observed with β = 1 / 2 .The most robust variant with respect to noise for both standard and optimal methods appears to be β = 1 / 2 , where d 1 = d e f f and d 2 = d e f f / 2 .In case of the optimal method, this can be explained by the fact that the local minima of both J 1 and J 2 cost functionals have maximal distance.For the other values of β, this distance is much smaller, we refer to Carrillo et al. (2019) for details.
In particular for β = 1 / 2 , between both methods, the optimal dual encoding appears to be more robust with respect to noise, especially when d e f f → | u true | , and slightly better than the standard dual encoding approach when d e f f > | u true | due to (A.10) .The possible explanation is that unwrapping and noise compensation are done simultaneously, and therefore, a more robust unwrapping method results.For other values of β, it appears that the standard dual approach performs better in terms of robustness when

Optimal multiple motion encoding (OMME)
We now propose a systematic method to include N measurements with dynamic ranges d 1 , . . ., d N in order to extend the effective dynamic range d e f f while keeping Var (u * ) as the one for the smallest d 1 , . . ., d N .The idea is that by doing so, we can increase the robustness of the unwrapping when σ ϕ increases.Therefore, such strategy can be of great utility when high quality images are needed e.g. at high spatial resolutions, for instance.
The optimal dual encoding formulation allows a straightforward extension to multiple MEGs, i.e.
The multiple motion encoding reconstruction u * is then the global minimum of smallest magnitude within [ −d e f f , d e f f ] , with d e f f the dynamic range of OMME.From the proof in Appendix A, J N has periodicity equal to the least common multiplier of 2 The variance of u * , is given by see detailed computation in Appendix A.2 .For instance for the case of β = 1 / 2 -which is the most robust as it was shown above -the respectively.Therefore, the gain in noise reduction with respect to the lowest dynamic range is only slightly reduced.At the same time, the computational complexity of the an exhaustive search of the global minimum of J N (u ) increases considerably with N, since the interval [ −d e f f , d e f f ] needs to be sampled according to d N .
Therefore, we propose here to just use J N (u ) to guide the unwrapping of u N , i.e. to find u * by solving: and then to set For image datasets as used in this work, e.g. for N = 3 solving Problem ( 21) is about 9 times faster than an exhaustive global minimum search of J N (u ) .
Fig. 2 shows the results OMME using a number of measurement combinations, values of β and σ φ .For β = 1 , the fast version of  OMME was used, while for β = 1 the standard version of OMME which averages the phases in the complex plane.
It can be appreciated that, when β = 1 , OMME with β = 1 / 2 provides the most robust unwrapping with respect to noise.It can be also seen that when further increasing the noise σ ϕ , OMME decreases its robustness.It can also appreciated that, for β = 1 / 2 , 1 , increasing N also increases the robustness when d e f f → u true for each level of noise.
Of course, the unwrapping capabilities of OMME fail when the noise grows.In that case, it is recommendable to just repeat the experiment with the same dynamic range and average the results, as done for instance in Fig. 2 f.This of course provides only a linear decrease in the variance with respect to the number of measurements N, as shown in Fig. 2 , instead of exponential as in OMME.Naturally, both strategies could be eventually combined (if enough scan time available) by averaging first for each MEG and then use OMME with β = 1 / 2 with a better SNR.

On the choice of β and N
If the user sets the desired dynamic range d e f f (e.g.20 μm in MRE), as presented above the best combination of β and N depends on σ ϕ .
If the dynamic range at the desired MNR is d N , and we fix the proportion β = a/b, the number of gradients can be found by using the relation and hence the d 1 image can still contain wraps while d e f f may not.However, what we recommend is first to take β = 1 / 2 , which provides the best unwrapping robustness when increasing SNR: in the next sections we show that in real MRE images β > 1 / 2 often fails already when N = 2 .

Subjects
In vivo MRE was performed in eight healthy men without a history of neurological diseases (mean age ± SD: 36 ± 9 years).The study was approved by the ethics committee of Charite-Universitaetsmedizin Berlin in accordance with the Ethical Principles for Medical Research Involving Human Subjects of the World Medical Association Declaration of Helsinki.Every participant gave written informed consent.

OMME-MRE experimental setup
All experiments were performed in a 3T MRI scanner (Siemens Magnetom Lumina, Erlangen, Germany).In order to separate the different contributions to the background phase as indicated in Eqs. ( 1) and (2) , four scans were consecutively acquired in each subject as summarized in Table 1 .
Measurements with harmonic vibrations sampled eight phase offsets equally spaced over a vibration period using pressurized air drivers as described elsewhere ( Schrank et al., 2020b ).The vibrations were induced with a forerun of 2 s before MRE data acquisition was started in order to establish a steady state of timeharmonic oscillations throughout the brain.Measurements with active MEGs were conducted for three spatial directions along head-feet, left-right and anterior-posterior consecutively and were repeated for each single MEG amplitude.The amplitude of the harmonic vibrations was the same in all individual acquisitions.It was tuned to avoid signal voids due to intra-voxel phase dispersion for the highest MEG amplitude and to not show any phase wraps for the MEG amplitude of 8 mT/m.In three subjects the measurement without vibration was repeated once with same MEG polarity (δ G + ) and once with opposite MEG polarity (δ G − ) to investigate the influence of MEG polarity on the induced static background phase and further detailed in Section 3.5 .
Using the same polarity has an advantage when doing dual multiple encoding.For a given dynamic range d G , either {−G/ 2 , G/ 2 } (symmetric) or { G ,background} (non-symmetric) need to be measured to perform the phase contrast and remove the background phase.In dual or multiple-encoding, this needs to be done before the reconstruction and not afterwards e.g. by the Fourier analysis.
With regard to Table 1 it can be appreciated that if, one full MEG measurement requires 1 scan time unit, the symmetric phase-contrast approach requires 2 while the latter would require 1 + 1 / 3 + 1 / 8 + 1 / 24 = 1 .5 , hence being more time-effective already for one single motion encoding.Of course, the symmetric approach has the advantage of being able to decrease the dynamic range with the same MEG amplitude, hence potentially being able to obtain better MNR (but more wraps).Also, it only performs two subtractions, instead of 4, leading to a better MNR.However, the MNR of the non-symmetric approach can be brought to the value of the symmetric case when φ 0 is smoothed and the imaging gradients are neglected for large MEGs as it is presented later in Section 3.8 .
In multiple encoding, for N dynamic ranges, the symmetric approach would lead to 2 N scan time units.The non-symmetric approach, would lead to N + 1 / 3 + N/ 8 + 1 / 24 = N(9 / 8) + 9 / 24 scan time units.This is a considerable advantage in terms of scalability of the non-symmetric approach with respect to the symmetric one.

OMME-MRE sequence
Single frequency MRE using a single-shot, spin-echo echoplanar imaging (EPI) sequence was performed for harmonic vibrations at 31.25 Hz. 17 axial slices with an interslice gap of twice the slice thickness were recorded using GRAPPA parallel acquisition ( Griswold et al., 2002 ) with acceleration factor of 2. Slice positioning was automatically done using the scanner build-in auto align function based on the localizer scan (head-brain).Further imaging parameters were: field of view 202 × 202 mm 2 , voxel size 1.6 × 1.6 × 1.6 mm 3 , echo time (TE) of 74 msec and repetition time (TR) of 2500 msec.Three components of the wavefield in orthogonal directions were acquired with first order flow-compensated MEGs of varying amplitude (32, 24, 16, 8, 4, 2 mT/m) and a fixed frequency of 34 Hz with slew rate of 125 mT/m/ms.The corresponding dynamic ranges were 7, 9, 13, 26, 52, 104 μm.The dynamic range of the imaging gradients was 149 μm.Each time the MEG amplitude was changed, one preparation scan was performed to reduce transient effects of eddy-currents.
Acquisition time for a set of 3D MRE data was approximately 6:55 min ( 70s per MEG amplitude, vibration on and MEG on).Additional acquisition time for the individual background phase contributions was 3:07 min.

Motion correction and segmentation
Complex MR images were corrected for stochastic head motion in the range of ± 2 mm using SPM12 Penny et al. (2011) .Moreover measurements without vibration were registered to corresponding measurements with vibration, since deflated actuators results in a vertical displacement of axial slices in the order of 1-2 mm with respect to the inflated actuators during vibration.Automatic segmentation of white matter (WM) and gray matter (GM) based on averaged MRE magnitude images was done using SPM12.
The tissue probability maps were thresholded at 0.8 for WM and 0.9 for GM to generate logical tissue-associated voxel masks.The GM threshold was higher to reduce boundary artifacts at cortical GM/fluid boundaries (see Fig. 3 ).

Reconstruction of phase contributions
The individual phase contributions in Eq. ( 2) were recovered by a number of subtractions.ϕ 0 (measurement 2, see Table 1 ) was subtracted from measurement 3 to determine the static MEGdependent phase δ G .Subtraction of measurement 3 ( ϕ 0 ) and measurement 4 gave the motion-dependent phase encoded by the imaging gradients m (u ) .The static background phase components ϕ 0 and δ G were smoothed using a Gaussian filter with 1 mm standard deviation in order to reduce noise enhancement by further subtraction of these components.This was justified since both static phases show low spatial variations within the brain.Finally we subtracted the individual phase contributions from measurement 1 to determine the tissue displacement encoded by the MEG only (1) .From the repeated measurements without vibrations and either same or opposite MEG polarity, we determined pixel-wise the relative absolute error (RAE) (RAE ) .The RAE was then averaged over WM and GM and subjects for each MEG amplitude.

Displacement reconstruction
Single encoding phase contrast images were computed for each MEG using Eq.(3) (assuming no noise).The background phase was obtained as detailed in Section 3.5 .Dual and tri-encoding phase images were computed using the OMME formula (21) .Dual encoding phase images were computed using the combinations of two single encoding images, namely 32 and 24 mT/m, 24 and 16 mT/m, 16 and 8 mT/m, 32 and 8 mT/m.In addition, OMME was used to combine three phase images acquired with MEG amplitudes of 32, 16, 8 mT/m and 32, 24, 16 mT/m.
As shown in the theory section all these combinations exhibit the same dynamic range d G given by the lowest encoding amplitude of 8 mT/m, which had no more phase wraps.Notice that due to the inclusion of m (u ) in the background phase for all MEGs, the phase difference measurements are not i.i.d. as assumed in the noise analysis.However, recall that Eq. ( 21) is used for the reconstruction and therefore the unwrapped image does not result in the combination of phase differences any more.Therefore, the measurements not being i.i.d.does not affect the variance of the reconstruction.
We determined the number of wrongly reconstructed voxel inside WM and GM tissue for each combination of MEG amplitudes in order to assess the noise sensitivity of the different combination possibilities in vivo as it was simulated before (see Fig. 1 ).We defined the single phase-contrast image with a MEG of 8 mT/m as our ground truth and calculated the voxel wise difference to the multiple MEG phase reconstructions.Based on the noise level in the image and the maximum encoded displacement, a threshold of 0.1 rad phase difference was used to identify wrongly reconstructed voxel in WM and GM.Relative error rates were determined by dividing the number of wrongly reconstructed voxels by the total number of voxels included in the GM and WM masks in all slices, timesteps and encoding directions.
To further investigate the noise sensitivity of the displacement reconstruction, we added complex Gaussian noise with a standard deviation of 15% of the mean absolute encoded phase in WM and GM to the single PC images and repeated the evaluations.
Furthermore wrapped single motion encoding phase contrast images for the highest MEG of 32 mT/m were unwrapped using Laplacian and Flynn based unwrapping algorithms.We chose to compare our proposed method with Flynn and Laplacian based unwrapping to include two common but different approaches which are publicity available at https://bioqic-apps.charite.de .We compared the different unwrapping approaches in terms of MNR as described further below and in terms of the visual quality of the reconstructed elastograms as outlined in the next section.

Shear wave speed reconstruction
Wrap-free phase images from unwrapping algorithms and from dual and multiple encoding methods were used for reconstruction of shear wave speed (SWS) maps based on phase-gradient wavenumber recovery to avoid noise amplification by the Laplacian operator which is inevitable in direct inversion techniques ( Hirsch et al., 2017;Mura et al., 2020 ).SWS is related to tissue stiffness and will be termed as such in the following.The principle of wavenumber (k-) based multi-component, elasto-visco (k-MDEV) inversion was originally introduced for liver MRE and is outlined in Tzschätzsch et al. (2016) .
It is important to note that each reconstructed voxel of the elastograms resulted directly from 24 individual voxels of the phase images (8 timesteps and three encoding directions) and indirectly from their surrounding voxels as well.If only one voxel in the phase images is wrongly reconstructed, the resulting elastogram voxel is corrupted.Therefore we analysed additionally the wrongly reconstructed voxels with respect to the elastograms for the comparison of different multiple encoding approaches.To calculate relative error rates, we divided again the number of wrongly reconstructed voxels by the number of all directly affecting voxels.In contrast to the phase images, the number was not multiplied by the amount of timesteps and encoding directions.k-MDEV inversion was adapted to the resolution of brain MRE as outlined in Herthum et al. (2021) .Compared to k-MDEV proposed for abdominal organs ( Tzschätzsch et al., 2016 ), smoothing the phase images prior to the unwrapping was omitted since this would have influenced the MNR estimations.Moreover, the linear radial filter in the spatial frequency domain was replaced by a radial bandpass Butterworth filter of third order with highpass threshold of 15 1/m and lowpass threshold of 250 1/m.
being the motion u encoded by the MEG with dynamic range d G and the background phase due to the static phase of the imaging gradients ϕ 0 , the static MEG-dependent phase δ G and the motion-dependent phase encoded by the imaging gradients m (u ) .Additionally the static-background phase δ −G for toggled MEG polarity is shown.MRE mean magnitude image and masks for white matter (WM) and gray matter (GM) are given as reference.The color scale of the phase images was adapted at each figure for better visualization.

Noise reduction by adding back imaging gradient's phase
In the OMME context, the subtraction of m (u ) is needed for the correct phase contrast when including several gradient strengths.However, it is a common practice to assume that the phase contribution m (u ) is small with respect to the contribution of the wave motion for the largest MEG (i.e.smallest dynamic range d N ).In such cases adding back the phase contribution of the imaging gradient's to the OMME reconstruction allows theoretically for a reduction factor 1 / √ 2 ≈ 0 .7 in the standard deviation of the noise of the displacement field.
Therefore, the displacements obtained with OMME are postprocessed by the following operation This effect is compared quantitatively in terms of MNR as outlined in the next section and qualitatively on the resulting elastograms in Fig. 5 .All other results, elastograms and wave fields are without re-added m(u).

Noise analysis and statistical tests
Signal power and MNR of the phase images are important parameters for the subsequent post-processing and final SWS reconstruction.According to our theory, OMME promises wrap-free phase images with MNR corresponding to the highest MEG used for OMME phase recovery.To calculate MNR for experimental data (unwrapped and unsmoothed phase images) we used the blind noise estimation method from Donoho et al. (1995) as outlined and previously applied to MRE data in Bertalan et al. (2019) .Noise estimation in the wavelet domain is expected to be well suited for wave images ( Barnhill et al., 2017;Selesnick et al., 2005 ).We estimated MNR from the dual-tree wavelet transformation of the displacement images with the median absolute deviation of the finest band of wavelet coefficients ( Donoho et al., 1995 ).The estimated signal power was derived from the L2-norm.Signal and noise levels were estimated from automatically segmented WM and GM regions (see Fig. 3 ) for all slices and components and averaged afterwards.
To test for significant differences in the number of reconstruction failures using OMME and dual encoding strategies, a linear mixed-effects model with varying intercept was employed.Error rates were used as dependent variables and the different methods as independent variables.Participants were assigned as random effect.To test for significant differences in the MNR of unwrapped phase images using OMME, Laplacian and Flynn unwrapping, a linear mixed-effects model with varying intercept was employed.MNR was used as dependent variables and the different methods as independent variables.Participants were assigned as random effect.All P -values were calculated using Tukey's post hoc test with Bonferroni correction for multiple comparisons.All statistical analysis was done in R (version 4.0.2).Unless otherwise stated, errors are given as standard deviation (SD).P -values below 0.05 were considered statistically significant.

Phase images
Fig. 3 shows the encoded phase of the complex MR signal, with the different contributions modeled by Eqs. ( 1) and (2) derived from the measurements listed in Table 1 .One central slice of the anterior-posterior encoding direction is displayed for one representative subject.The third column shows the static background phase induced by toggled MEGs.For reference the MRE mean magnitude and masks for WM and GM are given.Table 2 summarizes the encoding efficiency for the different MEG amplitudes and the imaging gradients.Group mean absolute displacement for all encoding directions averaged over WM tissue in rad is given.Furthermore the averaged RAE of the background phase for the same and opposite MEG polarity is tabled, RAE + and RAE − respectively.The encoded phase u π d G increased with increasing MEG amplitude and phase wraps occurred from 16 mT/m on.The static background phase induced by the MEG δ G decreased with amplitude until no difference compared to the background phase induced by the imaging gradients ϕ 0 was visible.Toggling the MEG resulted in a different background phase which is visible for MEG amplitudes.The quantitative analysis using the RAE showed that the difference between repeated measurements with same MEG polarity (32 mT/m: 2 .5 ± 2 .4% ) is lower than with toggled MEG polarity (32 mT/m: 10 .6 ± 7 .8% ).The error decreased with decreasing MEG amplitude.The displacement encoded by the imaging gradients m (u ) was small compared to the displacement encoded by the larger MEGs, although this estimate depends on the applied vibration frequency and is likely higher for higher frequencies.

Dual and multiple encoding unwrapping
Fig. 4 shows the different phase reconstructions of a single timestep for different MEG combinations with the same dynamic range of the MEG with 8 mT/m amplitude for one representative slice.In addition reconstructed elastograms of SWS are displayed.At the top, results are given for the original data, which is the phase encoded by the MEG only.Results with added noise are shown at the bottom.OMME with three MEGs (32, 16, 8 mT/m and 32, 24, 16 mT/m) was compared to dual encoding strategies using 32 and 24 mT/m, 24 and 16 mT/m, 16 and 8 mT/m, 32 and 8 mT/m.It is well visible that dual encoding with 32 and 24 mT/m performed worst in terms of reconstruction failures, which sub-  sequently corrupted the reconstructed elastogram.Despite no apparent reconstruction failures in the selected slice, also the other approaches showed defects in the final elastograms, which resulted from reconstruction failures on other timesteps or components.Moreover encoding approaches using higher MEGs showed less noise in the reconstructed phase image.Adding noise to the complex data before reconstruction increased the number of reconstruction failures and noise of the combined image in all ap-proaches.Consequently more corrupted voxels were visible in the final elastogram.Table 3 summarizes the findings as relative error rates for the phase images ( 3 a) and for the elastograms ( 3 b) compared to the total amount of voxels (mean ± SD: 43 , 639 ± 3 , 114 ) inside the GM and WM mask for each subject.Incorporating all timesteps and encoding directions resulted in a total 1,047,336 voxel which could possibly fail to be reconstructed properly.In addition MNR of re- constructed phase images is tabled ( 3 c).All numbers are given as group average and standard deviations in brackets.
In general, only little reconstruction failures ( < 1% ) were observed in comparison to all possible voxels in the phase images.Only the dual encoding with 32 and 24 mT/m with 15% added noise showed failures above 4% .However, due to the combinatorial nature of SWS reconstruction which combines up to 24 phase images to one SWS image, the error rates became substantial for the elastograms.The relative differences between the reconstruction approaches were conserved.For the original data, dual encoding with 32 and 24 mT/m performed significantly worse ( 3 .9 ± 3 .4% ) than OMME with 32, 16, 8 mT/m ( 0 .4 ± 0 .4% , p = 0 .001 ).There was no statistical difference between the other encoding strategies ( p > 0 .99 ).Nonetheless the MNR scaled with the highest MEG amplitude used, such that approaches with 32 mT/m had a MNR of 17 ± 3 dB, 24 mT/m gave 15 ± 3 dB and 16 mT/m gave 12 ± 3 dB.Adding noise to the original data inflated error rates in all approaches which became larger than 40% for the noise sensitive dual encoding approach with 32 and 24 mT/m.With increased noise, OMME with 32, 16, 8 mT/m ( 1 .5 ± 0 .9% ) also significantly outperformed the dual encoding with 24 and 16 mT/m ( 6 .9 ± 1 .3% , p = 0 .03 ) and 32 and 8 mT/m ( 16 .8 ± 7 .7% , p < 0 .0 0 01 ) in terms of reduced reconstruction failures.These findings are supported in the theory section.OMME with 32, 24, 16 mT/m shows higher error rates in the elastograms with added noise ( 2 .6 ± 1 .4% ) than the combination 32, 16, 8 mT/m, what is expected due to the poor performance of 32, 24 mT/m.Therefore, OMME with 32, 16, 8 mT/m is the most noise robust unwrapping approach tested here.Needless to say, that MNR was reduced when noise was added and MNR differences between the approaches were conserved.The last row of Table 3 indicates an increased MNR with additional postprocessing as described in Section 3.8 .MNR improvements with respect to the original data were achieved in all cases.Fig. 5 shows the corresponding elastograms in three volunteers for the comparison between the original data and the additional postprocessing.The increased MNR is again evident.

Comparison to other unwrapping methods
Fig. 6 shows representative results for the SWS maps reconstructed from wrap-free phase images.The unwrapping was either performed using OMME utilizing phase images from MEG amplitude of 32, 16 and 8 mT/m or by Laplacian and Flynn unwrapping algorithms applied to the PC image of 32 mT/m MEG amplitude.Anatomical reference images are based on T2 weighted MRE magnitude images.Red arrows indicate areas where OMME based SWS reconstruction visually outperforms the other two approaches.Overall the noise outside the brain was largely reduced using OMME and tissue/air interfaces were sharper.Especially the transition between the skull and the brain tissue was properly reconstructed, while the unwrapping methods smoothed that region which lead to spurious stiffness values and reduced contrast.
In the first subject, it was especially difficult to demarcate the tissue/air boundary in the area of the left superior temporal sulcus using SWS reconstruction based on Laplacian and Flynn unwrapping.Only OMME allowed a good boundary detection.A similar effect was visible at the lingula gyrus where the space between the two hemispheres was only preserved properly with OMME.For subject two, the central part of the right lateral ventrical showed spurious SWS values for Laplacian and Flynn unwrapping probably due to tissue/fluid boundary artifacts which were enhanced by the algorithms.OMME based reconstruction showed higher level of details by fully recovering the boundaries between brain tissue and either ventricles or gyri.In the magnitude image of subject three susceptibility artifacts are present.However, OMME based SWS reconstruction showed a good agreement with the anatomical reference and correctly reconstructs SWS values associated with tissue voxels in the area of the temporal pole.In contrast, Laplacian based SWS maps are heavily corrupted and no reference to the anatomical image is present.Flynn performs better but still with noisy SWS values and blurred CSF/tissue (solid/fluid) boundaries.Similar observations are visible in a more cranial area of the temporal pole in subject four.Noisy SWS values make the demarcation of the temporal pole difficult for Laplacian and Flynn unwrapping SWS reconstructions.
The MNR analysis based on wrap-free phase images with a MEG amplitude of 32 mT/m revealed for Laplacian unwrapping a group mean MNR of 15 .9 ± 2 .7 dB and for Flynn unwrapping 15 .6 ± 2 .4 dB.Both results are significantly lower than the MNR for OMME based unwrapping as listed in Table 3 ( p = 0 .02 ).

Discussion
We have developed, theoretically analyzed and assessed in numerical and human brain data a new method for combining an arbitrary number of motion-encoded PC-MRI images, the Optimal Multiple Motion Encoding method (OMME).
We compared the proposed method with dual motion encoding strategies and common phase unwrapping algorithms in terms of unwrapping success, MNR and quality of subsequently reconstructed SWS maps.To the best of the author's knowledge, this is the first reported method for combining a larger number of motion encoded images obtained from different MEGs.
For a fixed effective dynamic range of the encoded motion, OMME presents a superior performance with respect to noise compared to standard dual encoding unwrapping.This was assessed analytically and confirmed numerically in a "single voxel" experiment.The analysis on the in vivo data with respect to reconstruction failures and MNR confirm these findings.Additionally it was shown that inverting the MEG polarity affects the induced background phase of the MEG which is not critical for our proposed method, but it should be considered in classical PC approaches where a phase-difference image is calculated to remove contaminant phase information.
It was shown that unwrapping is most robust to noise when N images are combined which were measured in the dynamic ranges d 1 , . . ., d N such that d i = 2 −i +1 d 1 .This simplifies the acquisition protocol allowing the scanner operator to select the largest MEG and the number of measurements N only, as it is usually done when the MEG is kept fixed.
The OMME was compared against standard unwrapping methods (Laplacian and Flynn).Remarkably, OMME allows to improve the SWS maps by reducing the noise in the wave images without spatial smoothing as Laplace unwrapping does and without unwrapping failure as it may occur with Flynn unwrapping predominantly at boundaries.This showed that details can be preserved which are otherwise smoothed (out) by standard unwrapping methods.This can be relevant for higher resolution MRE in a variety of applications including tumor detection or characterization of lesion in multiple sclerosis (MS) Streitberger et al. (2012) .Moreover, we showed that standard unwrapping methods smear boundaries between fluid filled spaces and brain tissue.This not only affects cortical areas of the brain and their tissue/air boundaries but also interfaces between tissue and fluid filled ventricles.The importance of proper reconstruction of stiffness estimates for cortical areas has recently been addresses by Lilaj et al. (2021) .If tissue mechanical properties are altered at those boundaries, e.g. as a result of impaired CSF-brain barriers in MS Takeoka et al. (1983) OMME based wrap-free MRE phase im-ages could be sensitive to those alterations.Also other interfaces between tumor and healthy tissue could potentially be better resolved.Further, the increased dynamic range of OMME with good MNR properties could be utilized when high frequency vibrations induce heavy wraps near the source and are quickly damped towards small deflection amplitudes inside the tissue under investigation.The potential of OMME for higher frequency MRE needs to be further investigated.Without heavy wraps, Laplacian and Flynn unwrapping methods performed similar, which underlines that OMME might be very suitable for high MNR applications.
The postprocessing introduced in Section 3.8 allows for an improvement in the SWS image quality at no additional scan time and at negligible computational cost.Since this assumes large MEG encoding efficiencies compared with those of the imaging gradients, the applicability should be investigated with respect to each specific scanning protocol.
As a limitation of OMME, examination times are increased by additional measurements for multiple MEGs.Each applied MEG increases the total scan time by the acquisition time of one measurement.Moreover, the background phase (i.e.MEGs on and vibration off) needs to be measured, which still is required at one timestep only for all encoding directions, adding another 1/(number of timesteps) * acquisition time.Nevertheless, the time investment pays off when phase wraps can be avoided and maps are generated that are more detailed than standard methods.Even resolving wraps only partly supports unwrapping algorithms and permits higher encoding efficiencies than standard MRE towards measurement of damped waves without corrupting high amplitude regions.In addition, noise sources other than those mentioned could also affect OMME performance.For example, induced movements due to scanner table vibrations are encoded as a function of MEG amplitude and cannot be corrected.However, this was not a problem here because the wave amplitudes were probably strong enough.
OMME can be also applied to other PC-MRI methods like e.g.flow MRI.In that case the dynamic range will be the v enc parameter.However, some careful noise analysis may be needed when the then phase that does not depend on the motion is measured only once, as it is the case in 4D Flow, since then the phase differences for each v enc will be correlated.This might be investigated in a future work.

Conclusion
In this study, we proposed an optimal multiple motion encoding (OMME) method which is suitable for motion sensitive PC-MRI.A detailed theoretical analysis was provided to derive a rationale for choosing the combinations of motion encoding gradients that lead to robust unwrapping for a given effective dynamic range when SNR in the phase images decreases.We applied novel OMME to MRE measurements of in vivo human brain acquisitions.It was shown that OMME outperforms dual encoding strategies and allows to recover more tissue details due to its increased MNR ratio within a high dynamic range leading to SWS maps which preserve important details such as discontinuities in the stiffness.Especially for applications of high resolution MRE wrap-free images with proper MNR -as provided by OMME -are desired.
As we showed in the theory section, OMME will fail in unwrapping to the given dynamic range when further increasing the noise in the phase images.Depending on the value of beta, images may be wrap-free or not, what is given by the relation between the effective dynamic range and beta.

Fig. 3 .
Fig. 3. Measured phase ϕ G (u ) of the complex MR signal in one representative slice, encoding direction (anterior-posterior) and subject for different MEG ampltiudes ranging from 32 mT/m to 2 mT/m.The separated phase contributions correspond to the model ϕ G (u ) = ϕ 0 + δ G + m (u ) + u π

Fig. 4 .
Fig. 4. Different phase reconstructions of a single timestep with the same dynamic range d e f f using multiple motion encoding measurements.Here OMME with three MEGs (32, 16, 8 mT/m) and (32, 24, 16 mT/m) is compared to dual encoding using 32 and 24 mT/m, 24 and 16 mT/m, 16 and 8 mT/m, 32 and 8 mT/m.The bottom row shows the same reconstructions with added Gaussian noise with a standard deviation of 15% of the mean absolute encoded phase in WM.

Fig. 5 .
Fig.5.SWS maps based on wrap-free phase images using OMME (32, 16, 8 mT/m), without (mid column) and with (right column) applying Formula (23) for selected slices in three subjects.The anatomical reference image from T2 weighted MRE magnitude is included.

Fig. 6 .
Fig. 6.SWS maps based on wrap-free phase images using OMME (32, 16, 8 mT/m), Laplacian unwrapping and Flynn unwrapping for selected slices in four subjects.The anatomical reference image from T2 weighted MRE magnitude is included.Red arrows indicate areas where OMME shows more details and greater contrast in the SWS map.
and for β = 1 u pc is not defined and d 1 = d e f f .
that u true is a global minimum of J dual in the noise-free case.Now in Appendix A.1 , we prove that in the general case -and possibly with noisy dataunwrapping is produced by the fact d dual = lcm (d 1 , d 2 ) , with lcm dual (u ) = J dual (u + d dual ) when d 1 /d 2 ∈ Q ..Moreover, it can be also noted that d dual = d e f f .Indeed, writing β = a/b, a < b ∈ N we obtain The theoretical standard deviations curves are the ones following the derivations in the previous sections: √ • The curves are drawn by repeating this procedure in the interval d e f f ∈ [1 , 4] .•

Table 1
Measurement strategy for determining different phase contributions and for the dual motion encoding and OMME reconstructions.

Table 2
Encoding efficiency and mean absolute encoded displacement as group averages for different MEG am plitudes and the imaging gradients.In addition, averaged relative absolute differences RAE for repeated measurements with same and opposite polarity are given.Standard deviations are given in brackets.

Table 3
Number of wrongly reconstructed voxels in % in phase images and elastograms and MNR for OMME using three MEGs and different dual encoding strategies.Group mean values were averaged over WM and GM and tabled as group mean (sd).All combinations exhibit the same dynamic range d G with different noise sensitivities to the input image noise σϕ and different noise levels of the reconstructed phase images (MNR).In addition results with added Gaussian noise (15% of the mean absolute encoded phase in WM) and added m(u) are given.MEG amplitudes are given in mT/m.