Improved MR phase‐contrast velocimetry using a novel nine‐point balanced motion‐encoding scheme with increased robustness to eddy current effects

Phase‐contrast MRI (PC‐MRI) velocimetry is a noninvasive, high‐resolution motion assessment tool. However, high motion sensitivity requires strong motion‐encoding magnetic gradients, making phase‐contrast‐MRI prone to baseline shift artifacts due to the generation of eddy currents. In this study, we propose a novel nine‐point balanced velocity‐encoding strategy, designed to be more accurate in the presence of strong and rapidly changing gradients. The proposed method was validated using a rotating phantom, and its robustness and precision were explored and compared with established approaches through computer simulations and in vivo experiments. Computer simulations yielded a 39–57% improvement in velocity–noise ratio (corresponding to a 27–33% reduction in measurement error), depending on which method was used for comparison. Moreover, in vivo experiments confirmed this by demonstrating a 26–53% reduction in accumulated velocity error over the R–R interval. The nine‐point balanced phase‐contrast‐MRI‐encoding strategy is likely useful for settings where high spatial and temporal resolution and/or high motion sensitivity is required, such as in high‐resolution rodent myocardial tissue phase mapping. Magn Reson Med, 2013. © 2012 Wiley Periodicals, Inc.


INTRODUCTION
Phase-contrast MRI (PC-MRI) velocimetry is a well-established technique allowing noninvasive measurement of three-directional motion in all four dimensions (3Dþtime), with a spatial resolution identical to the acquisition matrix and a temporal resolution approach-ing that of echocardiography (1). PC-MRI is widely used in blood flow analysis (2)(3)(4) and is emerging as a potent and versatile tool in myocardial motion assessment, both in humans (5)(6)(7)(8)(9) and in rodents (10)(11)(12)(13). When used to assess tissue motion, PC-MRI is often referred to as tissue phase mapping (TPM).
In conventional MR imaging, the phase component of the signal is often discarded and only the magnitude images are used. However, in PC-MRI, motion-sensitizing gradients systematically encode the spin velocity into the signal phase, and the motion information can be reconstructed from the complex images. The range of velocities that can be measured without aliasing (the venc) must be set a priori and is inversely proportional to both the amplitude and the square of the duration of the motion-encoding gradients. Relative to blood flow that can reach velocities in the range of meter per second, myocardial velocities are low (in humans, <8 cm/s (1); in mice, <3 cm/s (10)). Consequently, cardiac TPM requires stronger (and/or longer) gradients than blood flow PC-MRI. Moreover, high-resolution (both temporal and spatial) imaging also requires strong gradients and short ramp times. Together, this can result in the generation of strong eddy currents in the tissue and any surrounding conductive material. These currents will subsequently create additional magnetic fields that perturbate the carefully controlled magnetic environment and introduce errors in the MR signal (14). In turn it can lead to significant errors in the velocity measurements (15). Hence, measuring low velocities, or with high spatial and temporal resolution, requires methods compensating for errors due to eddy current-induced baseline phase shifts. Unfortunately, eddy current-induced artifacts are not straightforwardly handled analytically (16) and must to a large extent be dealt with off-line relying on numerical approaches (17,18).
In this study, we develop and evaluate a novel PC-MRI technique, designed with the objectives to (1) increase robustness to baseline phase shifts, thus allowing stronger gradient amplitudes and/or shorter gradient durations and (2) provide a minimum of scan time penalty, both compared with established encoding schemes. Together, this makes high-resolution imaging with optimized velocity sensitivity and accuracy feasible, allowing detailed studies on, for example, the myocardial motion pattern in rodents.

MR Signal and PC Imaging
The MR signal from a group of spins with mean velocitỹ v 0 is (19) ðtÞtdt is the first moment of the linear gradient G(t) and U is the baseline phase from all other contributors, including higher order-of-motion components. c is the gyromagnetic ratio.
In PC-MRI,ṽ 0 is determined from Eq. [1] by the acquisition of several datasets with different first moments, but with all other imaging parameters kept unchanged (20). U, which is considered independent ofM 1 , is effectively eliminated by complex conjugate multiplication of two such datasets: where DM 1 ¼M 1;A ÀM 1;B . From the exponential of Eq. [2], the acquisition process can be written compactly (19) u where A is the acquisition matrix containing the values of gDM 1 used. A is of size N Â K, where N is the number of velocity-encoding steps and K is the number of encoding directions. The velocity vector is readily reconstructed from where A þ is the Moore-Penrose pseudoinverse of A (19,21). As phase information is only uniquely defined in the interval (Àp, p), aliasing will occur if the acquisition generates data outside this range. The highest velocity that can be measured without aliasing to occur (the venc) may be expressed directly in terms of the acquisition matrix (19): where k is the encoding direction and min is the minimum of all encoding sets (i.e., the rows of A).

Eddy Current Effects in PC-MRI
One challenge in PC-MRI is the creation of eddy currents in the system. Eddy currents are, from Faraday's law, created in conducting materials when subjected to a time-varying magnetic flux and will cause additional magnetic fields opposing the original inducing field (in accordance to Lenz's law). The eddy currents build when the magnetic field changes, i.e., during gradient ramps, with a rate proportional to the slew rate of the gradient (22). Thus, the use of stronger gradients creates stronger eddy currents.
The magnetic fields associated with the eddy currents enter as a part of F in Eq. [1]. However, the phase difference calculation [Eq. 2] does not cancel eddy currentinduced shifts, as they are not independent ofM 1 and thus additive in the reconstruction process (17,22). The contaminating phase contribution from eddy currents in a single voxel in the nth acquisition can be expressed as where (A) n is the nth row in the acquisition matrix. The baseline phase shift u n,ec is a function of several parameters, including the directions and amplitudes of the motion-encoding gradients. The presence of baseline phase shifts in PC-MRI data and their dependence on gradient configuration are seen in Fig. 1a-c. In addition to eddy currents, concomitant (Maxwell) gradient effects are in general not eliminated in the phase difference reconstruction. However, in contrast to eddy current errors, concomitant errors are readily corrected analytically (23).

Eddy Current Compensation
Several approaches for baseline phase shift correction in PC-MRI exists, hereby known as eddy current compensation (ECC). A set of techniques uses the phase data from a separate scan (in general of a stationary phantom) to estimate the distribution of baseline shifts (24). These approaches will, however, cause an increase in experimental time.
ECC may also be performed retrospectively without the need for additional scans. A commonly used strategy is a global subtraction of the phase from a static structure nearby the region of interest. This, however, will not capture the spatial variation in the baseline phase shift. An alternative method, and also our method of choice, is introduced by Walker et al (14). Here, the temporal standard deviation for the phase in each image pixel is calculated (this requires temporal data, e.g., a cine loop), and a mask is created whose pixels have standard deviation beneath a user-defined threshold value. These pixels are assumed to represent static regions in the field of view. In each cine set, a surface is fitted to the masked phase data, including all time points in the fit. The result is subtracted from the individual cine frames. The process is repeated for all N encoding steps.
In general, ECC modifies Eq. [6] to where e n ¼ u n,ec À u n,ECC is the residual phase error (fitting residual) from the nth acquisition. After least-squares fitting, the fitting residuals over all pixels in one acquisition can be assumed to be distributed with mean zero (25).  In the images shown, the motion-compensated reference scan has been subtracted. Furthermore, the ''best fit'' from the ECC protocol (which uses a surface fit of a 2D polynomial of degree 3) is shown (d-f), along with post-ECC phase data (g-i) where the fitted polynomials have been subtracted from the pre-ECC datasets. Velocity reconstruction both before and after ECC was performed following Eq. [4], and the spatial distribution of velocity magnitude is shown in both cases for all three encoding strategies (j-l). As no net motion is expected in the phantom, any measured non-zero velocity can consequently be treated as measurement error. The RMS and standard deviation of the reconstructed velocity magnitude is shown (m,n). A significant reduction in RMS error is seen in the nine-point method compared with the four-and five-point method (all data points included in analysis; P < 0.001 in all cases, N ¼ 10,313). Presumably, the spread in the error bars in m-n is largely nonrandom, due to the spatial variation in residual phase error across the phantoms (seen in j-l where N' is the number of non-zero elements in the kth column of A (i.e., the number of acquisitions with nonzero velocity encoding in direction k). When the fitting residuals are distributed around zero, equally likely to be positive or negative, the signs of the elements in (A þ ) k can be discarded. Thus, Eq. [9] essentially averages the residuals, and v k,err will be reduced following the concept of signal averaging.

PC-MRI-Encoding Strategies
Several different encoding schemes (i.e, different architecture of A) for three-directional PC-MRI have been described, among them referenced and balanced fourpoint (20,26) and balanced five-point (19). The referenced four-point method ( Table 1, first column and Fig. 2b,e) uses three motion-encoded scans (N 0 /N ¼ 1/3) withM 1 applied in each of the three Cartesian directions, in addition to a completely motion-compensated reference scan to eliminate U as in Eq. [2]. The balanced four-point method applies motion-encoding gradients in all four scans (N 0 /N ¼ 4/4), toggling the moments in pairs. The latter approach reduces the variance in the velocity data (20) but has a disadvantage in a rather complex spatial variance in dynamic range (19). The balanced five-point encoding scheme (Table 1, second column) is based on the balanced four-point method but adds one motion-compensated acquisition as a reference scan. Thereby, it achieves as high velocity-to-noise ratio (VNR) as balanced four-point method but manages to keep the aliasing space more symmetrical and thus efficient (Fig. 2c,f). The advantages and limitations of these different methods are excellently analyzed by Johnson and Markl (19). We propose an encoding approach extending the balanced five-point encoding to a balanced nine-point encoding ( Table 1, third column), where all eight permu-tations of the three-directionalM 1 are used in separate acquisitions (N 0 /N ¼ 8/8) in addition to a completely motion-compensated scan for use as reference in Eq. [2]. We propose that this scheme improves the accuracy in the velocity measurements because the number of different encoding steps N ¼ N 0 is larger, and thus more unique baseline phase shift maps are estimated using the ECC algorithm. Ultimately, this is expected to reduce the error in the final velocity reconstruction (Eqs. [8 and 9]). The dynamic range of the balanced nine-point method is identical to that of the balanced four-point ( Fig. 2d,g).
The three methods compared in this study are the referenced four-point, balanced five-point, and balanced nine-point methods. For simplicity, they will hereafter be referred to simply as the four-, five-, and nine-point method, respectively.

METHODS
Aiming to verify and investigate the performance of the nine-point method, computer simulation as well as phantom and in vivo TPM experiments were performed intending to allow comparison between the four-, five-, and nine-point method.

MRI Experiments Set-Up
All MRI experiments were performed on a 9.4 T/210 mm/ASR horizontal bore magnet (Agilent Technologies, Inc., Santa Clara, CA) with a high-performance actively shielded gradient coil (inner diameter 60 mm, rise time 180 ms, max strength 1000 mT/m). Motion-encoding gradients were incorporated into an radiofrequency-spoiled, three-directional motion-compensated gradient echo sequence (10). The acquisition matrices for the three encoding strategies are found in Table 1.

Static Phantom MRI Experiments
To investigate the effect of the gradient configuration on time-independent baseline phase shifts, a static agar phantom was imaged three times using PC-MRI, using the four-, five-, and nine-point encoding strategies. ECC was performed according to the procedure of Walker et al. by Table 1 Fundamentals of the Three Encoding Strategies Compared in the Computer Simulation and the TPM Experiment Venc: upper threshold for velocity that can be measured without aliasing. a The gradient moment used in the four-point method was twice as large to obtain roughly the same venc as in five-and nine-point.
Improved PC-MRI Using Nine-Point Motion Encoding fitting a two-dimensional (2D) polynomial of degree 3 to the phase data after subtraction of the motion-compensated reference dataset. Velocity reconstruction was performed (Eq. [4]) both without and with ECC.
The following imaging parameters were applied: echo time/pulse repetition time ¼ 2.60/3.54 ms; field of view ¼ 30.0 mm Â 30.0 mm; acquisition matrix 128 Â 128; slice thickness ¼ 1.00 mm; flip angle 25 ; receiver bandwidth ¼ 156.25 kHz; venc ¼ 5.0 cm/s (four-point) and 5.8 cm/s (five-/nine-point); NA ¼ 45/36/20 (four-/five-/ nine-point). The reason for choosing a rather high NA was to achieve exactly the same experiment time (i.e., signal-to-noise ratio (SNR)) for each acquisition. The datasets were not time resolved, meaning that they represent the time-independent baseline phase shifts.

Computer Simulation
For the computer simulation, a program was written in Matlab (The MathWorks, Natick, MA), where a rotating cylinder was modeled and a PC-MRI acquisition process was simulated using Eq. [3]. The operator chose all parameters including venc, encoding strategy, and signal averaging (NA). Subsequently, computer-generated baseline phase shift maps (Fig. 3a) and random pixel-wise noise were added to the simulated phase datasets individually. Then, the reference scan was subtracted and ECC was simulated according to the method of Walker et al., as a 2D polynomial of degree 3 was chosen for fitting a random selection of 10% of the stationary points (i.e., points outside the cylinder). An example of a typical computer-constructed baseline phase shift map, the ECC polynomial best fit and the resulting residual phase errors are shown in Fig. 3. Finally, the velocity components were reconstructed according to Eq. [4].
Several consecutive simulations were run applying the three different encoding algorithms (four-point, fivepoint, and nine-point) with NA ¼ 3/2/1, respectively, but all other parameters unchanged (including baseline phase shift amplitude). The relative acquisition times for and y (a). The maximum vector length jṽj that can be measured without aliasing to occur is calculated from Eq. [5] and displayed in grayscale as a function of vector orientation (b-d). The resulting venc corresponds in each figure to the lowest (darkest) value in each map. The map variance may be regarded as a measure of the spatial efficiency of the encoding. (e-g) Volume renderings of the dynamic ranges of the three methods intuitively illustrating the difference between the venc (the rendered spheres) and the actual directional dependency of the aliasing threshold (19). In the four-point referenced encoding strategy (b, e), aliasing occurs first if the velocity is directed parallel to any Cartesian direction. The five-and nine-point balanced-encoding schemes (c and d, f and g) share the same dynamic range, with aliasing occurring first at |v the simulated experiments were 1.33, 1.11, and 1.00 for the four-, five-, and nine-point method, respectively. Any difference in reconstruction accuracy between methods was investigated using both the error data and a VNR estimate from where v max is the highest velocity in the original data and s is the root-mean-square (RMS) of the pixel-wise deviation between the original data and the reconstructed data, taken over the whole disk.
To visualize the results of the computer simulations, two different approaches were used. First, the velocity magnitude of the data was plotted, illustrating change in the absolute value of the reconstructed velocity. Second, streamline plots were used to illustrate the change in direction of the reconstructed velocity, a phenomenon occurring due to the decomposition of the three-dimensional (3D) motion during the encoding procedure under different gradient configurations.

Rotating Phantom MRI Experiments
To evaluate the accuracy of the motion quantification methods, a rotating phantom was imaged using PC-MRI using the four-, five-, and nine-point approaches. An agar gel cylinder (Ø ¼ 2.68 cm) was placed in the isocenter of the magnet and was rotated by a motor (2.02 s/rot), resulting in a maximum speed of 4.17 cm/s at the edge of the phantom. Contrary to the computer simulation and in vivo experiments, no section in the images (except the center point of the phantom) was stationary and thus suitable for use in the ECC algorithm. Therefore, two separate datasets with identical parameters were acquired, one with motion and one without (i.e., the rod from the motor to the phantom was disconnected, while the motor was kept running). The latter was used in the ECC algorithm.
The imaging parameters were identical to the static phantom experiments parameters. The resulting velocity measurements were analyzed in Matlab and compared with the expected distribution of velocities in the phantom.

In Vivo MRI Experiments: Data Acquisition
In the in vivo MRI experiment, midventricular short-axis TPM cardiac cine loops from three healthy mice were recorded. Slice planning was performed as described by Schneider et al (27). All acquisitions were repeated three times using the four-point, five-point, and the nine-point method, respectively, but with all other parameters (except NA) unchanged, resulting in a total of nine cine loops. The acquisition order was varied in the three hearts. Although the acquisition was prospectively triggered by both electrocardiography (ECG) and respiration, steady-state magnetization was maintained by continuous radiofrequency and gradient activity also during respiration, while no data were acquired (28). The acquisition was noninterleaved, i.e., each encoding step was separately recorded with full temporal resolution (29,30).
All animals were cared for according to the Norwegian Animal Welfare Act, and the use of animals was approved by the Norwegian Animal Research Committee, conforming with the Guide for the Care and Use of Laboratory Animals published by the US National Institutes of Health (NIH Publication No. 85-23, revised 1996). Anesthesia was induced in an anesthesia chamber with a mixture of O 2 and 4% isoflurane, and further maintained by administration of a mixture of O 2 and 1.5-2.0% isofluran in freely breathing animals. Heart rate, body temperature, and respiration rate were constantly monitored during experiments. The level of anesthesia was regulated to keep the heart rate as close to 450 bpm as practically possible, and the body temperature was maintained by heated air.
The imaging parameters were as follows: echo time/ pulse repetition time ¼ 2.44/3.51 ms; field of view ¼ 25.6 mm Â 25.6 mm; acquisition matrix 128 Â 128; slice FIG. 3. The fitting error in the ECC procedure. a: An example of a computer-simulated nonpolynomial baseline phase shift map. These distributions were created by generating two low-resolution noise maps that were interpolated in 2D to high-resolution and summed, creating a random but smooth map of required resolution. The purpose of this was to create phase shift maps, unique for each encoding step that were spatially well behaved without discontinuities but not straightforwardly described by any polynomial. The outline of the simulated cylinder is shown. b: The polynomial best fit (using a polynomial of degree 3) of the actual map, calculated on the basis of a random selection of about 10% of the v ¼ 0 data points (outside the cylinder). This simulated the in vivo situation where stationary tissue nearby the region of interest has to be used for fitting the ECC map. c: The fitting error, i.e., the difference between the actual map (a) and the best fit (b), illustrating the inexactness in the ECC procedure and thus the source of eddy-current-induced errors in the reconstructed velocity (Eq. [8]).

In Vivo MRI Experiments: Postprocessing
The resulting datasets were semiautomatically segmented and analyzed in a purpose-written Matlab software. Required human inputs during postprocessing were definitions of (1) key time points and the corresponding left ventricular (LV) endocardial and epicardial wall, (2) the position of the papillary muscles at time t ¼ 0, and (3) regions in the image frame subject to fold-over artifacts (see Discussion section). Whenever human input was required, the operator was blinded to what encoding scheme was used. Regardless, potential human bias was evaluated by running a separate, supplementing postprocessing process, with one common segmentation mask for all datasets.
ECC was performed following the procedure of Walker et al. The 4% pixels with the lowest temporal standard deviation were used for fitting a 2D polynomial of degree 3.
The myocardium data points were furthermore automatically divided into 64 sectors (Fig. 4c). No data point was included in more than one sector. The whole segmentation was temporal dynamic, following both radial motion and in-plane rotation of the myocardium. The velocity vector in each pixel was decomposed into a cardiopolar coordinate system, constituting of in-plane radial and tangential components and a through-plane longitudinal component (13). The in-plane velocities were corrected for bulk cardiac movement by subtraction of the global (mean) in-plane motion (31). No signal filtering was used in postprocessing. To avoid temporal jitter, the individual datasets were normalized to end-systolic time defined by the first minimum peak in global radial motion (1).
The 3D motion path of each sector was subsequently calculated by integrating the 3D velocity data:x When no errors were present in the velocity data, these paths are expected to be closedxðt max Þ ¼ 0 ð Þover the course of an R-R interval. Conversely, any inaccuracies present in the motion data would through the integration accumulate to an error in the path. Therefore, we chose to usexðt max Þ, labeled residual error, as a parameter for determining the error in the velocity data. Velocity integration was performed and the mean 3D path and residual error for the nine datasets were calculated. In addition, the dependence of the error on the mean velocity in each sector was investigated through correlation analysis.
SNR was estimated from the magnitude images as the ratio of maximum signal intensity to the standard deviation of the signal in a region outside the object (i.e., in air).

Data Analysis and Statistics
Data postprocessing was performed in Matlab. Statistical analysis was performed in Matlab and Minitab (Minitab Inc., State College, PA). Two-tailed paired t tests and one-way unstacked analysis of variance (ANOVA) were used to determine statistical significance; the threshold for statistical significance was set at P < 0.05. All data presented are mean 6 SD.
In the computer simulation data, five consecutive runs were performed and results compared component-wise (total N ¼ 15). In the in vivo MRI experiments, three mice hearts were used and the slices were divided into 64 sectors. Thus, intra-animal comparison had N ¼ 64 and inter-animal comparison had N ¼ 192. Paired analysis was always performed animal and sector wise.

Static Phantom MRI Experiments
The PC-MRI images of the static phantom, illustrated in Fig. 1, clearly illustrates the presence of baseline phase shifts and their dependence on gradient configuration The location in the septal wall whose angular distance (with the LV in-plane center of mass as origin) to the papillary muscles are equal (a), defines the border between the first and the last sector (initial angle). b: The sectors are positioned counting clockwise in this view (observed from the base, with LV to the right), each covering an angle 2p/S where S is the number of sectors. c: The sector-wise segmentation in this study, where 64 sectors have been used. Each sector contains 1-8 pixels, depending on wall thickness and position in time. No data points are included in more than one sector. The sectors are dynamic in time and follows the rotational motion of the LV slice throughout the cardiac cycle. The rotation is calculated from integration of the tangential velocity.
( Fig. 1a-c) as well as the ECC surface fits (Fig. 1d-f) and post-ECC fitting residuals (Fig. 1g-i). Furthermore, as no net motion was expected in the phantom, the reconstructed velocities allowed comparison of measurement error between methods (Fig. 1m,n). The nine-point method exhibits a significant reduction in RMS error compared to the four-and five-point methods.

Computer Simulation
The computer simulation allowed for the systematical comparison between the original data and the data processed and reconstructed by use of the four-, five-, and nine-point methods. In Fig. 5, an example of a qualitative result of a single simulation is shown, where both error in magnitude and direction are illustrated. In the velocity magnitude images (Fig. 5a-d), the shape is noticeably more symmetrical and close to the original data in the nine-point data than with the two other encoding strategies. The actual error in the magnitude reconstruction is shown in Fig. 5e-g. The velocity direction data (Fig. 5i-k) also exhibit less geometrical distor-tion in the nine-point case compared with the four-and five-point case. Ultimately, both velocity magnitude and direction are clearly less distorted in the case with ninepoint encoding than in the four-and five-point case.
The mean deviations between original and reconstructed data, listed in Table 2, show a significant 33% and 27% (both P < 0.05) reduction in RMS error and a similarly significant 57% and 39% (both P < 0.05) increase in VNR in the nine-point method compared with the four-point and the five-point, respectively.

Rotating Phantom MRI Experiments
The accuracy of the encoding approaches was evaluated using a rotating phantom. The results of the experiment are shown in Fig. 6. Data from all the three encoding strategies demonstrate excellent correlation with the expected motion, with a linear regression coefficient close to 1.

In Vivo MRI Experiments
In the raw datasets, the in-plane myocardium was described by 253.5 6 25.6 data points depending on cycle The difference between original (a) and reconstructed (b-d) velocity magnitude, i.e., the error after reconstruction. These plots also suggests a less corrupted reconstruction using the nine-point approach, both with respect to spatial distribution and general amount. Note the difference in grayscale axis between the velocity images (a-d) and the error images (e-g). h-k: The error in velocity direction illustrated using streamline plots, with is the ideal case (h) and the four-(i), five-(j), and nine-point (k) results. There is visually less directional distortion in the nine-point plot (k) than in the other two (i, j). These images suggest that the nine-point method provides less degradation of velocity information, both with respect to magnitude and direction. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] phase. Accordingly, each of the 64 sectors (illustrated in Fig. 4c) contained an average four unique myocardial pixels throughout the cardiac cycle. Figure 7 shows the mean motion of the hearts (Fig.7ad) and resulting 3D motion paths (Fig.7e-g) calculated from Eq. [11], averaged over all sectors and animals. No significant difference in peak motion, positive or negative, was found in any direction between methods. Furthermore, the motion paths demonstrate the mean global radial, tangential, and longitudinal motion for the entire slice. The nine-point motion path in Fig. 7g seems less distorted and exhibits a markedly smaller residual error compared with the five-and four-point paths (Fig. 7e,f).
The result of the quantitative analysis is illustrated in Fig. 8. The residual error data from the individual hearts are shown as box plots in Fig. 8a-c, whereas the mean error is found in Fig. 8g. All datasets demonstrate individual significant reduction in residual error (P < 0.001 for all sets) in the nine-point datasets compared with the two others (Fig. 8a-c). On average, a 53% and 26% (both P < 0.001) reduction in mean error is observed in the nine-point encoded data compared with the four-and the five-point, respectively (Fig. 8g).
The SNR performance of the different methods is expected to be dependent on the number of individual datasets acquired. Indeed, the SNR was measured to be significantly smaller in the nine-point datasets (56.55 6 6.53) than in the four-point set (66.98 6 7.64) and five-point set (60.65 6 6.53; both P < 0.001 vs. the nine-point set).
Additionally, correlation analysis was performed to investigate the effect of mean velocity on residual error reduction (i.e., the difference in residual error between the methods; Fig. 8d-f). Significant positive correlation was found in the four-versus nine-point comparison.
In a supplementary postprocessing run to evaluate potential human bias, one common segmentation mask was used for all three encoding sets in each heart. This resulted in a somewhat smaller but still strongly significant reduction in residual error, from 1.51 6 0.76 mm (four-point) and 1.03 6 0.48 mm (five-point) to 0.72 6 0.50 mm (nine-point; all P < 0.001), corresponding to a 52% and 30% reduction, respectively.

DISCUSSION
The presented data suggest that the proposed nine-point method offers a substantial improvement in motion assessment accuracy compared with existing approaches. We have presented results from computer simulations as well as phantom and in vivo experiments, all demon-strating significant decrease in error and improvement in accuracy compared to established methods.

Cardiac Cycle Coverage and Effect of Discrete Sampling
In the in vivo experiments, the residual distance after a 3D velocity integration was used as a measure on errors in the velocity data and served as a parameter for comparison between methods. However, the accuracy of cyclic trajectories from velocity data can be substantially improved by a method described by Pelc et al (32) using both forward and backward temporal integration. Nevertheless, as we intended to examine the presence of errors in the velocity data before any such correction, this algorithm was not implemented in this study.
The length of the cine loop was 130 ms. Thus, at 450 bpm, 97.7% of the cardiac cycle was covered. However, the error in the velocity data was estimated expecting that the motion path from velocity integration should be closed, and any missing part of the cycle would make this assumption incorrect. Nevertheless, as only a few milliseconds of data are lost on average, even hypothetical large velocities in the final stages of the diastole would not be expected to alter the present results significantly. Besides, all experiments share this deficiency in cycle coverage, and should thus not compromise the basis for comparison.
Discrete sampling of the velocity inevitably results in low-pass filtering (29,30). Thus, the peak values of the velocities are always to a certain degree underestimated. This may have a considerable impact on the path calculation, as it will result in an underestimate of larger velocities. Correlation analysis of the residual error versus maximum velocity did indeed yield a significant positive correlation in all datasets (P < 0.001 for all nine sets, data not shown). However, this is also an effect that affects all experiments equally and is not expected to alter the comparison results significantly.

Segmentation
In contrast to LV segmentation used for evaluation of parameters such as ejection fraction and myocardial mass (33), in the present context the purpose was solely to exclude all pixels not containing myocardium. Thus, areas with signal loss artifacts (such as in the myocardium-lung border) were not included in the myocardium mask.
The segmentation masks were dynamic, designed to follow the contraction and rotation of the myocardium. However, the segmentation did not account for throughplane motion, so in reality slightly different parts of the myocardium were observed at different times. Specifically, the midventricular slices were visited (in late systole and early diastole) by faster moving basal myocardium (13). The effect of this, if any, should be a bias toward higher longitudinal velocities. However, as the basal myocardium travels back during diastole, the overall effect of this was not expected to affect the error estimation from the velocity integration. . This occurred despite that the maximum motion of the phantom (4.17 cm/s at the rim) was lower than the venc (5.00 cm/s), due to a combination of large baseline shift in the phase data and the dynamic range of the four-point strategy (Fig. 2c). Line profiles through the phantom were calculated for each experiment (e-g), representing the mean motion in 61 equally spaced points along the lines in Fig. 6b The different methods exhibit no significant difference in peak velocity, positive or negative, in either velocity component. All velocity tracings (a-d) have been adjusted for variation in heart rate. e-g: Plots showing the mean 3D motion paths, over all sectors and animals, calculated by integration of the 3D velocity. The path starting point is marked with an asterisk. Note that the coordinate system in the path plots is not Cartesian, but the cardiopolar coordinate system with radial (the velocity component toward the LV center of mass), the tangential (clockwise when seen from the base), and longitudinal (toward the base). 2D projections of the motion are shown for clarity. The difference in path between the three encoding strategies are clearly seen in the plots, and the residual error (the red line) is significantly smaller in the nine-point. As expected, all paths show an initial positive radial motion (systolic contraction) and a subsequent negative radial motion (diastolic relaxation). Additionally, a bulk negative longitudinal motion (toward the apex) is initially seen in all datasets, followed by a positive trend, consistent with published data (13). However, the diastolic longitudinal motion seems to be grossly overestimated in the four-point dataset (e), contributing substantially to the residual error. Moreover, the subtle tangential motion (rotation) is visibly heavily corrupted in both the four-(e) and five-point (f) datasets. As the tracking of the LV sectors (see Fig. 4) depends directly on the tangential velocities, such errors will potentially affect the accuracy of the LV segmentation substantially. In this case, the motion paths (e-g) are calculated form the mean velocities over the whole LV, and are thus unaffected by segmentation.

ECC Algorithm
The ECC algorithm uses stationary tissue for estimation of the spatial distribution of baseline phase shifts (14). In our short-axis images, the heart was always placed close to the middle of the image frame, often resulting in foldover artifacts in the phase-encoding direction from remote tissue (e.g., as seen in the upper right parts of Fig. 4). Care was taken during slice planning so that the aliased tissue signal did not overlap the myocardium. However, as the location of the folded-over tissue in the image frame did not correspond to the true origin of the signal, it had to be excluded from use in the ECC algorithm. This was achieved by adding an early step in the postprocessing procedure where an exclusion mask was created, defining the misplaced tissue. Another issue with the ECC was that we chose to use a single, mean surface fit for each cine set, thus estimating the temporal mean of the baseline phase shift throughout the cardiac cycle. In cine PC-MRI, the baseline phase shift can in general be considered to be time independent (32). However, the small gap in gradient activity between two consecutive cardiac cycles may result in a transient variation in eddy currents. Nevertheless, we strived to keep this gap in excitation as short as possible (3.76 ms on average), to be able to get the best possible estimate for the temporal mean of the phase shift. A slight modification in the ECC procedure, in which each time frame would be treated individually (14), would capture these changes but at the cost of substantially decreased fitting accuracy due to lower number of data points in each fit.

Limitations
In the computer simulations, random baseline phase shift maps were created to simulate eddy currentinduced phase distortions in a PC-MRI experiment. Obviously, this does not correspond to the actual situation, where these maps are nonrandom functions of the gradient waveforms. However, we believe that the random maps fulfill their purpose in (1) being unique for each velocity-encoding step and (2) producing fitting residuals after ECC. Both these effects are clearly seen in the actual case in Fig. 1a-c, g-i.
The proposed method has an obvious drawback in needing nine separate acquisitions. However, as signal averaging is nearly always used in MRI, we demonstrate through computer simulations and in vivo experiments that applying the nine-point scheme with NA ¼ n and comparing the five-point scheme with NA ¼ 2 Â n, and the four-point scheme with NA ¼ 3 Â n results in a decrease in scan time (and thus SNR) but still improves measurement accuracy. In rodent cardiac TPM, NA > 2 is often used (10,11,34), making the nine-point method beneficial with respect to acquisition time. However, applying the four-point method with NA ¼ 2 (13) results in a shorter scan time than is possible with the ninepoint method, but at the cost of lower SNR. When compared with encoding strategies with considerably shorter acquisition time (such as the four-or five-point method with NA ¼ 1), the SNR in the nine-point method will be higher, but factors such as drift in cardiac rate could impair the measurements, primarily in diastole (35,36).
The argumentation for the expected benefit of the nine-point method is based on an assumption that the ECC algorithm is unable to create ''perfect'' fits of the true baseline phase shift distributions (i.e., the post-ECC phase images will contain non-zero fitting residuals). When, however, the ECC was indeed able to produce zero-residual fits, all values for e in Eq. [8] would be zero, and the nine-point method would offer no improvement over the five-point method. Nevertheless, Fig. 1g-i shows that there indeed were non-zero post-ECC fitting residuals, and that they varied depending on gradient configuration. This was supported by the results of both the moving phantom experiment (Fig. 6) and the in vivo experiments.

CONCLUSIONS
We have presented a novel nine-point balanced motionencoding strategy for use in PC-MRI. The encoding scheme is shown from computer simulations as well as phantom and in vivo experiments to provide higher accuracy in velocity reconstruction than existing methods. We argue that this is due to increased robustness to baseline phase shifts in the MR data, which is due to eddy current effects. The proposed method offers better performance when strong gradients are necessary, such as when high sensitivity and/or high resolution is required. The method has no penalty in acquisition time when compared to methods with comparable SNR, where signal averaging NA > 2 is used.

ACKNOWLEDGMENTS
Many thanks to Øyvind Nordbø at the Norwegian University of Life Sciences, Å s, Norway for insightful comments on the postprocessing software and Roy Trondsen and William Louch for invaluable help with designing/ constructing the MR phantoms.

APPENDIX
In all PC-MRI-encoding schemes treated in this article, the acquisition matrices A have (1) the same absolute value (gDM 1 ) in every non-zero element and (2) mutually orthogonal columns with equal norm h ¼ ffiffiffiffiffi ffi N 0 p gDM 1 ð Þ. Here, N' is the number of non-zero elements in each column of A, corresponding to the number of acquisition steps with non-zero velocity encoding in each direction (thus, N' N). In this case, A has full rank (37) and its pseudoinverse is given by (21) Introducing the orthonormal matrix U ¼ A/h whose column norms are equal to unity, we have (37) and thus