Venc Design and Velocity Estimation for Phase Contrast MRI

A time-varying magnetic gradient field may be used to encode spin velocity into image phase. Here, we consider encoding the velocity component in one direction. Consider a spin moving through a magnetic field, the Taylor series expansion of spin position $p(t)\in\mathbb{R}$ at $t=0$ yields

where $v$ is instantaneous velocity and $a$ is acceleration. Let $\gamma$ be the gyromagnetic ratio, $B_{0}$ be the main static magnetic field strength, and $g(t)\in\mathbb{R}$ be the time-varying magnetic field gradient. Then to first order approximation of $p(t)$ [1], the phase accumulated from $t=0$ to echo time TE is

where $m_{0}$ and $m_{1}$ denote the zeroth and first moments of $g(t)$. The zeroth moment, $m_{0}$, encodes the spin position into phase and combines with $\gamma B_{0}\text{TE}$ to make the background phase, $\phi_{0}$. The first moment, $m_{1}$, 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

where $\alpha\in\{1,...,N_{e}\}$ indexes over all encodings, and $\beta\in\{1,...,N_{c}\}$ indexes over all coils. The resulting signal amplitude is $A_{\alpha}\in\mathbb{R}$; $S_{\beta}\in\mathbb{C}$ is the coil sensitivity; $v$ is the resultant instantaneous velocity; $m_{1\alpha}\in\mathbb{R}$ is the first moment; $\delta_{\alpha\beta}\in\mathbb{C}$ 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., intra-voxel dephasing effect). Generally, $A_{\alpha}$ decreases as $|m_{1\alpha}|$ increases [10].

Let $\widetilde{\bm{Y}}$ denote an ${N_{e}\times N_{c}}$ complex-valued measurement matrix with $(\alpha,\beta)^{\text{th}}$ entry $\widetilde{y}_{\alpha\beta}$ for encoding $\alpha$ and coil $\beta$. Observe that the unambiguous range of velocities for the $N_{e}$ encodings is

where $\text{LCM}(\cdot)$ 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+\Omega^{\prime}$ are indistinguishable given data, $\widetilde{\bm{Y}}$. Considering noise, the MLE of $v$ involves $N_{e}+2N_{c}+2$ real-valued unknowns, $\left\{v,\phi_{0},A_{1},\cdots,A_{N_{e}},|S_{1}|,\angle S_{1},\cdots,|S_{N_{c}}|,\angle S_{N_{c}}\right\}$, and is a nonlinear least-squares fit to the data. Here, $\angle(\cdot)$ denotes angle of a complex number. The optimization task can be reduced to

where $\widetilde{\bm{R}}\in\mathbb{C}^{N_{e}\times N_{e}}$ and the “steering vector” $\bm{s}\in\mathbb{C}^{N_{e}\times 1}$ are

Derivation of (5) is a straightforward extension of [11, p. 288] to accommodate unequal amplitudes, $A_{\alpha}$. Here, $(\cdot)^{\intercal}$ and $(\cdot)^{\mathsf{H}}$ denote transpose and conjugate transpose. For a given $v$, the amplitudes $A_{\alpha}$ may be found by solving (5) as a generalized eigenvalue problem [12, p. 461]; 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. Data are simulated using $\left[\frac{|A_{1}S_{1}|}{\sigma(\delta_{11})},\frac{|A_{2}S_{1}|}{\sigma(\delta_{21})},\frac{|A_{3}S_{1}|}{\sigma(\delta_{31})}\right]=[2.5,5,2.5]$ to represent a moderate 50% loss of signal amplitude due to intra-voxel dephasing. Throughout the paper, we use same 50% 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\approx-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.

From (5) we see that $\widetilde{\bm{R}}$ is a sufficient statistic for $v$ in (3). We depart from the multi-variate optimization in (5) to instead work with the phases of the off-diagonal entries of $\widetilde{\bm{R}}$, and we use amplitudes of $\widetilde{\bm{Y}}$ to construct the phase difference noise covariance matrix. In so doing, we gain four advantages: (1) characterization of the unambiguous set of estimated velocities; (2) characterization of the probability of unwrapping errors; (3) ability to design the encodings $\left[m_{11},...,m_{1N_{e}}\right]$ to minimize mean squared estimation error subject to guarantees on the probability of unwrapping error and required unambiguous range; (4) fast, grid-free parameter estimator for velocity, $v$. Further, the proposed surrogate estimator is asymptotically efficient, providing a good MLE approximation. In this section, we derive an approximate covariance matrix for the phase measurements in $\widetilde{\bm{R}}$ to lay the groundwork for building the proposed estimator of $v$.

Observe that $\widetilde{\bm{R}}=\widetilde{\bm{Y}}\widetilde{\bm{Y}}^{\mathsf{H}}$ performs coil combining and phase differencing. Denote the $(a,b)^{\text{th}}$ entry of $\widetilde{\bm{R}}$ as $\widetilde{r}_{ab}$, so

where $(\cdot)^{*}$ denotes complex conjugation. The noisy phase difference $\widetilde{\theta}_{ab}$ for encodings $a>b\in\{1,\cdots,N_{e}\}$ is

This phase differencing results in venc value

We have $\frac{N_{e}(N_{e}-1)}{2}$ such combinations. Throughout the paper, we let $\text{venc}_{ab}$ have units cm/s. Multiplying the vencs with phases, we obtain a (possibly wrapped) noisy velocity $\widetilde{v}_{ab}$

Phases $\widetilde{\theta}_{ab}$ are unambiguous on any interval of length $2\pi$, and for convenience we use $[0,2\pi)$; so, $\widetilde{v}_{ab}\in[0,2\text{venc}_{ab})$. Thus, we have a set of noisy congruence equations for $v$:

where $n_{ab}$ is additive zero mean noise. Define $k_{ab}$ as the wrapping integer for $\widetilde{v}_{ab}$. Momentarily assume that the true wrapping integer, $k_{ab}$, is known; then, we can rewrite (11) as a set of linear equations:

where $\circ$ is the Hadamard (element-wise) product, and

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

This also implies that $\Omega$ is the smallest repetition period of $v$ to construct the same $\widetilde{\bm{R}}$. Recall from (4) that $\Omega^{\prime}$ is a repetition period of $v$ to construct the same $\widetilde{\bm{Y}}$, as well as the same $\widetilde{\bm{R}}$. So, $\Omega^{\prime}$ is an integer multiple of $\Omega$.

Similar to the assumption $\widetilde{\theta}_{ab}\in[0,2\pi)$, we assume $v\in[0,\Omega)$ for convenience of derivation. This interval could also be shifted to $\left[-\frac{\Omega}{2},\frac{\Omega}{2}\right)$, for example, for bidirectional flow in MRI.

Let $\widehat{v}\in[0,\Omega)$ be the linear unbiased estimator of $v$ with smallest root mean squared error (RMSE). To compute $\widehat{v}$ from the noisy $\widetilde{\bm{v}}$ in (12), we only need the covariance matrix of the noise in the remainders, $\bm{\Sigma}\left(\bm{n}\right)$. Define $\sigma$ to be the standard deviation of the i.i.d. noise $\delta_{\alpha\beta}$ in (3). The mean and variance of $\widetilde{r}_{ab}$ are [13]:

Then from (15)-(16), the squared signal-to-noise ratio (SNR), or Rician factor, for $\widetilde{r}_{ab}$ is

where $s_{\alpha\beta}=\frac{|A_{\alpha}S_{\beta}|}{\sigma}=\operatorname{SNR}\left(\widetilde{y}_{\alpha\beta}\right)$. As $\operatorname{SNR}\left(\widetilde{r}_{ab}\right)\rightarrow\infty$, $\widetilde{\theta}_{ab}$ converges in distribution to a Gaussian random variable. On the contrary, as $\operatorname{SNR}\left(\widetilde{r}_{ab}\right)\rightarrow 0$, $\widetilde{\theta}_{ab}$ converges in distribution to a uniformly distributed random variable on $[0,2\pi)$. Because $\lim_{x\rightarrow 0}\frac{x}{\sin(x)}=1$, and $\sigma\left(\widetilde{\theta}_{ab}\right)$ conditioned on no wrapping is inversely proportional to $\operatorname{SNR}\left(\widetilde{r}_{ab}\right)$,

Similarly, we can obtain covariance given no wrapping

and $\operatorname{cov}\left(\widetilde{\theta}_{ab},\widetilde{\theta}_{cd}\right)=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_{\alpha\beta}$. To ameliorate estimation difficulty and use complex measurements $\widetilde{y}_{\alpha\beta}$ only, we further approximate and scale the variance and covariance:

where $\mathrel{\vbox{ \offinterlineskip\halign{\hfil$#$\cr\propto\cr\kern 2.0pt\cr\sim\cr\kern-2.0pt\cr}}}$ denotes “approximately proportional to.” Here, the scaling factors are the same for (22) and (25). Let

Then, with the elementwise approximations above, we can formulate the scaled approximated $\bm{\Sigma}\left(\widetilde{\bm{\theta}}\right)$ directly from observed voxel magnitudes. This scaling in (22) and (25) does not affect the estimator $\widehat{v}$, as seen from (36) below. Finally, because the true covariance matrix is close to rank deficient, the element-wise approximations in (22) and (25) can potentially violate the positive semi-definite property of a covariance matrix. Accordingly, we follow the approximation step by projection to the closest positive semi-definite matrix [14]. This projection operator, $\Pi(\cdot)$, for a symmetric matrix $\mathbf{M}$ with eigen-decomposition $\mathbf{M}=\mathbf{V}\mathbf{S}\mathbf{V}^{-1}$ is given by

where $\max(\mathbf{S},0)$ is applied element-wise to the eigenvalues.

To illustrate accuracy of covariance modeling, we adopt the cosine similarity metric, which is scale invariant. For modeled covariance $\Pi\left(\bm{\Sigma}\left(\widetilde{\bm{\theta}}\right)\right)$ and sample covariance $\widetilde{\bm{\Sigma}}\left(\bm{\theta}\right)$, the cosine similarity is

Results are computed for the case of a single-coil, symmetric encoding, and intra-voxel dephasing $2s_{11}=s_{21}=2s_{31}$; $10^{6}$ random draws are used at each $s_{21}$ to provide a sample covariance matrix and to compute the mean cosine similarity. Fig. 2 displays the similarity metric results, from which we see excellent agreement for $s_{21}>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 [7].

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, $\bm{n}$, in (12) conditioned on no wrapping is

Based on the covariance derivation in 2.2, we can now formulate the estimator. We adopt two approximations suitable when the SNR is not extremely low. Our first approximation is that $\bm{n}$ follows a joint Gaussian distribution with covariance matrix given in (30),

Denote the velocity estimate $\widehat{v}$ given wrapping integers $\bm{k}$ as $\widehat{v}_{\bm{k}}$. Then, due to (31), $\widehat{v}_{\bm{k}}$ and its resulting RMSE given no wrapping, $\sigma(\widehat{v}_{\bm{k}})$, are

where

and $\langle\bm{a}\rangle_{\bm{b}}$ denotes remainders after elementwise modulo $\bm{a}$ by $\bm{b}$. Thus, the estimate $\widehat{v}_{\bm{k}}$ is a weighted sum of unwrapped noisy velocities.

The second assumption we adopt is

where $\preccurlyeq$ is element-wise less than or equal. This approximation yields the likelihood $\mathbb{P}(\widetilde{\bm{v}}|v)$

where $\operatorname{d}_{\bm{z}}(\bm{x},\bm{y})$ is a “wrapped displacement” between $\bm{x}$ and $\bm{y}$ with respect to $\bm{z}$,

We use $-,\lfloor\cdot\rceil,\oslash$ to denote element-wise subtraction, rounding, and division. One could perform search over all possible $\bm{k}$ to minimize the negative log likelihood,

In 2.4, we instead present a fast method to detect the best wrapping integers, $\bm{k}^{\star}$. In 2.6 and 2.8, we explicitly consider three-point encoding, $N_{e}=3$, to provide concrete results and to optimize the design of venc and underlying first moments.

We introduce a fast estimator based on (32) and detection of the wrapping integers, $\bm{k}$. We refer to this estimator as Phase Recovery from Multiple Wrapped Measurements (PRoM). PRoM extends the prior signal processing result in [15] 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\in[0,\Omega)$, define the set $\mathcal{K}^{\prime}(\widetilde{\bm{v}})$ of wrapping integers

By (35), $\mathbb{P}(\bm{k}\in\mathcal{K}^{\prime}(\widetilde{\bm{v}}))\approx 1$. So, from (32) and (36), we can minimize the negative log likelihood

with probability $\approx 1$. Next, we leverage the equality constraint in (12) combined with the second approximation in (35) to decrease the cardinality of $\mathcal{K}^{\prime}(\widetilde{\bm{v}})$, denoted as $|\mathcal{K}^{\prime}(\widetilde{\bm{v}})|$. We have

which is equivalent to

where $\lceil\cdot\rceil$ is element-wise ceiling function. Then, the pruned search set for $\bm{k}$ can be expressed

So, the parsimonious construction considers only $v\in[0,\Omega)$ such that $-\tfrac{1}{2}+\tfrac{v-\widetilde{v}_{ab}}{2\text{venc}_{ab}}$ is integer for any $2\text{venc}_{ab}$. The cardinality of the search set

Thus, the number of searches is bounded by the summation of $\bm{h}$ instead of product of $\bm{h}$. Then, the minimization in (41) can be reduced to

with probability $\approx 1$. Together, the construction of the pruned set, $\mathcal{K}(\widetilde{\bm{v}})$, of candidate wrapping integers and the efficient search over $\{\widehat{v}_{\bm{k}}|\bm{k}\in\mathcal{K}(\widetilde{\bm{v}})\}$ to minimize $\mathcal{L}(\widetilde{\bm{v}},v)$ comprise PRoM, an approximate MLE of $v$. For a general $N_{e}$-point encoding, PRoM pseudo-code is given in Alg. 1, and PRoM code is available at https://github.com/Zhao-Shen/PRoM.

Fig. 3 illustrates $\mathcal{L}(\widetilde{\bm{v}},\widehat{v})$ for $\textbf{venc}=[35,10,14]^{\intercal}$ cm/s. The searched candidates $\{\widehat{v}_{k}|\bm{k}\in\mathcal{K}(\widetilde{\bm{v}})\}$ are marked by superimposed red dots. For this case, $|\mathcal{K}^{\prime}(\widetilde{\bm{v}})|=252$ and $|\mathcal{K}(\widetilde{\bm{v}})|=14$; thus, 94.4% of the search space $\mathcal{K}^{\prime}$ is bypassed via the proposed construction of $\mathcal{K}(\widetilde{\bm{v}})$. If $\bm{n}$ is known to concentrate in a smaller volume compared to the second assumption (35), or $v$ is known and restricted to a range less than $\Omega$, then $\mathcal{K}(\widetilde{\bm{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 spacing used in [7] entails computation for $7000$ candidate velocities. In contrast, PRoM only requires search over only $|\mathcal{K}(\widetilde{\bm{v}})|=14$ candidates, yielding a $500$-times reduction.

PRoM admits a simple geometric interpretation. Observe that the noisy velocity measurement $\widetilde{\bm{v}}$ resides in a hyper-rectangle $\{\widetilde{\bm{v}}|\bm{0}\preccurlyeq\widetilde{\bm{v}}\preccurlyeq 2\textbf{venc}\}$. The vector of noiseless velocity measurements $\langle v\rangle_{2\textbf{venc}}$ for $v\in[0,\Omega)$ is a point in the hyper-rectangle lying on wrapped line segments parallel to $\bm{1}$. Then, $\widehat{v}$ is found by an oblique projection of the noisy $\widetilde{\bm{v}}=\langle v+\bm{n}\rangle_{2\textbf{venc}}$ to the closest line segment. Here, the “oblique projection” is determined by the $\bm{\Sigma}^{-1}(\bm{n})$ weighted distance. The search for the closest line segment is reduced to search for $\bm{k}\in\mathcal{K}(\widetilde{\bm{v}})$.

In this section, we derive $\mathbb{P}(\widehat{v}|v)$, the distribution of the estimated velocity given the true velocity; this somewhat technical derivation, in turn, enables the optimized design of venc and the underlying first moments. To derive the distribution, we first establish two lemmas. The first lemma tells us that adding the same constant, $\eta$, to all noise realization components does not affect the error in detecting the wrapping integers and simply adds $\eta$ to the PRoM velocity estimate, modulo $\Omega$.

Fig. 4 provides visualisation of wrapping integers $\bm{k}^{\star}(\widetilde{\bm{v}})$ for the case $\textbf{venc}=[35,10,14]^{\intercal},[s_{11},s_{21},s_{31}]=[2.5,5,2.5]$. All $\widetilde{\bm{v}}$ along direction $\bm{1}$ inside the hyper-rectangle share the same $\bm{k}^{\star}$, which is a consequence of Lemma 1.

In addition, $\mathcal{L}(\widetilde{\bm{v}},v_{\bm{k}})$ as a function of $\widetilde{\bm{v}}$ is piece-wise quadratic with the same curvature for all $\bm{k}$, so the decision boundaries of $\bm{k}^{\star}(\widetilde{\bm{v}})$ are linear, which is also illustrated in Fig. 4.

By Lemma 1, the error in wrapping integers, $\langle\bm{k}^{\star}(\widetilde{\bm{v}})-\bm{k}(\widetilde{\bm{v}})\rangle_{\bm{h}}$, remains constant for all noise realizations $\bm{n}$ along any line parallel to $\bm{1}$. So, we can divide the space of all possible noise realizations, $\mathbb{R}^{\frac{N_{e}(N_{e}-1)}{2}}$, into “tubes” $\mathcal{T}(\bm{x})$ parallel to $\bm{1}$, based on the difference, $\bm{x}$, of estimated wrapping integers $\bm{k}^{\star}$ and true wrapping integers $\bm{k}$.

We can, as seen below, integrate over these tubes to arrive at error probabilities for detecting the wrapping integers. Next we establish a second lemma describing the orthogonality between pairwise noise differences and the error in the estimated velocity.

Armed with the two lemmas, we can specify the distribution.

Let $f(\widehat{v}|\mathcal{T}(\bm{x}),v)$ denote the conditional Gaussian probability density function (pdf) in Theorem 1 without wrapping by modulo $\Omega$. Thus, by wrapping the pdf function and invoking the law of total probability, we have

To complete (51), we need only to calculate $\mathbb{P}(\bm{n}\in\mathcal{T}(\bm{x})|v)$ 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}(N_{e}-1)}{2}-1$ variables.

Fig. 5 illustrates the histogram of $\widehat{v}|v=0$ using $10^{5}$ trials, $N_{e}=3$, $N_{c}=1$, $\textbf{venc}=[99,18,22]^{\intercal}$ cm/s at $s_{21}=5,10$. Invoking Theorem 1, we predict the five most probable wrapping integers, which correspond to detection regions $\mathcal{T}([0,0,0]^{\intercal}),\mathcal{T}([5,4,1]^{\intercal}),\mathcal{T}([6,5,1]^{\intercal}),\mathcal{T}([1,1,0]^{\intercal})$, and $\mathcal{T}([10,8,2]^{\intercal})$. These five predicted detections correspond to $f(\widehat{v}|\mathcal{T}(\bm{x}),v)$ centered at $0$ (the true velocity), $\pm 177.81$, and $\pm 40.37$ cm/s; these five predictions are validated in the histogram. For the higher SNR case of $s_{21}=10$, the probability of unwrapping errors goes very small, and the $10^{5}$ trials are insufficient to encounter an unwrapping error to $\pm 40.37$ cm/s.

We consider here three-point encoding for velocity in one direction. Due to (9, 14), every $\text{venc}_{ab}$, and hence the unambiguous range, $\Omega$, depends only on the differences between first moments; thus, $\text{venc}_{ab}$ and $\Omega$ are unaffected by adding the same constant to every first moment. We assume the following ordering for the first moments:

Thus, the three vencs are ordered: $\text{venc}_{31}<\text{venc}_{32}<\text{venc}_{21}$. Let $\xi=\frac{\text{venc}_{32}}{\text{venc}_{31}}$, noting that (9) and (52) imply $\xi\in(1,2)$. For rational $\xi=p/q$ with co-prime integers $p$ and $q$, the unambiguous range $\Omega$ in (14) is

Thus, by jointly unwrapping multiple vencs, one can construct an unaliased velocity range that is larger than the highest venc, $\text{venc}_{21}$, by a factor of $2(p-q)$. The covariance matrix for three-point encoding can be formulated via (30). Correspondingly, we can calculate the combination weights $\bm{w}$ and the resulting RMSE given wrapping integers (32, 33). Armed with the explicit error variance and the probability of unwrapping errors derived below, we present in 2.8 an optimized design of venc 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.

In the notation of (10) and (13), the approaches in [3], [4] use the unaliased high venc measurement, $\widetilde{v}_{21}$, to unwrap the low venc measurement, $\widetilde{v}_{31}$, while $\text{venc}_{32}$ goes unused. The estimator in [4], which we denote as standard dual-venc (SDV), is given by

In [7], two potentially aliased measurements $\widetilde{v}_{31},\widetilde{v}_{32}$ are jointly unwrapped by minimizing