VENC DESIGN AND VELOCITY ESTIMATION FOR PHASE CONTRAST MRI

Abstract

Systems and methods to implement Phase Recovery from Multiple Wrapped Measurements (PRoM) as a fast, approximate maximum likelihood estimator of velocity from multi-coil data with possible amplitude attenuation due to dephasing. The estimator can recover the fullest possible extent of unambiguous velocities, which can exceed four times the highest venc. The estimator uses all pairwise phase differences and the inherent correlations among them to minimize the estimation error. Derivation of the estimator yields explicit probabilities of unwrapping errors and the posterior probability distribution for the velocity estimate; this in turn allows for optimized design of the phase-encoded acquisition.

Description

FIELD OF THE INVENTION

The present invention relates to magnetic resonance imaging (MRI) and more specifically to a fast, approximate maximum likelihood estimator of velocity from multi-coil data and optimized data acquisition.

BACKGROUND

Phase-contrast magnetic resonance imaging (PC-MRI) is a quantitative, non-invasive technique to measure hemodynamics in vivo. PC-MRI also enables higher dimensional velocimetry such as 4D flow imaging. Spin velocity is encoded into voxel phase via a time-varying gradient field. Lowering the first moment of the encoding gradient extends the unaliased range but degrades the velocity-to-noise ratio (VNR). To address this issue, multi-point acquisitions such as “dual-venc” have been proposed, whereby a high-venc measurement is used to unwrap a potentially wrapped, but less noisy, low-venc velocity measurement. In contrast, multiple (possibly wrapped) pairwise phase differences can be jointly processed, leading to a potentially larger unambiguous range of velocities and a reduced estimation error. Existing estimators usually employ a computationally expensive grid search and cannot guide the optimized design of the acquisition.

SUMMARY

Accordingly, the present disclosure discloses Phase Recovery from Multiple Wrapped Measurements (PRoM), an approximate maximum likelihood estimator (MLE) for a set of linear congruence equations with additive correlated noise together with implementations thereof. The estimator is used in multi-point PC-MRI to jointly process all pairwise phase differences. The PRoM estimator first constructs a set of candidate tuples of wrapping integers. The probability that the true tuple of wrapping integers is in this set is arbitrarily close to 1. For each candidate tuple, the corresponding candidate velocity is found without grid searching as a simple weighted combination of the noisy measurements. The final velocity estimate is chosen among the small set of candidate velocities to maximize the likelihood function. The search for the tuple of wrapping integers can be accelerated as a k-d tree search. The approximate MLE explicitly accommodates both intravoxel dephasing and the inherent noise correlation present among phase differences. Compared to simplified grid-search MLE, the computation of PRoM is orders of magnitude faster.

The probability distribution of the velocity estimate from noisy data is derived in closed form, which allows for the optimized design of the phase encoding. Additionally, the same estimation problem appears in other array processing applications, where PRoM extends current art to provide: a fast, grid-free estimator; accommodation of correlated noise; and, principled design of sparse array geometry. The probability distribution also enables spatial processing to exploit spatial correlations among velocity values, providing further robustness to noise.

The foregoing illustrative summary, as well as other exemplary objectives and/or advantages, and the manner in which the same are accomplished, are further explained within the following detailed description and its accompanying drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

A detailed description of certain aspects of the present disclosure in accordance with various example implementations will now be provided with reference to the accompanying drawings. The drawings form a part hereof and show, by way of illustration, specific implementations and examples. In referring to the drawings, like numerals represent like elements throughout the several figures.

FIG. 1 illustrates the sum of squared error versus velocity using a single-coil, symmetric three-point acquisition;

FIG. 2 illustrates similarity metric results, which shows agreement for s₂₁>2;

FIG. 3 illustrates ({tilde over (v)},{circumflex over (v)}) for venc=[35,10,14]^T cm/s;

FIG. 4 provides a visualization of wrapping integers k*({tilde over (v)}) for the case venc=[35,10,14]^T, [s₁₁,s₂₁,s₃₁]=[2.5,5,2.5];

FIGS. 5A and 5B illustrate a histogram of {circumflex over (v)}|v=0 using 10⁵ trials, N_e=3, N_c=1, venc=[99,18,22]^T cm/s at s₂=5,10;

FIG. 6 illustrates an example method of implementing PRoM for N_e-point Encoding;

FIG. 7 illustrates that PRoM very closely approximates the maximum likelihood estimator;

FIGS. 8A-8B show the RMSE results averaged over 10⁵ trials at each true velocity value with grid 0.1 cm/s;

FIGS. 9A-9B show the RMSE improvement enabled by optimized venc design;

FIGS. 10A-10H show the result for simulation of intravoxel dephasing;

FIGS. 11A-11C show measured data from the spinning phantom and illustrate the square root of sum of squared coil images for m₁₁, m₁₂, m₁₃;

FIGS. 11D-11G show measured data from the spinning phantom and illustrate velocity estimates from SDV, ODV, PRoM, and PRoM+;

FIGS. 11H-11K are zoomed-in versions of FIGS. FIGS. 11D-11G;

FIGS. 12A-12C show measured data from a healthy volunteer to validate that PRoM and PRoM+ can unwrap a larger range of velocities than existing techniques, given the same encodings and illustrate the square root of sum of squared coil images for m₁₁, m₁₂, m₁₃;

FIGS. 12D-12G show measured data from the healthy volunteer and illustrate velocity estimates from SDV, ODV, PRoM, and PRoM+;

FIGS. 12H-12K are zoomed-in versions of FIGS. FIGS. 12D-12G; and

FIG. 13 is a view illustrating a structure of an example magnetic resonance imaging (MRI) apparatus that may be used to acquire image data.

DETAILED DESCRIPTION

1. Introduction

The present disclosure is directed to a Phase Recovery from Multiple Wrapped Measurements (PRoM), which is a fast, approximate maximum likelihood estimator of velocity from multi-coil data with possible amplitude attenuation due to dephasing. The estimator can recover the fullest possible extent of unambiguous velocities, which can exceed four times the highest venc. The estimator uses all pairwise phase differences and the inherent correlations among them to minimize the estimation error. Derivation of the estimator yields explicit probabilities of unwrapping errors and the posterior probability distribution for the velocity estimate; this in turn allows for optimized design of the phase-encoded acquisition. Simulation, phantom, and in vivo results for three-point PC-MRI acquisitions validate the benefits of fast computation, reduced estimation error, increased recovered velocity range, and optimized acquisition. The examples presented demonstrate two orders of magnitude reduction in computational complexity and 10% decrease in root mean squared error versus conventional “dual-venc” techniques.

II. Theory

A. Phase Encoding

A time-varying magnetic gradient field may be used to encode spin velocity into image phase. Here, the velocity component is encoded in one direction. Consider a spin moving through a magnetic field; the Taylor series expansion of spin position p(t)ϵ at t=0 yields

$\begin{array}{cc}p\left(t\right)=p\left(0\right)+\frac{v}{1!}t+\frac{a}{2!}{t}^2+\dots ,& \left(1\right)\end{array}$

where v is instantaneous velocity and α is acceleration. Let γ be the gyromagnetic ratio, B₀ be the main static magnetic field strength, and g(t)ϵ be the time-varying magnetic field gradient. Then to first order approximation of p(t), the phase accumulated from t=0 to echo time TE is

$\begin{array}{cc}\begin{array}{c}\varphi ={\int }_0^{\mathrm{TE}}\gamma B_a+\gamma g\left(t\right)p\left(t\right)\mathrm{dt}\\ \approx \stackrel{{\varphi }_a}{\overbrace{\gamma B_0\mathrm{TE}+\gamma p\left(0\right)\underset{m_0}{\underbrace{{\int }_0^{\mathrm{TE}}g\left(t\right)\mathrm{dt}}}}}+\mathrm{\gamma \upsilon }\underset{m_1}{\underbrace{{\int }_0^{\mathrm{TE}}g\left(t\right)\mathrm{tdt}}}\\ ={\varphi }_0+\gamma m_1\upsilon ,\end{array}& \left(2\right)\end{array}$

where m₀ and m₁ denote the zeroth and first moments of (t). The zeroth moment, m₀, encodes the spin position into phase and combines with γB₀TE to make the background phase, ϕ₀. The first moment, m₁, encodes spin velocity into phase. So, for N_e-point encoding and N_c coils, the integral of all spins in a voxel yields the measurement

$\begin{array}{cc}{\tilde{y}}_{\alpha \beta }=A_{\alpha }{S}_{\beta }{e}^{i\left({\varphi }_0+\gamma m_{1\alpha }v\right)}+{\delta }_{\mathrm{\alpha \beta }},& \left(3\right)\end{array}$

where αϵ{1, . . . , N_e} indexes over all encodings, and βϵ{1, . . . , N_c} indexes over all coils. The resulting signal amplitudes are A_αϵ; S_βϵ is the coil sensitivity weight; v is the resultant instantaneous velocity; m_1αϵ is the first moment; δ_αβϵ is independent and identically distributed (i.i.d.) complex circularly symmetric Gaussian noise. This i.i.d. assumption can be aided by pre-whitening along the coil dimension. The existence of heterogeneous spin velocities and proton density can make v in (3) differ from the mean moving spin velocity (e.g., partial volume effect) and can reduce the amplitude (e.g., intravoxel dephasing effect). Generally, A_α decreases as |m_1α| increases.

Let {tilde over (Y)} denote an N_e×N_c complex-valued measurement matrix with (α,β)^th entry {tilde over (y)}_αβ, and observe that the unambiguous range of velocities for the N_e encodings is

$\begin{array}{cc}{\Omega }^{\prime }=\mathrm{LCM}\left(\frac{2\pi }{\gamma m_{11}},\dots ,\frac{2\pi }{\gamma m_{1N_e}}\right),& \left(4\right)\end{array}$

where LCM(·) denotes least common multiple, which is the smallest positive real number that is an integer multiple of all input numbers. It follows that v and v+Ω′ are indistinguishable given data, {tilde over (Y)}. Considering noise, the MLE of v involves N_e+N_c+2 unknowns, {v, ϕ₀, A₁, . . . , A_Ne, S₁, . . . , S_Nc}, and is a nonlinear least-squares fit to the data. The optimization task can be reduced to

$\begin{array}{cc}\underset{s}{\arg \max }\frac{s^H\tilde{R}s}{s^{H}s}s.t.s=1,& \left(5\right)\end{array}$

where Ŕϵ^Ne×Ne and the “steering vector” sϵ^Ne×1 are

$\begin{array}{cc}\begin{array}{c}\tilde{R}=\tilde{Y}{\tilde{Y}}^H\\ s={\left[A_1{e}^{i\gamma m_{11}v},\dots ,A_{N_e}{e}^{i\gamma m_{1N_e}v}\right]}^{\top }.\end{array}& \left(6\right)\end{array}$

Derivation of (5) is a straightforward extension of known results to accommodate unequal amplitudes, A_α. Here, (·)^T and (·)^H denotes transpose, and conjugate transpose. The steering vector s is determined up to scale. For a given v, the amplitudes may be found by solving (5) as a generalized eigenvalue problem; yet the optimization over v nonetheless encounters a difficult cost surface with many local minima.

A simple illustration of the sum of squared error versus velocity using a single-coil, symmetric three-point acquisition is shown in FIG. 1, which shows a sum of squared error versus velocity for a single coil, symmetric three-point acquisition.

$\gamma \left[m_{11},m_{12},m_{13}\right]=\left[-\frac{\pi }{20},\frac{\pi }{70},\frac{\pi }{20}\right]s/\mathrm{cm},\left[\frac{❘A_1{S}_1❘}{\sigma \left({\delta }_{11}\right)},\frac{❘A_2{S}_1❘}{\sigma \left({\delta }_{21}\right)},\frac{❘A_3{S}_1❘}{\sigma \left({\delta }_{31}\right)}\right]=\left[2.5,5,2.5\right]$

The fit error to the noisy data at the true velocity v=0 cm/s is shown by the red diamond. The global minimum at v≈−1.01 cm/s for the given noise realization is shown by the green marker. Competing local minima, shown as light blue markers, may become the global minimum due to noise, thereby causing unwrapping errors.

Data are simulated using

$\left[\frac{❘A_1{S}_1❘}{\sigma \left({\delta }_{11}\right)},\frac{❘A_2{S}_1❘}{\sigma \left({\delta }_{21}\right)},\frac{❘A_3{S}_1❘}{\sigma \left({\delta }_{31}\right)}\right]=\left[2.5,5,2.5\right]$

to represent a moderate 50% loss of signal amplitude due to intravoxel dephasing. Throughout the disclosure, the same 50% is used in simulation. The fit error at the true velocity, v=0 cm/s, is shown by the red diamond, and the global minimum at v≈−1.01 cm/s for the given noise realization is shown by the green marker. Competing local minima, shown as blue markers, may become the global minimum due to noise, thereby causing unwrapping errors.

B. Estimation of Phase Noise Covariance

From (5), it can be understood that the {tilde over (R)} is a sufficient statistic for v in (3). Departing from the multi-variate optimization in (5) to instead work with the phases of the off-diagonal entries of {tilde over (R)}, the amplitudes of {tilde over (Y)} may be used to inform the phase-to-noise ratio. In so doing, there are advantages: (1) fast, grid-free parameter estimator for velocity, v; (2) characterization of the unambiguous set of estimated velocities; (3) characterization of the probability of unwrapping errors; (4) ability to design the encodings [m₁₁, . . . , m_1Ne] to minimize mean squared estimation error subject to guarantees on the probability of unwrapping error and required unambiguous range. Further, the proposed surrogate estimator is asymptotically efficient, providing a very good approximation to the MLE. In this section, an approximate covariance matrix is derived for the phase measurements in {tilde over (R)} to lay the groundwork for building the proposed estimator of v.

Observe that {tilde over (R)}={tilde over (Y)}{tilde over (Y)}^H performs coil combining and phase differencing. Denote the (a, b)^th entry of {tilde over (R)} as {tilde over (r)}_ab, so

$\begin{array}{cc}{\tilde{r}}_{\mathrm{ab}}={\sum }_{\beta }{\tilde{y}}_{a\beta }{\tilde{y}}_{b\beta }^{*},& \left(7\right)\end{array}$

where (·)* denotes complex conjugation. The noisy phase difference {tilde over (θ)}_ab for encodings a>bϵ{1, . . . , N_e} is

$\begin{array}{cc}{\tilde{\theta }}_{\mathrm{ab}}=\angle {\tilde{r}}_{\mathrm{ab}},& \left(8\right)\end{array}$

where ∠(·) denotes angle of a complex number. This phase differencing results in venc value

$\begin{array}{cc}{\mathrm{venc}}_{\mathrm{ab}}=\frac{\pi }{\gamma \left(m_{1a}-m_{1b}\right)}.& \left(9\right)\end{array}$

There are

$\frac{N_e\left(N_e-1\right)}{2}$

such combinations. Throughout this disclosure, venc_ab has units cm/s. Multiplying the vencs with phases, a (possibly wrapped) noisy velocity {tilde over (v)}_ab may be obtained:

$\begin{array}{cc}{\tilde{v}}_{\mathrm{ab}}=\frac{{\tilde{\theta }}_{\mathrm{ab}}}{\pi }{\mathrm{venc}}_{\mathrm{ab}}.& \left(10\right)\end{array}$

Phases {tilde over (θ)}_ab are unambiguous on any interval of length 2π, and for convenience [0,2π) may be used; so, {tilde over (v)}_abϵ[0,2venc_ab). Thus, a set of noisy congruence equations for v are as follows:

$\begin{array}{cc}{\tilde{v}}_{\mathrm{ab}}\equiv v+n_{\mathrm{ab}}\mathrm{mod}2{\mathrm{venc}}_{\mathrm{ab}},& \left(11\right)\end{array}$

where n_ab is additive zero mean noise. Define k_ab as the wrapping integer for {tilde over (v)}_ab. Momentarily assume that the true wrapping integer, k_ab, is known; then, (11) may be rewritten as a set of linear equations:

$\begin{array}{cc}\tilde{v}+k·2\mathrm{venc}=v1+n,& \left(12\right)\end{array}$

where ∘ is the Hadamard (elementwise) product and

$\begin{array}{cc}\tilde{v}={\left[{\tilde{v}}_{21},\dots ,{\tilde{v}}_{N_e\left(N_e-1\right)}\right]}^{\tau }& \left(13\right)\end{array}$ $k={\left[k_{21},\dots ,k_{N_e\left(N_e-1\right)}\right]}^{\tau }$ $\mathrm{venc}={\left[{\mathrm{venc}}_{21},\dots ,{\mathrm{venc}}_{N_e\left(N_e-1\right)}\right]}^{\tau }$ $1={\left[1,\dots ,1\right]}^{\tau }$ $n={\left[n_{21},\dots ,n_{N_e\left(N_e-1\right)}\right]}^{\tau }.$

From the Chinese remainder theorem, the smallest Ω>0 such that v+Ω satisfies (11) is

$\begin{array}{cc}\Omega =\mathrm{LCM}\left(2\mathrm{venc}\right).& \left(14\right)\end{array}$

This also implies that Ω is the smallest repetition period of v to construct the same {tilde over (R)}. Recall from (4) that Ω′ is a repetition period of v to construct the same {tilde over (Y)}, as well as the same {tilde over (R)}. So, Ω′ is an integer multiple of Ω.

Similar to the assumption {tilde over (θ)}_abϵ[0,2π), assume vϵ[0,Ω) for convenience of derivation. This interval could also be shifted to

$[-\frac{\Omega }{2},\frac{\Omega }{2}),$

for example, for bidirectional flow in MRI.

Let {circumflex over (v)}ϵ[0, Ω) be the linear unbiased estimator of v with smallest root mean squared error (RMSE). To compute {circumflex over (v)} from the noisy {tilde over (v)} in (12), only the covariance matrix of the noise in the remainders, Σ(n) is need. Define a to be the standard deviation of the i.i.d. noise δ_αβ in (3). The mean and variance of {tilde over (r)}_ab are:

$\begin{array}{cc}𝔼\left({\tilde{r}}_{\mathrm{ab}}\right)={\sum }_{\beta }{A}_a{A}_b{e}^{i\gamma \left(m_{1a}-m_{1b}\right)v}{❘S_{\beta }❘}^2& \left(15\right)\end{array}$ $\begin{array}{cc}{\sigma }^2\left({\tilde{r}}_{\mathrm{ab}}\right)={\sum }_{\beta }\left({\sigma }^2\left(A_a^2+A_b^2\right){❘S_{\beta }❘}^2+{\sigma }^4\right).& \left(16\right)\end{array}$

Then from (15)-(16), the squared signal-to-noise ratio (SNR), or Rician factor, for {tilde over (r)}_ab is

$\begin{array}{cc}\begin{array}{c}{\mathrm{SNR}}^2\left({\tilde{r}}_{\mathrm{ab}}\right)=\frac{{❘𝔼\left({\tilde{r}}_{\mathrm{ab}}\right)❘}^2}{{\sigma }^2\left({\tilde{r}}_{\mathrm{ab}}\right)}\\ =\frac{{\left({\sum }_{\beta }{A}_a{A}_b{❘S_{\beta }❘}^2\right)}^2}{{\sum }_{\beta }\left({\sigma }^2\left(A_a^2+A_b^2\right){❘S_{\beta }❘}^2+{\sigma }^4\right)}\\ =\frac{{\left({\sum }_{\beta }{s}_{a\beta }{s}_{b\beta }\right)}^2}{{\sum }_{\beta }\left(s_{a\beta }^2+s_{b\beta }^2+1\right)},\end{array}& \left(17\right)\end{array}$

where

$s_{\alpha \beta }=\frac{❘A_{\alpha }{S}_{\beta }❘}{\sigma }=\mathrm{SNR}\left({\tilde{y}}_{\alpha \beta }\right).$

As SNR({tilde over (r)}_ab)→∞, {tilde over (θ)}_ab converges in distribution to a Gaussian random variable. On the contrary, as SNR({tilde over (r)}_ab)→0, {tilde over (θ)}_ab converges in distribution to a uniformly distributed random variable on [0,2π). Because

${\lim }_{x\to 0}\frac{x}{\sin \left(x\right)}=1,$

and σ({tilde over (θ)}_ab) conditioned on no wrapping is inversely proportional to SNR({tilde over (r)}_ab),

$\begin{array}{cc}{\sigma }^2\left({\tilde{\theta }}_{\mathrm{ab}}\right)\approx \frac{1}{2{\mathrm{SNR}}^2\widetilde{(}{r}_{\mathrm{ab}})}=\frac{{\sum }_{\beta }\left(s_{a\beta }^2+s_{b\beta }^2+1\right)}{2{\left({\sum }_{\beta }{s}_{a\beta }{s}_{b\beta }\right)}^2}.& \left(18\right)\end{array}$

Similarly, covariance can be obtained given no wrapping

$\begin{array}{cc}\mathrm{cov}\left({\tilde{\theta }}_{\mathrm{ab}},{\tilde{\theta }}_{\mathrm{cb}}\right)\approx \frac{❘\mathrm{cov}\left({\tilde{r}}_{\mathrm{ab}},{\tilde{r}}_{\mathrm{cb}}\right)❘}{{\sum }_{\beta }2A_a{A}_b^2{A}_c{❘S_{\beta }❘}^4}=\frac{N_c}{{\sum }_{\beta }2s_{b\beta }^2}& \left(19\right)\end{array}$ $\mathrm{cov}\left({\tilde{\theta }}_{\mathrm{ab}},{\tilde{\theta }}_{\mathrm{bd}}\right)\approx \frac{-N_c}{{\sum }_{\beta }2s_{b\beta }^2}$

and cov({tilde over (θ)}_ab, {tilde over (θ)}_cd)=0 for two phase differences that do not share a common encoding.

In practice, it may be difficult to accurately estimate the noise power for each voxel, and thus hard to estimate s_αβ. To ameliorate estimation difficulty and use complex measurements {tilde over (y)}_αβ only, the variance and covariance may be approximated and scaled:

$\begin{array}{cc}{\sigma }^2\left(\text{?}\right)\approx \frac{{\sum }_{\beta }\left(\text{?}+\text{?}\right)}{2{\left({\sum }_{\beta }\text{?}\right)}^2}\stackrel{\propto }{\sim }\frac{{\sum }_{\beta }\left(\text{?}+\text{?}\right)}{2{\left({\sum }_{\beta }\text{?}\right)}^2}& \left(20\right)\end{array}$ $\begin{array}{cc}\mathrm{cov}\left({\tilde{\theta }}_{\mathrm{ab}},{\tilde{\theta }}_{\mathrm{bd}}\right)\stackrel{\propto }{\sim }\frac{N_c}{\text{?}2\text{?}},& \left(21\right)\end{array}$ $\text{?}\text{indicates text missing or illegible when filed}$

where denotes “approximately proportional to.” Here, the scaling factors are the same for (20) and (21). Let

$\tilde{\theta }={\left[{\tilde{\theta }}_{21},\dots ,{\tilde{\theta }}_{N_e\left(N_e-1\right)}\right]}^{\tau }.$

Then, with the elementwise approximations above, the scaled approximated Σ({tilde over (θ)}) can be formulated directly from observed pixel magnitudes. This scaling does not affect the estimator {circumflex over (v)}, as seen from (30) below. Finally, because the true covariance matrix is close to rank deficient, the element-wise approximations in (20) and (21) can potentially violate the positive semi-definite property of a covariance matrix. Accordingly, the approximation step can be followed by projection to the closest positive semi-definite matrix. This projection operator, Π(·), for a symmetric matrix M with eigen-decomposition M=VSV⁻¹ is given by

$\begin{array}{cc}\Pi \left(M\right)=V\max \left(S,0\right)V^{-1},& \left(22\right)\end{array}$

where max(S, 0) is applied elementwise to the eigenvalues.

To illustrate accuracy of covariance modeling, the cosine similarity metric is adopted, which is scale invariant. For modeled covariance Π(Σ({acute over (θ)})) and sample covariance {tilde over (Σ)}(θ), the cosine similarity is

$\begin{array}{cc}\frac{\mathrm{Trace}\left({\Pi \left(\sum \left(\tilde{\theta }\right)\right)}^H\sum ^{~}\left(\theta \right)\right)}{{\Pi \left(\sum \left(\tilde{\theta }\right)\right)}_F{\sum ^{~}\left(\theta \right)}_F}.& \left(23\right)\end{array}$

Results are computed for the case of a single-coil, symmetric encoding, and intravoxel dephasing 2s₁₁=s₂₁=2s₃₁; 10⁶ random draws are used at each s₂₁ to provide a low-variance sample covariance matrix and to compute the mean cosine similarity. FIG. 2 shows the similarity metric results, which shows excellent agreement for s₂₁>2. Thus, the covariance model is accurate at modest SNR. Also shown in FIG. 2 is the similarity metric for the (scaled) identity covariance model, which is implicitly employed when using least-squares estimation with phase differences.

From (10), the wrapped velocity measurements are linearly related to the phase differences; thus, the approximated scaled positive semi-definite covariance matrix for the noise, n, in (12) conditioned on no wrapping is

$\begin{array}{cc}\sum \left(n\right)\stackrel{\propto }{~}\frac{1}{{\pi }^2}\mathrm{diag}\left(\mathrm{venc}\right)\Pi \left(\sum \left(\theta \right)\right)\mathrm{diag}\left(\mathrm{venc}\right).& \left(24\right)\end{array}$

C. Best Linear Unbiased Estimator

Based on the covariance derivation in II-B, the estimator can now be formulated. Two suitable approximations can be adopted when the SNR is not extremely low. The first approximation is that n follows a joint Gaussian distribution with covariance matrix given in (23),

$\begin{array}{cc}n~𝒩\left(0,\sum \left(n\right)\right).& \left(25\right)\end{array}$

Denote the velocity estimate {circumflex over (v)} given wrapping integers k as {circumflex over (v)}_k. Then, due to (25), {circumflex over (v)}_k and its resulting RMSE given no wrapping, σ({circumflex over (v)}_k), are

$\begin{array}{cc}{\tilde{v}}_k={\langle w^{\tau }\left(\tilde{v}+k·2\mathrm{venc}\right)\rangle }_{\Omega }& \left(26\right)\end{array}$ $\begin{array}{cc}\sigma \left(v_k\right)=w^{\tau }\sum \left(n\right)w,& \left(27\right)\end{array}$ $\mathrm{where}$ $\begin{array}{cc}{w}^{\tau }=\frac{1^{\tau }{\sum }^{-1}\left(n\right)}{1^{\tau }{\sum }^{-1}\left(n\right)1},& \left(28\right)\end{array}$

and a_b denotes remainders after elementwise modulo α by b. Thus, the estimate {circumflex over (v)}_k is a weighted sum of unwrapped noisy velocities.

The second assumption is

$\begin{array}{cc}\mathbb{P}\left(-\mathrm{venc}\preccurlyeq n\preccurlyeq \mathrm{venc}\right)\approx 1,& \left(29\right)\end{array}$

where is elementwise less than or equal. This approximation yields the likelihood ({tilde over (v)}|v)

$\begin{array}{cc}\mathbb{P}\left(\tilde{v}|v\right)\approx \frac{e^{-\frac{1}{2}{d}_{2\mathrm{venc}}^{\tau }\left(\tilde{v},v1\right){\sum }^{-1}\left(n\right)d_{2\mathrm{venc}}\left(\tilde{v},v1\right)}}{\det \left(2\pi {\sum }^{-1}\left(n\right)\right)}& \left(30\right)\end{array}$

where d_z(x, y) is a “wrapped displacement” between x and y with respect to z,

$\begin{array}{cc}{d}_z\left(x,y\right)x-y-\mathrm{round}\left(\left(x-y\right)z\right)·z,& \left(31\right)\end{array}$

with denoting elementwise division. A search over all possible k may be performed to minimize the negative log likelihood,

$\begin{array}{cc}ℒ\left(\tilde{v},{\hat{v}}_k\right)=\frac{1}{2}{d}_{2\mathrm{venc}}^{\tau }\left(\tilde{v},{\hat{v}}_{k}1\right){\sum }^{-1}\left(n\right)d_{2\mathrm{venc}}\left(\tilde{v},{\hat{v}}_{k}1\right).& \left(32\right)\end{array}$

In II-D below, a fast method to detect the best wrapping integers, k* is presented. In II-F and II-H, below, three-point encoding, N_e=3, is considered in order to provide concrete results and to optimize the design of venc and underlying first moments.

D. PRoM

Below is introduced a fast estimator based on (26) and detection of the wrapping integers, k. The detector and (26) are referred to together as Phase Recovery from Multiple Wrapped Measurements (PRoM). PRoM extends the prior signal processing result to accommodate correlated phase errors and to provide a fast computation of the wrapping integers. Moreover, PRoM provides a grid-free alternative to grid search over v.

Assuming vϵ[0,Ω], define the set ′({tilde over (v)}) of wrapping integers

$\begin{array}{cc}{𝒦}^{\prime }\left(\tilde{v}\right)=\left\{k❘-1\preccurlyeq k\preccurlyeq h\right\}& \left(33\right)\end{array}$ $\begin{array}{cc}h=\Omega 12\mathrm{venc}.& \left(34\right)\end{array}$

By (29), (kϵ′({tilde over (v)}))≈1. So, from (26) and (30), the negative log likelihood can be minimized

$\begin{array}{cc}\hat{v}=\underset{v\in [0,\Omega )}{\arg \min }ℒ\left(\tilde{v},v\right)=\underset{v\in \left\{{\hat{v}}_k❘k\in {𝒦}^{\prime }\left(\tilde{v}\right)\right\}}{\arg \min }ℒ\left(\tilde{v},v\right),& \left(35\right)\end{array}$

with probability ≈1. Next, the equality constraint in (12) is leveraged and combined with the second approximation in (29) to decrease the cardinality of ′({tilde over (v)}), denoted as |′({tilde over (v)})|. We have

$\mathbb{P}\left(-\frac{1}{2}\preccurlyeq k+\left(\tilde{v}-v\right)\left(2\mathrm{venc}\right)\preccurlyeq \frac{1}{2}\right)\approx 1,$

which is equivalent to

$\begin{array}{cc}\mathbb{P}\left(k=\lceil -\frac{1}{2}-\left(\tilde{v}-v\right)\left(2\mathrm{venc}\right)\rceil \right)\approx 1& \left(36\right)\end{array}$

where ┌·┐ is element-wise ceiling function. Then, the pruned search set for k can be expressed

$\begin{array}{cc}𝒦\left(\tilde{v}\right)\{\lceil -\frac{1}{2}-\left(\tilde{v}-v\right)\left(2\mathrm{venc}\right)\rceil |v\in [0,\Omega )\}.& \left(37\right)\end{array}$

So, the parsimonious construction considers only vϵ[0,Ω] such that

$-\frac{1}{2}+\frac{v-{\tilde{v}}_{\mathrm{ab}}}{2{\mathrm{venc}}_{\mathrm{ab}}}$

is integer for any 2venc_ab. The cardinality of the search set

$❘\left(\tilde{v}\right)❘\ge 1^{T}h.$

Thus, the number of searches is bounded by the summation of h instead of product of h. Then, the minimization in (35) can be reduced to

$\begin{array}{cc}\hat{v}=\underset{n\in \left\{{\stackrel{⌢}{\upsilon }}_k|k\in \left(\tilde{v}\right)\right\}}{\arg \min }ℒ\left(\tilde{v},\upsilon \right),& \left(38\right)\end{array}$

with probability ≈1. Together, the construction of the pruned set, ({tilde over (v)}), of candidate wrapping integers and the efficient search over {{tilde over (v)}_k|kϵ({tilde over (v)})} to minimize ({tilde over (v)}, v) comprise PRoM, an approximate MLE of v. For a general N_e-point encoding, PRoM pseudo-code is below:

Algorithm 1 PRoM for N∈-point Encoding
Require: Measurements, {tilde over (Y)}. First moments, m , ..., m .
1: Calculate  , {tilde over (ν)}, Ω,   ({tilde over (ν)}) via (9, 11, 14, 38).
2: Calculate scaled Σ (n) and w via (24, 28).
3: Calculate {circumflex over (ν)} via (26, 39).
Ensure: {circumflex over (ν)}
indicates data missing or illegible when filed

FIG. 3 illustrates ({tilde over (v)}, {circumflex over (v)}) for venc=[35,10,14]^T cm/s. The searched candidates {{circumflex over (v)}_k|kϵ({tilde over (v)})} are marked by superimposed red dots. For this case, |′({tilde over (v)})|=252 and |′({tilde over (v)})|=14; thus, 94% of the search space ′ is bypassed via the proposed construction of ({tilde over (v)}). If n is known to concentrate in a smaller volume compared to the second assumption (29), or v is known and restricted to a range less than Ω, then ({tilde over (v)}) may be further pruned accordingly.

To illustrate the reduction of computation complexity in PRoM, consider the case in FIG. 3. Grid search MLE over velocity with standard spacing of Venc₂/1000 entails computation for 7000 candidate velocities. In contrast, PRoM only requires search over only |({tilde over (v)})|=14 candidates, yielding a 500-times reduction.

PRoM admits a simple geometric interpretation. Observe that the noisy velocity measurement {tilde over (v)} resides in a hyper-rectangle {{tilde over (v)}|0{tilde over (v)}2venc}. The vector of noiseless velocity measurements v1_2venc for vϵ[0, Ω) is a point in the hyper-rectangle lying on wrapped line segments parallel to 1. Then, {circumflex over (v)} is found by an oblique projection of the noisy {tilde over (v)}=v1+n_2venc to the closest line segment. Here, the “oblique projection” is determined by the Σ⁻¹(n) weighted distance. The search for the closest line segment is reduced to search for kϵ({tilde over (v)}).

E. Conditional Distribution of the Estimate

Below, ({circumflex over (v)}|v) is derived, which is the distribution of the estimated velocity given the PP-92,ART true velocity; this somewhat technical derivation in turn enables optimized design of venc and the underlying first moments. To derive the distribution, two lemmas are established. The first lemma indicates that adding the same constant, η, to all noise realization components does not affect the error in detecting the wrapping integers and simply adds η to the PRoM velocity estimate, modulo Ω.

Lemma 1 For ηϵ, let n′=n+η1 and {tilde over (v)}′=1+n′_2vec. Then

$\begin{array}{cc}\hat{v}\left({\tilde{v}}^{\prime }\right)={\langle \hat{v}\left(\tilde{v}\right)+\eta \rangle }_{\Omega }& \left(39\right)\end{array}$ $\begin{array}{cc}{\langle k^{*}\left(\tilde{v}\right)-k\left(\tilde{v}\right)\rangle }_h={\langle k^{*}\left({\tilde{v}}^{\prime }\right)-k\left({\tilde{v}}^{\prime }\right)\rangle }_h.& \left(40\right)\end{array}$ $\mathrm{Proof}.\mathrm{By}\left(31\right)$ $\begin{array}{c}{d}_{2\mathrm{venc}}\left({\tilde{v}}^{\prime },\left(v+\eta \right)1\right)=d_{2\mathrm{venc}}{\langle (\tilde{v}+\mathrm{\eta 1}\rangle }_{2\mathrm{venc}},\left(v+\eta \right)1)\\ =d_{2\mathrm{venc}}\left(\tilde{v}+\mathrm{\eta 1},\left(v+\eta \right)1\right)\\ =d_{2\mathrm{venc}}\left(\tilde{v},v1\right)\end{array}.$

Thus ({tilde over (v)}, v)=({tilde over (v)}′, v+η) and {tilde over (v)}({tilde over (v)}′)=η_Ω. Below estimates using k* versus the true k are estimated along the 1 direction. By (26)

$\begin{array}{cc}{\langle \hat{v}\left(\tilde{v}\right)-w^T\left(v1+n\right)\rangle }_{\Omega }={\langle w^T\left({\langle k^{*}\left(\tilde{v}\right)\rangle }_h·2\mathrm{venc}\right)\rangle }_{\Omega },& \left(41\right)\end{array}$ $\mathrm{and}$ $\begin{array}{cc}{\langle \hat{v}\left({\tilde{v}}^{\prime }\right)-w^T\left(v1+n^{\prime }\right)\rangle }_{\Omega }={\langle w^T\left({\langle k^{*}\left({\tilde{v}}^{\prime }\right)-k\left({\tilde{v}}^{\prime }\right)\rangle }_h·2\mathrm{venc}\right)\rangle }_{\Omega }.& \left(42\right)\end{array}$

By the previously derived {tilde over (v)}({tilde over (v)}′)={circumflex over (v)}({tilde over (v)})+η_Ω, results in

${\langle \hat{v}\left({\tilde{v}}^{\prime }\right)-w^T\left(v1+n^{\prime }\right)\rangle }_{\Omega }={\langle \hat{v}\left(\tilde{v}\right)+\eta -w^T\left(v1+n\right)-\eta \rangle }_{\Omega }={\langle \hat{v}\left(\tilde{v}\right)-w^T\left(v1+n\right)\rangle }_{\Omega }.$

Thus (41) and (42) are equal for all w, and

${\langle k^{\star }\left(\tilde{v}\right)-k\left(\tilde{v}\right)\rangle }_h={\langle k^{\star }\left({\tilde{v}}^{\prime }\right)-k\left({\tilde{v}}^{\prime }\right)\rangle }_h.$

FIG. 4 provides a visualization of wrapping integers k*({tilde over (v)}) for the case venc=[35,10,14]¹, [s₁₁, s₂₁, s₃₁]=[2.5,5,2.5]. All {tilde over (v)} along direction 1 inside the hyper-rectangle share the same k*, which is a consequence of Lemma 1. In FIG. 4, the colored tiles for each region are drawn on the surfaces of the hyper-rectangle, and regions extend unchanged parallel to the [1,1,1]^T direction, which is the direction of the line of sight in FIG. 4. PRoM detects wrapping integers for fast, accurate estimation of velocity.

In addition, ({tilde over (v)}, v_k) as a function of {tilde over (v)} is piece-wise quadratic with the same curvature for all k, so the decision boundaries of k*({tilde over (v)}) are linear, which is also illustrated in FIG. 4.

By Lemma 1, the error in wrapping integers, (k*({tilde over (v)})−k({tilde over (v)})_h, remains constant for all noise realizations n along any line parallel to 1. So, the space of all possible noise realizations,

${ℝ}^{\frac{N_e\left(N_e-1\right)}{2}},$

can be divided into “tubes” (x) parallel to 1, based on the difference, x, of estimated wrapping integers k* and true wrapping integers k.

$\begin{array}{cc}\left(x\right)=\left\{n|{\langle k^{\star }\left(\tilde{v}\right)-k\left(\tilde{v}\right)\rangle }_h=x\right\}.& \left(43\right)\end{array}$

As seen below, integration over these tubes can be performed to arrive at error probabilities for detecting the wrapping integers. Next, a second lemma is established which describes the orthogonality between pairwise noise differences and the error in the estimated velocity.

$\mathrm{Lemma}2$ $\begin{array}{cc}{n}_{\mathrm{ab}}-n_{cd}\perp w^{T}n,\mathrm{ab}\ne cd& \left(44\right)\end{array}$ $\mathrm{Proof}.$ $\mathrm{cov}\left(w^{T}n,n_{\mathrm{ab}}-n_{cd}\right)=w^T𝔼\left(n\left(n_{\mathrm{ab}}-n_{cd}\right)\right)=\frac{1^T{\sum }^{-1}\left(n\right)\sum \left(n\right)\left(e_{\mathrm{ab}}-e_{cd}\right)}{1^T\sum {\left(n\right)}^{-1}1}=0,$

• • where e_ab is the standard basis equal to 1 at one position corresponding to n_ab in n, and 0 otherwise.

With these two lemmas, the distribution can be specified.

Theorem 1. Given v, ∀nϵ(x), {circumflex over (v)}({tilde over (v)}) follows wrapped normal distribution (v+w^T(x·2venc)_Ω, w^TΣ(n)w).

Proof. Note that the event nϵ(x) is only determined by the pairwise difference of n_ab−n_cd thus uncorrelated, thereby independent of w^Tn by the Gaussian assumption (25). So, for all nϵ(x)

$\begin{array}{c}\hat{v}\left(\tilde{v}\right)={\langle w^T\left(v1+x·2\mathrm{venc}+n\right)\rangle }_{\Omega }\\ ={\langle v+w^T\left(x·2\mathrm{venc}\right)+w^{T}n\rangle }_{\Omega }\end{array}.$

Let f({circumflex over (v)}|(x), v) denote the conditional Gaussian probability density function (pdf) in Theorem 1 without wrapping by modulo Ω. Thus, by wrapping the pdf function and invoking the law of total probability, results in

$\begin{array}{cc}\begin{array}{c}\mathbb{P}\left(\widehat{v}|v\right)={\sum }_x\mathbb{P}\left(n\in \left(x\right)|v\right)f_{\Omega }\left(\hat{v}|\left(x\right),v\right)\\ f_{\Omega }\left(\hat{v}|\left(x\right),v\right)=\{\begin{array}{cc}{\sum }_{l\in \mathbb{Z}}f\left(\hat{v}+l\Omega |\left(x\right),v\right),& \hat{v}\in [0,\Omega )\\ 0,& \mathrm{otherwise}\end{array}\end{array}.& \left(45\right)\end{array}$

To complete (45), (nϵ(x)|v) need to be calculated by integration of a multivariate normal distribution. Leveraging Lemma 1, the integration can be simplified to a definite integration of a normal distribution in

$\frac{N_e\left(N_e-1\right)}{2}-1$

variables.

FIGS. 5A-5B illustrate the histogram of {circumflex over (v)}|v=0 using 10⁵ trials, N_e=3, N_c=1, venc=[99,18,22]^T cm/s at s₂₁=5,10. Invoking Theorem 1, the five most probable wrapping integers are predicted, which correspond to detection regions ([0,0,0]_T), ([5,4,1]), ([6,5,1]^T), ([1,1,0]^T), and ([10,8,2]^T). These five predicted detections correspond to f(({circumflex over (v)}|(x), v) centered at 0 (the true velocity), ±177.81, and ±40.37 cm/s; these five predictions are validated in the histogram. For the higher SNR case of s₂₁=10, the probability of unwrapping errors goes very small, and the 10⁵ trials are insufficient to encounter an unwrapping error to ±40.37 cm/s. In FIG. 5, the logarithmic scale for the vertical axis covering four (FIG. 5A) and five (FIG. 5B) orders of magnitude. The histogram illustrates both noise sensitivity via the spread of each Gaussian component and the probability of unwrapping errors via the presence of multiple components.

F. Three-point Encoding

Below, a three-point encoding for velocity in one direction is described. Due to (9, 14), every venc_ab, and hence the unambiguous range, Ω, depends only on the differences between first moments; thus, venc_ab and Ω are unaffected by adding the same constant to every first moment. For the first moments it is assumed:

$\begin{array}{cc}{m}_{11}<m_{12}<m_{13},m_{12}-m_{11}<m_{13}-m_{11}.& \left(46\right)\end{array}$

Thus, the three vencs are ordered:

${\mathrm{venc}}_{31}<ven{c}_{32}<ven{c}_{21}.\mathrm{Let}\xi =\frac{{\mathrm{venc}}_{32}}{{\mathrm{venc}}_{31}},$

noting that (9) and (46) imply ξϵ(1,2). For rational ξ=p/q with co-prime integers p and q, the unambiguous range Ω in (14) is

$\begin{array}{cc}\Omega =2\left(p-q\right){\mathrm{venc}}_{21}.& \left(47\right)\end{array}$

Thus, by jointly unwrapping multiple vencs one can construct an unaliased velocity range that is larger than the highest venc, venc₂₁, by a factor of 2(p−q). The covariance matrix for three-point encoding can be formulated via (24). Correspondingly, the combination weights w can be calculated and the resulting RMSE given wrapping integers (26, 27). Armed with the explicit error variance and the probability of unwrapping errors derived below, an optimized design of venc is presented in II-H for three-point encoding with the constraint on the largest first moment, given a desired unambiguous range of velocities to be observed and SNR level for each encoding.

G. Existing Estimators for Three-point Encoding

In the notation of (10) and (13), the unaliased high venc measurement, {tilde over (v)}₂₁, to unwrap the low venc measurement, {tilde over (v)}₃₁, is used while venc₃₂ goes unused. The estimator in, is denoted as standard dual-venc (SDV), is given by

$\begin{array}{cc}\hat{v}=\{\begin{array}{cc}{\tilde{v}}_{21}-4{\mathrm{venc}}_{31},& \frac{{\tilde{v}}_{21}-{\tilde{v}}_{31}}{2{\mathrm{venc}}_{31}}\in \left(-2.4,-1.6\right)\\ {\tilde{v}}_{21}-2{\mathrm{venc}}_{31},& \frac{{\tilde{v}}_{21}-{\tilde{v}}_{31}}{2{\mathrm{venc}}_{31}}\in \left(-1.2,-0.8\right)\\ {\tilde{v}}_{21}+2{\mathrm{venc}}_{31},& \frac{{\tilde{v}}_{21}-{\tilde{v}}_{31}}{2{\mathrm{venc}}_{31}}\in \left(0.8,1.2\right)\\ {\tilde{v}}_{21}+4{\mathrm{venc}}_{31},& \frac{{\tilde{v}}_{21}-{\tilde{v}}_{31}}{2{\mathrm{venc}}_{31}}\in \left(1.6,2.4\right)\end{array}.& \left(48\right)\end{array}$

Two potentially aliased measurements {tilde over (v)}₃₁, {tilde over (v)}₃₂ are jointly unwrapped by minimizing

$\begin{array}{cc}\hat{v}=\underset{v\in [-\Omega /2,\Omega /2)}{\arg \min }{\sum }_{l=\left(31,32\right)}\left(1-\cos \left(\frac{\pi v}{{\mathrm{venc}}_l}-{\tilde{\theta }}_l\right)\right).& \left(49\right)\end{array}$

This prior approach is called “optimal dual-venc (ODV)” herein, and the minimization is accomplished by searching vϵ[−Ω/2, Ω/2) with grid spacing

$\frac{{\mathrm{venc}}_{31}}{1000}.$

The cost adopted in ODV is equivalent to

$\begin{array}{cc}\frac{1}{2}{❘e^{i\frac{\pi v}{{\mathrm{venc}}_{31}}}-e^{i{\tilde{\theta }}_{31}}❘}^2+\frac{1}{2}{❘e^{i\frac{\pi v}{{\mathrm{venc}}_{32}}}-e^{i{\tilde{\theta }}_{32}}❘}^2,& \left(50\right)\end{array}$

which intrinsically assumes no correlation between the noisy phase differences, {tilde over (θ)}₃₁ and {tilde over (θ)}₃₂. The ODV approach recommends

${\mathrm{venc}}_{32}=\frac{3}{2}{\mathrm{venc}}_{31},$

yielding an unambiguous velocity range of length 2venc₂₁, which is twice the highest venc. The choice 3/2 is a heuristic to lessen the probability of unwrapping errors when minimizing (49) in the presence of noise.

The non-convex optimization (NCO) in prior work iteratively minimizes a cost similar to (50) with weights to accommodate a lower SNR in presence of intravoxel dephasing:

$\begin{array}{cc}{❘{\tilde{r}}_{31}❘}^2{❘e^{i\frac{\pi v}{{\mathrm{venc}}_{31}}}-e^{i{\tilde{\theta }}_{31}}❘}^2+{❘{\tilde{r}}_{32}❘}^2{❘e^{i\frac{\pi v}{{\mathrm{venc}}_{32}}}-e^{i{\tilde{\theta }}_{32}}❘}^2.& \left(51\right)\end{array}$

Both ODV and NCO can be applied to any number of encodings; further, the NCO algorithm also incorporates a spatial regularization across voxels in the form of the Laplacian of the velocity map.

H. Design for Three-Point Encoding

The selection of the pair {venc₃₁, ξ} defines first moments [m₁₁, m₁₂, m₁₃] up to translation. To minimize the worst intravoxel dephasing, symmetric encoding m₁₁=−m₁₃ may be selected; to further ameliorate intravoxel dephasing, a user-defined upper bound on the largest first moment is provided for, m₁₃≤m_τ. The choice of venc entails the interplay of noise sensitivity, probability of correct unwrapping, and the range of reliably unaliased velocities. A performance guarantee strategy is adopted for navigating these competing objectives. The design inputs are: the maximum range of velocities to be reliably detected, Ω_ϵ; a lower bound of operating measurement SNR, s_αβ; an upper bound, m_τ, on the magnitude of the largest first moment; and bounds on the per-pixel probability of an unwrapping error. The design outputs are the venc and an underlying [m₁₁, m₁₂, m₁₃]. To formalize the notion of optimality, four definitions may be made.

• • Definition 1 (Unwrapping error). If k*({tilde over (v)})−k({tilde over (v)})_h≠0, i.e., wrapping integers are incorrectly detected, then an Unwrapping Error occurs.
  • Definition 2 (Aliasing error) Given no unwrapping error, if {circumflex over (v)}_k* is aliased, then an Aliasing Error occurs.
  • Definition 3 (ϵ-Reliable Range) For given measurement SNR, s_αβ, and two small numbers ϵ₁, ϵ₂, the e-reliable range Ω_ϵ is the range of the velocities for which (Unwrapping Error)≤ϵ₁ and (Aliasing Error)≤ϵ₂.
    • Definition 3 allows specifying a reliable velocity range Ω_ϵ<Ω to guard against aliasing. Armed with these three definitions, a precise meaning of optimality for three-point design can be stated.
  • Definition 4 (Optimal venc Design) Given the SNR for the complex-valued data s_αβ, the desired maximum velocity range of length Ω_ϵ, an upper bound m_τ on the largest first moment, and unwrapping error bounds ϵ₁, ϵ₂, the optimal venc minimizes the RMSE among all designs for which the unwrapping and aliasing errors satisfy (Unwrapping Error)≤ϵ₁ and (Aliasing Error)≤ϵ₂ across the entire range of length Ω_ϵ.

Given s_αβ and venc, (Unwrapping Error) can be calculated through Monte Carlo simulation by setting v=0 and counting the trials for which k*−k_h≠0. Independence of (x) on v allows simulation of v=0 to be sufficient. The number of trials is selected as 100/ϵ₁.

Bounding the Aliasing Error can be achieved via explicitly designing Ω different than Ω_ϵ. Let the normal cumulative distribution function be Φ(·). Then (Aliasing Error)≤ϵ₂ when

$\begin{array}{cc}\Omega -{\Omega }_{ϵ}\ge 2{\Phi }^{-1}\left(1-{ϵ}_2\right)\sqrt{w^T\sum \left(n\right)w}.& \left(52\right)\end{array}$

The design procedure in Alg. 2 is an offline finite search to optimally design venc and the corresponding first moments. To explore design options for ξϵ(1,2), rational values ξ=p/q are searched among

$\begin{array}{cc}\left\{\frac{p}{q}❘\gcd \left(p,q\right)=1,\frac{p}{q}\in \left(1,2\right),p\le P,q\le Q\right\},& \left(53\right)\end{array}$

where gcd(·,·) is greatest common divisor, and P, Qϵ⁺ are predefined upper bounds on the positive integers p, q.

The output [m₁₁, m₁₂, M₁₃] of Alg. 2 is a symmetric encoding that can be translated by ±m₁₃ to yield a referenced encoding; however, the referenced encoding suffers an increased risk of severe intravoxel dephasing, owing to the doubling of the largest first moment.

I. Post-Processing Using Spatial Information: PROM+

Below is presented an effective post-processing strategy paired with PRoM. The PRoM per-voxel estimation can benefit from leveraging spatial correlations among the per-voxel phase unwrapping integers. It is assumed that the noiseless velocity map u(ρ) belongs to a surface class U where ρ denotes spatial position. For example, a polynomial model has been used for the brain image phase and the Hagen-Poiseuille equation has been used for laminar blood flow throughout most of the circulatory system.

When the model is accurate, the difference between the noisy unbiased estimated and true velocity map at each location should be at the noise level, and we assume at each location the difference follows i.i.d. normal distribution with variance 1/(2λ).

From the negative log likelihood, the spatial post-processing can be expressed as minimizing the negative log likelihood

$\begin{array}{cc}\underset{u\in U,\hat{v}}{\arg \min }ℒ\left(\tilde{v}\left(\rho \right),\hat{v}\left(\rho \right)\right)+{\lambda \left(\hat{v}\left(\rho \right)-u\left(\rho \right)\right)}^2& \left(54\right)\end{array}$

Here, to avoid over-smoothness due to the regularization using uϵU, {circumflex over (v)}(ρ)ϵ{{circumflex over (v)}_k|kϵK({circumflex over (v)}(ρ))} may be restricted, i.e., only allow spatial information to affect k.

To optimize this spatially regularized cost, an alternating minimization strategy may be adopted. For current {circumflex over (v)}(ρ), it may be fit with the best u(ρ) via surface fitting. For current u(ρ), the choice of {circumflex over (v)}(ρ) is updated, per voxel to minimize the cost. These two steps guarantee convergence in terms of the cost function. Iterations continue to convergence, which for the integer-valued k simply means no change; no convergence threshold is required, as would be with real-valued variables. Convergence is observed in two iterations in all experiments reported below. To reduce computation, only a few most likely velocity candidates are considered. Moreover, air regions are masked-out through magnitude thresholding to reduce computation.

The PRoM estimator, together with the spatial post-processing, is denoted “PRoM+”. In the section below, U is adopted for a basic non-parametric local quadratic regression for both phantom and in vivo experiments.

III. Methods

FIG. 6 illustrates an example method 600 of implementing PRoM for N_e-point encoding. The method 600 is a fast algorithmic procedure to estimate through-plane velocity at a voxel from three or more encoding gradients. At 602, the proposed estimator first constructs a set of candidate tuples of wrapping integers; the probability that the true tuple of wrapping integers is in this set is arbitrarily close to one. At 604, at each candidate tuple of wrapping integers, an approximate conditional maximum likelihood estimator performs a linear combination that accounts for noise variance and correlation in the relative phases of image voxels across all encodings. At 606, the estimator provides an unbiased, grid-free alternative to searching a dense grid of possible velocity estimates.

A. Simulation

To validate the two assumptions (25, 29) used in PRoM, the RMSE values for PRoM and for grid search MLE are compared to the square root of the Cramér-Rao lower bound (CRLB) derived from complex measurements in (3). Results are computed for venc=[15,6,10]^T cm/s and 50% intravoxel dephasing of amplitudes for high first moments: 2s₁₁=s₂₁=2s₃₁. RMSE values for both the MLE from the complex measurements and PRoM are each calculated using 10⁵ trials, where the grid search of MLE on v has spacing 0.006 cm/s to avoid bias from gridding.

Algorithm 2 PRoM Optimal Design for Three-Point Encoding
Require: P, Q, sαβ, Ω  , γ1, γ2, mr
1:  σ ← ∞, p ← 2
2:  while p ≤ P do
3:  q ← 1
4:  while q ≤ Q do
5:  if gcd(p, q) = 1 and p/q ϵ (1, 2) then
6:  v ← 0, venc ← [pq, q(p − q), p(p − q)].
7:  Simulate {tilde over (v)} with 100γ1−1 trials, apply PRoM.
8:  if (Frequency of (k*({tilde over (v)})−k({tilde over (v)}))h ≠ 0) < γ1 then
9:  c⟵ΩϵΩ-2⁢Φ-1(1-γ2)⁢wT⁢∑(n)⁢w
10: c⟵max[c,π2⁢γ⁢mγ⁢q⁡(p-q)]T
11: if c{square root over (wTΣ(n)w)} ≤ σ then
12: venc ← evene, σ ← c{square root over (wTΣ(n)w)}
13: end if
14: end if
15: end if
16: q ← q + 1
17: end while
18: p ← p + 1
19: end while
20: m                        11                                            ⟵                                                                        -                          π                                                                          2                          ⁢                          γ                          ⁢                                                      venc                            31                                                                                                                ,                                                                                            m                          12                                                ⟵                                                  π                                                      2                            ⁢                            γ                            ⁢                                                          venc                              31                                                                                                                          -                                              π                                                  2                          ⁢                          γ                          ⁢                                                      venc                            32                                                                                                                ,                                                                  m                        13                                            ⟵                                              -                                                  m                          11
Ensure: venc, [m11, m12, m13]
indicates data missing or illegible when filed

To compare the performance of PRoM versus SDV and ODV, the same measurements are processed using different estimators. Simulation parameters include: venc=[15,6,10]^T cm/s, [s₁₁, s₂₁, s₃₁]=[10,20,10], N_c=1 coil. The ODV estimation algorithm uses {tilde over (v)}₃₁, {tilde over (v)}₃₂, while SDV uses {tilde over (v)}₃₁, {tilde over (v)}₂₁. RMSE is calculated averaged over 10⁵ trials at each true v on [−30,30] cm/s sampled every 0.1 cm/s.

To assess the optimized encoding design in Alg. 2, a required velocity range of [−150,150] cm/s is set and the recommended choices of symmetric three-point encoding are considered for each algorithm. Simulation parameters include: N_c=1 coil, s₂₁=20. SDV suggests using venc=[150,60,100]^T cm/s. ODV recommends

$\xi =\frac{3}{2}$

and specifies the three vencs to be venc=[150,50,75]^T cm/s. For PRoM, it is assumed intravoxel dephasing leads to [s₁₁, s₂₁, s₃₁]=[10,20,10]; other input includes [P,Q]=[10,10], [ϵ₁, ϵ₂]=[10⁻⁷, 10⁻⁷], Ω_ϵ=300 cm/s,

$\gamma m_{\tau }=\frac{\pi }{50}s/\mathrm{cm}.$

The design procedure in Alg. 2 results in venc=[153.73, 25.62, 30.75]^T cm/s. Because PRoM uses a 95.2% larger m₁₃ thereby potentially leading to more intravoxel dephasing, ODV and SDV are advantaged by assuming no intravoxel dephasing: [s₁₁, s₂₁, s₃₁]=[20,20,20]. RMSE is calculated averaged over 10⁵ trials at each true v on [−150,150] cm/s sampled every 0.5 cm/s.

To simulate the complex intravoxel dephasing and assess estimator performance in this case, vessels are simulated with circularly symmetric parabolic velocity profiles, as seen in prior works. Parameters include: 0.1 mm³ isotropic resolution and N_c=1 coil. The five vessels share the same peak velocity 60 cm/s but have different diameters of 5.5,3.9,3.2,2.7,2.4 mm. The proton density is set to be 30% in the background region and 50% in the static tissue region compared to that in the vessel regions. The complex signal is generated using (3) with symmetric three-point encoding such that venc=[60,20,30]^T cm/s. Regions of 5×5×5 voxels are merged into one to generate intravoxel dephasing and 0.5 mm³ isotropic resolution. Then i.i.d. white complex Gaussian noise is added to make the maximum s_αβ for all voxels in the vessel regions reach 30. No post-processing is used with PRoM to allow for pure comparison of per-voxel estimation performance in various dephasing scenarios.

B. Phantom

A phantom experiment allows for a controlled comparison of estimation performance. An agarose gel-filled cylindrical container is used to generate the MRI signal. An air-coupled propeller rotates the container. The rotational rate is counted with a photomicrosensor. The phantom was scanned on a 1.5T scanner (Siemens MAGNETOM Avanto). The in-plane bottom to top velocity increases linearly with the distance from the center of the container. In this experiment, the vertical component of the in-plane velocity is encoded using a symmetric three-point acquisition, and the velocity component ranges from −240 to 240 cm/s. The acquisition parameters include:

$N_c=16\mathrm{coils};\mathrm{venc}={\left[200,100,\frac{500}{3}\right]}^T\mathrm{cm}/s;$

field-of-view (FOV) 520×260 mm; flip angle 5°; TR 4.38 ms; TE 2.66 ms; and, matrix size 192×125. There are 15 repeated acquisitions. The averaged k-space is used as a reference to calculate the RMSE and aliasing error for all voxels except the background region across 15 scans. For per-frame post-processing of PRoM, voxels with less than 30% maximum voxel magnitude are masked, and only the two most likely velocity estimates are considered. Locally weighted quadratic surface class, U, is adopted using a span of 25% closest points and λ=1.

C. In Vivo

An in vivo experiment is used to verify that PRoM can unwrap velocity on an interval Ω larger than twice the largest venc, venc₂₁, as claimed in (47). A healthy volunteer was scanned on a 3T scanner (Siemens MAGNETOM Vida). For the recruitment and consent of human subject used in this study, the ethical approval was given by an Internal Review Board (2005H0124) at The Ohio State University. The venc scouting scan showed the maximum absolute value of velocity to be above 90 cm/s and less than 100 cm/s. A breath-held, N_c=30 coils, symmetric three-point encoding PC-MRI dataset was collected, with the imaging plane intersecting both the ascending aorta and descending aorta; through-plane velocity was encoded. Other acquisition parameters include: FOV 360×270 mm; flip angle 15°; TR 5.56 ms; TE 3.69 ms; matrix size 192×108; and, cardiac phases, 13. The three-point acquisition is designed using Alg. 2 for: [P,Q]=[10,10], [s_1β, s_2β, s_3β]=[5,10,5], [ϵ₁, ϵ₂]=[10⁻⁶, 10⁻⁶], Ω_ϵ=200 cm/s,

$\gamma m_{\tau }=\frac{\pi }{40}s/\mathrm{cm}.$

The design results in ξ=5/3 with highest first moment

$\frac{\pi }{\gamma m_{13}}=42\mathrm{cm}/s.$

Restricted by input precision number of the scanner interface, venc=[52.5,21,35]^T cm/s; the resulting unambiguous range is Ω=210 cm/s, which is double the range [−52.5,52.5) cm/s of SDV processing. For per-frame post-processing of PRoM, after square-root sum-of-squares coil combination, voxels less than 30% maximum image were masked, and only the two most likely velocity estimates were considered. Due to the complex velocity map across the FOV, locally weighted quadratic surface class, U, was adopted using a span of 3% closest points and λ=1.

IV. Results

A. Simulation Results

FIG. 7 numerically explores the approximations adopted in PRoM. The square root of the CRLB is plotted versus s₂₁ for the model in (3) and provides a lower bound on the RMSE for any unbiased estimator. The bound is from a local analysis of the likelihood function, and thus optimistically does not consider unwrapping errors. Superimposed are the RMSE values for both the MLE from the complex measurements and PRoM. For low SNR, both estimators diverge from the bound due to unwrapping errors. Both the MLE and PRoM asymptotically coincide with CRLB^0.5. Moreover, throughout the entire range of noise powers considered, the RMSE difference between MLE and PRoM is negligible, hence justifying the characterization of PRoM as an approximate MLE and validating the assumptions adopted in the derivation of the PRoM estimator.

FIGS. 8A-8B show the RMSE results for the same acquisition as in FIG. 7 with 10⁵ trials at each true velocity value with grid 0.1 cm/s. Both ODV and PRoM can unwrap a large range of velocities. FIG. 8B provides a zoomed view. PRoM improves RMSE by modeling the non-zero noise correlation between phase difference measurements.

FIGS. 9A-9B illustrate the benefit of optimized design. Results are computed for the same acquisition as in FIG. 7 with 10⁵ trials at each true velocity value with grid 0.1 cm/s. Here, ODV and SDV are advantaged by assuming no intra-voxel dephasing for their acquisition. For PRoM, intra-voxel dephasing is present with 2s₁₁=s₂₁=2s₃₁. The PRoM design uses ξ=6/5, which explicitly suppresses both Unwrapping Errors and Aliasing Errors to ensure reliable estimation across the full range, [−150,150] cm/s. Further, despite the handicap of simulated 50% intravoxel dephasing, PRoM still provides a 10.5% decrease in RMSE compared to ODV and a 25.1% decrease versus SDV used with its recommended acquisition strategy.

FIGS. 10A-10H show the result for simulation of intravoxel dephasing. FIG. 11A is the velocity profile simulated with refined resolution and FIG. 11B is the lower acquired resolution which leads to intravoxel dephasing. FIG. 11C, FIG. 11D and FIG. 11E illustrate intravoxel dephasing for m₁₁, m₁₂, m₁₃. More serious intravoxel dephasing is observed in m₁₁, m₁₃ and at voxels close to the boundaries of the simulated vessels. FIG. 11F, FIG. 11G, and FIG. 11H show the recovered velocity. Aliasing error is observed in the SDV estimated velocity, but not in ODV and PRoM. The RMSE values for SDV, ODV and PRoM are 120.27, 10.24 and 10.12 cm/s, respectively. Moreover, the RMSE given no Aliasing Error for SDV, ODV, and PRoM are 10.23, 10.24, and 10.12 cm/s, respectively.

B. Phantom Results

FIGS. 11A-11K use measured data from a spinning phantom to evaluate RMSE in velocity estimation. In addition, the phantom data validate that both ODV and PRoM paired with simple post-processing can reliably unwrap a larger range of velocities than SDV, given the same encodings. The number of aliased voxels and RMSE values are reported in Table 1. The PRoM+ iteration is observed for all frames to converge in only two iterations. In this instance, PRoM+ eliminates aliasing errors and reduces RMSE by 25.8% versus ODV, and 48.5% compared to SDV.

TABLE 1
Metric	SDV	ODV	PRoM	PRoM+
Number of aliased voxels	27	5	5	0
RMSE, all voxels (cm/s)	6.08	4.22	3.85	3.13
RMSE, excluding aliasing (cm/s)	3.22	3.59	3.13	3.13

C. In Vivo Results

FIGS. 12A-12K show in vivo results. FIGS. 12A-12C show more serious intra-voxel dephasing for m₁₁ and m₃ compared to m₁₂. Because the largest absolute value of true velocity is larger than the largest venc of 52.5 cm/s, significant aliasing is observed in SDV recovery from FIG. 12E. However, FIG. 12E and FIG. 12F illustrate that both ODV and PRoM can recover velocities larger than the largest venc. Here, the acquisition designed for PRoM departs from the ξ=3/2 heuristic to use ξ=5/3, resulting in an unambiguous range of velocities four times that for standard dual-venc for the same highest venc. FIG. 12G illustrates that PRoM+ can incorporate spatial correlations to improve unwrapping performance. For the in vivo example using coil-combined image: ODV computation takes 9.106 seconds, and PRoM takes 2.246 seconds, with an additional 1.717 seconds for PRoM+. The iteration of PRoM+ is observed to converge in only two steps for all frames.

V. Discussion

For the simulation and phantom studies, where the ground truth is available, PRoM offers a significant RMSE advantage over standard dual venc processing, i.e., 25.1% for the simulation study and 48.5% for the phantom study. Although PRoM offers a fourfold computation advantage over ODV, its RMSE advantage over ODV, when both methods use the same venc values, is marginal. However, there are two features that distinguish PRoM from ODV, NCO, and other dual-venc methods. First, PRoM allows for an optimized venc design, which can translate to a more significant reduction in RMSE, as evident by 10.5% reduction in RMSE over ODV (FIG. 9).

Second, PRoM can leverage the conditional probabilities of different wrapping integers to enable a new mechanism for phase unwrapping.

The presentation here for PRoM is limited to a single component of velocity; extension to the estimation of all three velocity components, and hence congruence equations in multiple variables, is ongoing work.

Several three-point encoding have been proposed and validated for PC-MRI aiming to improve VNR or unambiguous velocity range. Although performance depends strongly on vencs, the selection of vencs has been based on heuristics. PRoM, for the first time, provides an avenue to optimize vencs. FIGS. 12A-12K demonstrates that the PRoM-inspired design procedure in Alg. 2 can provide an unambiguous velocity range more than four times as large as the highest venc. The acquisition in FIGS. 12A-12K illustrates the derivation in (47) and the associated design procedure in Alg. 2. Indeed, the design procedure allows the range to grow to the greatest extent allowed by the presumed noise floor, which is given as an input to the design. The design then minimizes the predicted RMSE while meeting constraints on unwrapping errors, reliable range of velocities, and maximum first moment. In the regime of low SNR, the PRoM design yields (=3/2, coinciding with a conventional heuristic. The Alg. 2 design provides unwrapping and velocity range guarantees, given a noise floor and bound on the highest first moment. For any higher SNR encountered, the guarantee still holds, and the RMSE reduces according to (27).

For volumetric imaging applications with vast numbers of voxels, the processing speed of PRoM may provide a desirable benefit. The careful pruning of the set of candidate wrapping integers results in a fast estimator without expensive grid search or gradient-based iterative optimization.

PRoM+ provides a simple but effective post-processing strategy for PRoM, which only affects the choice of k from a Bayesian perspective. It differs from conventional unwrapping algorithms that typically only allow ±2π adjustment in the possibly wrapped phases. And, the processing incorporates the relative conditional probabilities of different wrapping integers, which are available as a byproduct of the PRoM algorithm. This strategy can be easily generalized to multiple velocity components and high-dimensional imaging.

FIG. 13 is a view illustrating a structure of an example magnetic resonance imaging (MRI) apparatus 1300 that may be used to acquire image data. Generally, the MRI apparatus 1300 is a low-field or high-field system equipped with a high-end data processing unit (DPU) to enable the implementation of structure aware (SA) recovery in clinically relevant times. The DPU may be comprised of multiple GPUs running, e.g., a Gadgetron framework, which was developed at the National Heart, Lung, and Blood Institute. The SA recovery methods are based on Bayesian inference where the structure in the image is expressed probabilistically, as a set of priors or conditional priors. For implementation of SA recovery, fast iterative image recovery methods based on Generalized Approximate Message Passing (GAMP) may be employed, yielding fast algorithms with self-tuning parameters.

The MRI apparatus 1300 includes a scanner 1303 that generates magnetic fields used for the MR examination within a measurement space 1304 having a patient table 1302. In accordance with the present disclosure, the scanner 1303 may include a wide bore 70 cm superconducting magnet having a field strength of approximately 0.5-3.0 Tesla (T).

A controller 1306 includes an activation unit 1311, a receiver device 1312 and an evaluation module 1313. During a phase-sensitive flow measurement, MR data are recorded by the receiver device 1312, such that MR data are acquired in, e.g., a measurement volume or region 1315 that is located inside the body of a patient 1305. The MRI apparatus 1300 may: include a coil array (e.g., arranged as two 3×3 grids); support parallel imaging using SPIRiT, GRAPPA, SENSE, VISTA, AMP, FISTA, SCoRE, and/or Bayesian Inference; and perform analog-to-digital conversion (ADC) at a gantry of the MRI apparatus 1300.

An evaluation module 1313 prepares the MR data such that they can be graphically presented on a monitor 1308 of a computing device 1307 and such that images can be displayed. In addition to the graphical presentation of the MR data, a three-dimensional volume segment to be measured can be identified by a user using the computing device 1307. The computing device may include a keyboard 1309 and a mouse 1310. The computing device may include a Xeon central processing unit (CPU) or better; 16 GB of random-access memory (RAM); Multi-GPU, K20 or Titan Z reconstruction hardware; support DiCOM 3.0; and support simultaneous scan and reconstruction.

Software for the controller 1306 may be loaded into the controller 1306 using the computing device 1307. Such software may implement a method(s) to process data acquired by the MRI apparatus 1300, as described below. It is also possible for the computing device 1307 to operate such software. Yet further, the software implementing the method(s) of the disclosure may be distributed on removable media 1314 so that the software can be read from the removable media 1304 by the computing device 1307 and be copied either into the controller 1306 or operated on the computing device 1307 itself.

Although the subject matter has been described in language specific to structural features and/or methodological acts, it is to be understood that the subject matter defined in the appended claims is not necessarily limited to the specific features or acts described above. Rather, the specific features and acts described above are disclosed as example forms of implementing the claims. The present disclosure is capable of other implementations and of being practiced or carried out in various ways.

It must also be noted that, as used in the specification and the appended claims, the singular forms “a,” “an” and “the” include plural referents unless the context clearly dictates otherwise. Ranges may be expressed herein as from “about” or “approximately” one particular value and/or to “about” or “approximately” another particular value. When such a range is expressed, other exemplary implementations include from the one particular value and/or to the other particular value.

By “comprising” or “containing” or “including” is meant that at least the named compound, element, particle, or method step is present in the composition or article or method, but does not exclude the presence of other compounds, materials, particles, method steps, even if the other such compounds, material, particles, method steps have the same function as what is named.

In describing example implementations, terminology will be resorted to for the sake of clarity. It is intended that each term contemplates its broadest meaning as understood by those skilled in the art and includes all technical equivalents that operate in a similar manner to accomplish a similar purpose. It is also to be understood that the mention of one or more steps of a method does not preclude the presence of additional method steps or intervening method steps between those steps expressly identified. Steps of a method may be performed in a different order than those described herein without departing from the scope of the present disclosure. Similarly, it is also to be understood that the mention of one or more components in a device or system does not preclude the presence of additional components or intervening components between those components expressly identified.

As discussed herein, a “subject” (or “patient”) may be any applicable human, animal, or other organism, living or dead, or other biological or molecular structure or chemical environment, and may relate to particular components of the subject, for instance specific organs, tissues, or fluids of a subject, may be in a particular location of the subject, referred to herein as an “area of interest” or a “region of interest” (ROI).

Claims

1. A method for Phase Recovery from Multiple Wrapped Measurements (PRoM), comprising: constructing a set of candidate tuples of wrapping integers;at each candidate tuple of wrapping integers, performing a linear combination using an approximate conditional maximum likelihood estimator; andproviding an estimate of through-plane velocity at an image voxel.

2. The method of claim 1, wherein a probability that a true tuple of wrapping integers is in this set is arbitrarily close to one.

3. The method of claim 1, wherein the approximate conditional maximum likelihood estimator accounts for noise variance and correlation in the relative phases of image voxels across all encodings.

4. The method of claim 1, further comprising producing an efficiently pruned list of candidate wrapping integers.

5. The method of claim 1, further comprising producing a spatially adaptive, data-driven estimator of the scaled noise covariance matrix for phase differences.

6. The method of claim 5, wherein the method is performed without reliance on direct measurement of per-voxel noise power.

7. The method of claim 5, further comprising incorporating effects of spin dephasing.

8. The method of claim 1, further comprising producing a posterior probability distribution of the velocity estimate at each voxel.

9. The method of claim 8, further comprising an optimized design of data acquisition to minimize mean-squared estimation error subject to constraints on at least one of minimum range of unaliased velocities; upper bound on first moment of the time-varying magnetic field gradient; upper bound on per-voxel probability of unwrapping error.

10. The method of claim 8, further comprising spatial post-processing to reduce bias due to phase unwrapping errors.

11. An MRI apparatus, comprising: a scanner that generates magnetic fields used for the MR examination;a measurement space having a patient table; anda controller having an evaluation module,wherein the evaluation module executes instructions for performing a method for Phase Recovery from Multiple Wrapped Measurements (PRoM) that includes: constructing a set of candidate tuples of wrapping integers;at each candidate tuple of wrapping integers, performing a linear combination using an approximate conditional maximum likelihood estimator;spatially post-processing to further remove bias from phase unwrapping errors; andproviding an estimate of through-plane velocity at each voxel in an array of voxels.

12. The apparatus of claim 11, wherein a probability that a true tuple of wrapping integers is in this set is arbitrarily close to one.

13. The apparatus of claim 11, wherein the approximate conditional maximum likelihood estimator accounts for noise variance and correlation in the relative phases of image voxels across all encodings.

14. The apparatus of claim 11, wherein an efficiently pruned list of candidate wrapping integers is produced.

15. The apparatus of claim 11, wherein a spatially adaptive, data-driven estimator of the scaled noise covariance matrix for phase differences is produced.

16. The apparatus of claim 15, wherein the method executed by the apparatus is performed without reliance on direct measurement of per-voxel noise power.

17. The apparatus of claim 15, wherein effects of spin dephasing are incorporated.

18. The apparatus of claim 11, wherein a posterior probability distribution of the velocity estimate at each voxel is produced.

19. The apparatus of claim 18, wherein a design of data acquisition is optimized to minimize mean-squared estimation error subject to constraints on at least one of minimum range of unaliased velocities; upper bound on first moment of the time-varying magnetic field gradient; upper bound on per-voxel probability of unwrapping error.

20. The apparatus of claim 18, wherein spatial post-processing is used to reduce bias due to phase unwrapping errors.