Quantification of Blood Velocity with 4D Digital Subtraction Angiography Using the Shifted Least-Squares Method

BACKGROUND AND PURPOSE: 4D-DSA provides time-resolved 3D-DSA volumes with high temporal and spatial resolutions. The purpose of this study is to investigate a shifted least squares method to estimate the blood velocity from the 4D DSA images. Quantitative validation was performed using a flow phantom with an ultrasonic flow probe as ground truth. Quantification of blood velocity in human internal carotid arteries was compared with measurements generated from 3D phase-contrast MR imaging. MATERIALS AND METHODS: The centerlines of selected vascular segments and the time concentration curves of each voxel along the centerlines were determined from the 4D-DSA dataset. The temporal shift required to achieve a minimum difference between any point and other points along the centerline of a segment was calculated. The temporal shift as a function of centerline point position was fit to a straight line to generate the velocity. The proposed shifted least-squares method was first validated using a flow phantom study. Blood velocities were also estimated in the 14 ICAs of human subjects who had both 4D-DSA and phase-contrast MR imaging studies. Linear regression and correlation analysis were performed on both the phantom study and clinical study, respectively. RESULTS: Mean velocities of the flow phantom calculated from 4D-DSA matched very well with ultrasonic flow probe measurements with 11% relative root mean square error. Mean blood velocities of ICAs calculated from 4D-DSA correlated well with phase-contrast MR imaging measurements with Pearson correlation coefficient r = 0.835. CONCLUSIONS: The availability of 4D-DSA provides the opportunity to use the shifted least-squares method to estimate velocity in vessels within a 3D volume.

2 D digital subtraction angiography and a 3D rotational angiography acquisition followed by 3D-DSA reconstruction are the standards for vascular morphology assessment. There is an increasing demand for hemodynamic information, including blood flow rate and velocity for diagnosis, treatment planning, and evaluation. [1][2][3] Many algorithms have been proposed to estimate blood flow in arteries from 2D-DSA; these are divided into 2 major classes: bolus-tracking and computational methods. [4][5][6][7][8][9][10][11] A thorough review of these techniques can be found in Shpilfoygel et al. 12 However, 2D-DSA based methods are challenged by vessel overlap and vessels being aligned with the primary x-ray beam direction.
Recently, a new reconstruction method [13][14][15][16] has been applied to a single 3D rotational angiography acquisition to generate time-resolved 3D-DSA at frame rates up to 30 frames/second (4D DSA). The availability of the geometric and temporal data in a 4D-DSA reconstruction provides the opportunity to estimate velocity and flow more accurately than has previously been possible with 2D DSA. In this article, a shifted least-squares based technique was used with the 4D-DSA data to estimate blood velocity. The proposed algorithm was first validated using flow phantom studies, in which the velocity could be documented using an ultrasonic flow probe. Blood velocity was also determined in 4D-DSA datasets from 14 ICAs of volunteers who also underwent 3D fast phase-contrast with vastly undersampled isotropic projection reconstruction (PC VIPR) studies. [17][18][19]

4D DSA Acquisition and Data Preparation
Noncontrast (mask) and contrast-enhanced (fill) rotational acquisitions were performed using a conventional flat panel detector angiographic system (Artis zee; Siemens, Erlangen, Germany) with the following settings: 70 kV(peak), 0.36 Gy/frame, 304 images, 260°arc. Contrast injection (iohexol, Omnipaque 300; GE Healthcare, Piscataway, New Jersey) into the subject was started at the beginning of the fill run (3 mL/second, 7-second injection duration). The rotational acquisition protocol yields xray projection images, which are used to automatically reconstruct both 3D-DSA and 4D-DSA with an isotropic spatial resolution of 0.46 ϫ 0.46 ϫ 0.46 mm 3 .
The centerline of each vascular segment was determined from the static 3D volume using an efficient 3D parallel-thinning algorithm. 20 The contrast waveform map, which is the time concentration curves (TCCs) along the centerline path (curvilinear length z), was extracted from the 4D-DSA and stored for each vascular segment for velocity calculation.

Shifted Least-Squares Algorithm
In an arterial contrast injection, there is a temporal oscillation in iodine concentration that arises from the mixing of contrast medium, which is injected at a fixed rate, and nonopacified blood, which flows at a variable rate driven by the cardiac cycle. This temporal variation in contrast, referred to as pulsatility, appears at points downstream of the injection with a time delay. The time delay is related to the distance downstream and the blood velocity. Therefore, velocity can be estimated by measuring the distance along the vessel centerline and the time delay.
For any 2 points, i and j, along the centerline, the shifted leastsquares difference between their TCCs, c i (t) and c j (t), can be calculated as The value ij 0 , ij 0 ʦ [0, T] that minimizes the shifted least-squares difference, ( ij ), is regarded as the time of bolus transport between these 2 points. The spatial distance z between the 2 points is calculated from the 3D path length along the vessel centerline. The time-shift 0 as a function of the spatial distance z can then be fit to a linear relation: where the slope a is the inverse of the velocity v.

Optimizing Waveform Selection
It has been found that the shifted least-squares algorithm provides the most reliable velocity calculation in vessel segments where pulsatility is strong and consistent. To automatically select the waveform regions with strong and consistent pulsatility, we generated a sideband ratio (SBR) map. 21 For each point in the waveform map, a short-time Fourier transform was applied to generate the local power spectrum PS fn , where f n is the n th characteristic frequency. The SBR at this point can be calculated as where the denominator summation includes the neighboring frequency points and excludes the fundamental frequency f 0 . From this SBR map, a low-pass filter was applied and a threshold of one-fourth of the median SBR was enforced to generate a mask. This mask excludes low pulsatile signals from consideration in the velocity calculation.

Phantom Study
To validate the proposed method for velocity estimation, we used a flow phantom containing a loop of plastic tubing to perform 4D-DSA studies. Pulsatile water flow with different mean flow rates was generated using a perfusion pump with a controllable flow setting (S3; Storkcert, Freiburg, Germany). A 5F catheter was introduced into the tubing for injection of iodinated contrast medium (Omnipaque 300 diluted with 25% water) with a power injector. 4D-DSA-derived flow was compared with flow measurements made with an ultrasonic flow probe (400-Series Multi-Channel Flowmeter; Transonic, Ithaca, New York) attached to the tubing 5 inches away from the catheter tip. Rotational projection images were acquired using the mask and iodine-filled scans with 12-second acquisition over a 260°arc. The frame rate was 30 frames per second. The 4D-DSA volume data were reconstructed using the vendor's prototype software (Siemens Healthineers, Forchheim, Germany).
A time-resolved flow profile was recorded during each fill scan using the ultrasonic probe. The reference flow measurement was taken as the average of the flow profile over the same time window that was used for 4D-DSA-derived flow.
Five different flow settings were selected to provide average flow rates ranging from 285 to 1140 mL/min (corresponding to the mean velocity of 15-60 cm/s for the tubing with a 0.64-cm diameter). Three 4D-DSA scans and reference flow measurements were made for each flow setting.

Clinical Study
Nine subjects who had both DSA and PC VIPR examinations were selected from a data base generated under University of Wisconsin review board approval. Because 5 of 9 subjects had bilateral studies, a total of 14 ICAs were available for evaluation. Scanning parameters for postcontrast PC VIPR were the following: FOV ϭ 22 ϫ 22 ϫ 22 cm 2 , TR/TE ϭ 12.5/4.8 ms, velocity-encoded gradient-echo imaging ϭ 80ϳ100 cm/s, bandwidth ϭ 83.3 kHz, readout matrix ϭ 320 points per projections, spatial resolution for the complex difference image ϭ 0.68 ϫ 0.68 ϫ 0.68 mm 3 . The complex difference images from the PC VIPR scan were segmented for vessel content using Materialise Mimics (Materialise NV, Leuven, Belgium). The segmented MR imaging data were then registered to the 3D-DSA using a correlated point registration with manual point placement, followed by an affine rigid registration. Centerline coordinates of the ICAs generated from the 3D-DSA volumes were used for velocity estimation. An automated workflow for the PC VIPR analysis was developed in Matlab (MathWorks, Natick, Massachusetts). This software performed cross-sectional plane placement for every voxel along the 3D centerline path, segmentation of the vessel boundaries in each plane, and averaging of the normal velocities for all points inside the segmented boundaries. The reported velocity value for a vessel was then the average result over all the centerline points that were present in both the PC VIPR and 3D-DSA volumes.

Statistical Analysis
The validation process was performed using linear regression fit between flow probe (V fp ) and V dsa (calculated using 4D-DSA) for the phantom study. Correlation analysis and Bland-Altman plots were used to compare the velocity measurements using PC VIPR versus 4D-DSA for the clinical study. The Pearson r was used for correlation analysis, and a P value Ͻ .05 for the correlation coefficient was considered statistically significant. Figure 1 shows sample data from the phantom study. First, the centerline was generated from the static 3D volume (Fig 1A). In the displayed format, higher blood velocities produce ridges that are closer to vertical, whereas lower blood velocities produce ridges closer to horizontal. The highlighted region is the optimized waveform region for the flow calculation determined using the SBR method. TCCs of 2 selected voxels along the centerline (marked as blue and red stars on    Figure 1D shows the least-squares differences of these 2 signals as a function of the time-shift, where the minimal appears at 0 ϭ 5 is considered as the time of bolus transport from the blue voxel to the red voxel. Figure 1F demonstrates the linear fitting between the time-shift of each voxel along the centerline and its corresponding centerline position. Figure  1E is the flow profile recorded from the flow probe (downsampled to 1/30 second to correspond to the 4D-DSA timeframes). Flow measurement was taken as an average over the same temporal window used in the 4D-DSA flow calculation. The linear regression fit between the flow probe measurement and the 4D-DSA calculation is shown in Fig 2. The linear regression equation is as follows: V dsa ϭ 1.0069V fp ϩ 2.8482, where V dsa is the velocity calculated from 4D-DSA and V fp is the measurement using the flow probe. The determination coefficient R 2 of the fit was 0.9677. The 95% confidence interval slope was 0.8966 -1.1172 and the 95% CI intercept was Ϫ1.4917-7.1881. These evaluations have been obtained with significance levels of P ϭ 4.521E-11 and P ϭ 0.18, respectively. The relative root mean square error (normalized by mean velocity) was 11.3%.

Clinical Study
The average ICA flow rates were calculated from 13 of the 14 ICAs that had 4D-DSA examinations. One subject was excluded due to the poor pulsatility in the TCCs. Figure 3 shows an example of the ICA studies. Average ICA flow rates were also measured with the PC VIPR studies. Correlation between 2 measurements using 4D-DSA and PC VIPR is shown in Fig 4. Average ICA flow rates calculated using 4D-DSA and measured using PC VIPR are well-correlated, with the Pearson r ϭ 0.835 and a significance level P ϭ .0002. Figure 5 is the Bland-Altman analysis of 4D-DSA compared with PC VIPR, which shows all points within 2 SDs, with a bias of 6.5 cm/s and 2-SD limits of agreement of Ϫ3.56 cm/s and 16.5 cm/s.

DISCUSSION
The shifted least-squares algorithm is a simple-but-robust method to calculate flow rates by finding the bolus transport time between any 2 points along the centerline and performing linear regression with the corresponding positions. The combination of the spatial and temporal data available in the 4D-DSA reconstruction provides the opportunity to apply this technique to measurement of blood flow in routine 4D-DSA acquisitions. The flow phantom study demonstrated that by using this method, we were able to calculate the flow rates with 11% relative root mean square error compared with ultrasonic flow probe measurements. Results from the 13 human ICAs showed that the calculated velocities from 4D-DSA studies correlated well with the measurements  using the 3D phase-contrast MR imaging method. We believe that the reason for the higher velocities in the 4D-DSA studies compared with the PC VIPR measurements (6.5 cm/s bias) is the increment in flow caused by the intra-arterial power injection of contrast medium.
The proposed method calculates the mean flow rate within a selected window defined by the SBR map. This truncated average of pulsatile flow varies depending on the window size and location. This variation can be demonstrated in Fig 6A and B, where the average of the flow profile was performed with a gradually expanding average window. The truncated average oscillates depending on where the truncation window ends and eventually converges to the mean signal as the window expands. Similar oscillation has been observed in the flow calculation using 4D-DSA. Figure 6D shows the estimated average velocity changes as the right side of the window gradually expands as shown in Fig 6C. The current algorithm automatically selects a window purely on the basis of the SBR map so that the selected window range may fall anywhere between 4 and 6 pulsatile cycles. Further criteria could be included to trim the window to an integer number of pulsatile cycles to minimize this error.
The proposed algorithm worked best when there was good and consistent pulsatility in the 4D-DSA TCCs. With each cardiac cycle, as an intra-arterial contrast bolus is injected, there is an oscillation in the ratio of contrast medium to nonopacified blood. This pulsatility can be identified and extracted from the contrast waveform map by thresholding the SBR map. For vessels more distal to the contrast injection site where the pulsatility can be low, injection protocols should be modified to enhance the pulsatility strength or contrast kinetics. The present work focused on the evaluation of blood velocity within the ICAs. Investigation of more downstream segments, including the MCAs and basilar arteries, should be performed in the future.

CONCLUSIONS
4D DSA is a novel technique that provides both the geometric and temporal information required to measure blood flow velocity and flow. In this study, calculated velocity values from 4D-DSA TCCs correlated well with human ICA measurements derived from phase-contrast MR imaging and flow phantom measurements using an ultrasonic flow probe.