Quantification of Internal Carotid Artery Flow with Digital Subtraction Angiography: Validation of an Optical Flow Approach with Doppler Ultrasound

BACKGROUND AND PURPOSE: Digital subtraction angiography is the reference standard technique to evaluate intracranial vascular anatomy and used on the endovascular treatment of vascular diseases. A dedicated optical flow-based algorithm was applied to DSA to measure arterial flow. The first quantification results of internal carotid artery flow validated with Doppler sonography are reported. MATERIALS and METHODS: We included 22 consecutive patients who underwent endovascular procedures. To assess the sensitivity of the algorithm to contrast agent-blood mixing dynamics, we acquired high-frame DSA series (60 images/s) with different injection rates: 1.5 mL/s (n = 19), 2.0 mL/s (n = 18), and 3.0 mL/s (n = 13). 3D rotational angiography was used to extract the centerline of the vessel and the arterial section necessary for volume flow calculation. Optical flow was used to measure flow velocities in straight parts of the ICAs; these data were further compared with Doppler sonography data. DSA mean flow rates were linearly regressed on Doppler sonography measurements, and regression slope coefficient bias from value 1 was analyzed within the 95% confidence interval. RESULTS: DSA mean flow rates measured with the optical flow approach significantly matched Doppler sonography measurements (slope regression coefficient, b = 0.83 ± 0.19, P = .05) for injection rate = 2.0 mL/s and circulating volumetric blood flow <6 mL/s. For injection rate = 1.5 mL/s, volumetric blood flow <3 mL/s correlated well with Doppler sonography (b = 0.67 ± 0.33, P = .05). Injection rate = 3.0 mL/s failed to provide DSA–optical flow measurements correlating with Doppler sonography because of the lack of measurable pulsatility. CONCLUSIONS: A new model-free optical flow technique was tested reliably on the ICA. DSA-based blood flow velocity measurements were essentially validated with Doppler sonography whenever the conditions of measurable pulsatility were achieved (injection rates = 1.5 and 2.0 mL/s).

E ven though digital subtraction angiography has traditionally been confined to standard anatomy assessment, quantification of blood flow based on DSA is becoming an important topic that could help neurointerventionists in making adequate periprocedural decisions. Some reports have described new techniques based on DSA that are able to assess flow or flow changes during treatment of stented aneurysms. 1,2 However, the development of clinically useful tools based on the integration of engineering, hemodynamic and physiologic knowledge still requires improved translation of biofluid mechanical information into clinical applications. 3 X-ray video densitometry, based on the detection of the displacement of radiopaque contrast material through the vascular system, has been studied since the early 1960s and has been divided into 2 main classes: tracking and computational methods. 4 Sarry et al 5 estimated the flow by using an inverse advection model. Bogunović and Loncarić 6 proposed the combination of DSA and 3D rotational angiography (3DRA), using an analysis of the time-attenuation curves. Rhode et al 7 developed a model-based and weighted-optical field (OF) approach to improve already existing techniques and compared the results with simulation data, while Imbed et al 8,9 developed a similar approach on the femoral artery and simulated angiographic data. Waechter et al 10 developed a model-based approach to measure flow in the cerebral arteries.
These approaches were all limited in case of fast flows, required longer straight-vessel segments, and were affected by low signal-to-noise ratios. A dedicated algorithmic scheme was developed to reduce instabilities due to temporal and spatial noise and to cope with fast flows, to address these issues. 11 Essentially, we used only the modulation by the cardiac cycle to extract flowvelocity values. The OF principle was then applied to the pulsating component of dye-concentration signal. 11 In this article, the first clinical results of the proposed algorithm applied to the ICA are reported. The clinical implementation and verification against Doppler ultrasound (USD) data in a consecutive cohort of patients are described, and the limits of the current technique are discussed.

Clinical Data
The study protocol was approved by our institutional ethics committee (NEC 07-056). Between 2008 and 2009, all patients undergoing angiographic investigation or endovascular treatment for intracranial aneurysms in our department were considered for recruitment. Patients with intraventricular drainage, hydrocephalus, or hematoma were excluded because of a possible unstable intracranial pressure that could interfere with intracranial flow during the procedure and image acquisition. We included prospectively 22 consecutive patients with ruptured or unruptured intracranial aneurysms. All patients were evaluated under general anesthesia. A 5F 0.038-inch diagnostic catheter (Cook, Bloomington, Indiana) was placed selectively in the ICA (3 cm after the common carotid bifurcation) by using a femoral approach. The angiographic examinations were performed with a monoplane angiographic C-arm Allura FD20 system (Philips Healthcare, Best, the Netherlands). Every patient underwent 3DRA to acquire the 3D vascular anatomy with an injection rate (IR) ϭ 3 mL/s during 6 seconds. An optimized projection, offering the best view of the vessel segment, was chosen for flow-rate measurement with the fewest overlapping vessels. High-speed DSA (60 frames/s with the 1024 ϫ 1024 detector pixel matrix) was acquired in 10 seconds. The radiation exposure was equivalent to 2-frames-persecond sequences in terms of dose level. Two or 3 DSA acquisitions by using different contrast agent IRs were obtained for each patient. The contrast agent (CA) was iopamidol (Iopamiro 300; Bracco, Milan, Italy), and the injection duration was 3 seconds. CA was injected with an Imaxeon Avidia angiographic contrast injector (Medrad, Indianola, Pennsylvania). IRs were randomly selected before each acquisition among 3 protocols: 1.5 mL/s (n ϭ 19), 2.0 mL/s (n ϭ 18), and 3.0 mL/s (n ϭ 13).

Optical Flow Principle
The OF-based algorithm used in this study to measure internal carotid flows was based on both 2D DSA and 3DRA image data. 11 The blood-contrast flow-velocity fields were measured from 2D DSA sequences. The injected contrast material was diluted in the blood stream while being transported through the vascular network. The pulsating contrast patterns generated by dye injection under the effect of the cardiac cycle were captured in DSA sequences. The main improvement provided by this approach consisted of extracting only the cardiac modulation from the signal, while discarding the low-frequency wash-in and washout compo-nents. The OF approach is a mass-conserving technique, which was applied in an image-processing algorithm to extract velocity and flow curves. The second step consisted of projecting the DSA sequence on the arterial centerline extracted from the 3DRA data to generate the 2D "contrast wave map" or "flow map," which defines the effective progression of the contrast wave along the arterial axis. Then, a 1D OF scheme was applied iteratively on the contrast wave map to extract the velocity profile. The iterative process was stopped when the quadratic error was no longer decreasing. The high-frame-rate acquisition was allowed to reduce constraints on the required arterial segment length and, therefore, to address more complex arterial geometries such as in cerebral vessels. For average velocities as high as 60 cm/s for instance, the required arterial segment length was approximately 20 mm to allow a proper warping, which was largely achieved at a frame rate of 60 Hz (2 times 10 mm). If C (t, z) is the physiologically modulated contrast attenuation at instant time t, and z is the curvilinear axial position along the centerline, displacement velocity V (t, z) satisfies the differential contrast wave equation:
If S (z) is the arterial cross-section at z, the volume flow curve is given by where L is the curvilinear length of the arterial segment. S(z) is measured from the segmented vasculature in the 3DRA dataset. Figure 1 summarizes the whole velocity-extraction process, 11 starting from the DSA and 3DRA acquisitions ( Fig 1A). The centerline of the artery was used to determine the curvilinear path represented by z. The contrast-level intensity along the curvilinear length z, namely the flow map or contrast wave map, 10,12 is represented in Fig 1B. The velocity profile extracted from the flow map along with the volumetric flow is represented in Fig 1C. To ensure good mixing between CA and blood, we selected the region of interest in the ICA at least 10 arterial diameters from the injection site. 13 The region of measurement was reachable by submandibular USD examination.
The presence of physiologic pulsatility in the contrast material motion pattern is effectively required for our OF algorithm 11 and was correlated with injection rate quotient (IQR) metric I q , where Q i is the IR (milliliters per second) and Q D is the timeaverage USD volume flow measure (milliliters per second). It has been shown, by comparing cardiac cycles measured from both electrocardiogram and time-intensity curves, that for I q Յ 1.2, the physiologic pulsating character is conserved in a CA motion pattern. 14

Doppler Measurement (USD)
The ICAs were examined by an experienced investigator with a 5to 11-MHz color-coded duplex linear probe connected to an Aplio 770A scanner (Toshiba Medical Systems, Tokyo, Japan) by using a standard technique. 15 The sample volume, adjusted to the size of the artery, was placed in the ICA with a submandibular approach, to obtain velocity waveforms. Its placement was adjusted around the site of measurements to acquire the higher color-coded signal. The angle-corrected mean velocity, peak systolic velocity, and end-diastolic velocity were determined by tracing the maximum frequency envelope of the waveform. The reference flow curves, based on Womersley solutions, were extracted with a Matlab (MathWorks, Natick, Massachusetts) program from images of velocity profiles acquired with USD. The vessel diameter needed to estimate the volumetric flow rate was measured on the 3D reconstruction. For cases in which several USD measures were acquired, repeatability was estimated.

Statistics and Regression Analysis
The validation process was performed by using linear regression fit between Q D (Doppler) and Q x (x-ray). Q D was either a single USD time-average flow rate measure (n ϭ 6) or an average of several measures (n ϭ 16). Q x was the time-average measure provided by the OF method. A robust fitting technique was used to reduce the impact of outliers by iteratively reweighting least squares with a bisquare weighting function (r ϭ 4.6) 16,17 4) The slope b and offset a were the regression coefficients of the robust fit model. The analysis was made for every injection rate group separately (1.5, 2.0, and 3.0 mL/s) and for several Q D thresholds (Q D Յ 3.0, Յ4.0, Յ5.0, Յ6.0, Յ7.0 mL/s) to test the volume flow limits of the current OF technique prototype. Ideally, the 95% CI of the regression slope estimate should contain the value 1. The y-intercept a coefficient was proportional to the CA injected volume flow and was not taken into account in the validation process because of the lack of knowledge of this factor. 18 Nonparametric Wilcoxon rank sum U tests and the Kruskal-Wallis H test were used to test statistics for the continuous data: b-slope regression coefficient, injection rate quotient, and Doppler measurement precision. If the H test was significant among 3 groups, the U test was performed to find out which of the differences were statistically significant. The significance level was fixed to ␣ ϭ .05.

RESULTS
The On-line

Sampling Effect on the Regression Slope Coefficient
The comparison of the regression slope coefficient among Q i groups ( Fig 3A) indicated that the only significant difference was between the Q i ϭ 2.0 mL/s and Q i ϭ 3 mL/s groups (P ϭ .0015). The mean slopes, if sampling errors were ignored, were 0.48, 0.75, and 0.35 for Q i ϭ 1.5, 2.0, and 3.0 mL/s, respectively, with SE ϭ 0.06. The corresponding slope upper values for the respective Q i values were 0.84, 1.21, and 0.72, respectively, with a SE of 0.12 ( Fig  3B). The Q i ϭ 3 mL/s group was significantly biased compared with both Q i ϭ 1.5 mL/s and 2.0 mL/s (P ϭ .039). For each CA injection rate Q i , the sampling effect on slope driven by Q D cutoffs at 3.0, 4.0, 5.0, 6.0, and 7.0 mL/s is plotted in Fig 4A. The error bars are superimposed for the group Q i ϭ 2.0 mL/s. For this group, the slope bias was not significant up to Q D ϭ 6 mL/s. Above Q D ϭ 5 mL/s, flow quantification was not sensitive to sampling. For Q i ϭ 1.5 mL/s, slope bias was significant (P Ͻ .05) for sampling cutoffs of Ͼ3 mL/s; and for Q i ϭ 3.0 mL/s, the slope bias was effective for the whole sample (P ϭ .0004).

Injection Rate Quotient
The ratio between the injection rate and the mean blood flow is illustrated in Fig 4B. The groups Q i ϭ 1.5 and 2.0 mL/s were not different (P ϭ .85), while both of them were different from the  A, Slope versus Q D for injection rates Q i ϭ 1.5 (blue), 2.0 (red), and 3.0 mL/s (green). The 95% CI limits have been superimposed to Q i ϭ 2.0 mL/s points. The graph shows that the slope bias calculated is not significantly different from 1 up to Q D ϭ 6 mL/s (P Ͼ .05). The Q i ϭ 1.5 mL/s group (blue points) follows the same evolution trend as Q i ϭ 2.0 mL/s, though with accentuated bias. A nonsignificant slope bias is likely for the Q D range below 3 mL/s (P ϭ .06). B, Injection rate quotients for Q i ϭ 1.5 mL/s (group 1), 2.0 (group 2), and 3.0 mL/s (group 3). The blue bars correspond to quotients Ͻ1.2, and the red bars, to quotients Ͼ1.2. The fact that the Q i ϭ 1.5 and 2.0 mL/s groups are reconstructed in similar ways with minimum slope bias, unlike the Q i ϭ 3.0 mL/s group, is partly due to differences in injection rate quotients (P ϭ .01).
Q i ϭ 3.0 mL/s group (P ϭ .01). In the first 2 groups, 90% of the data points had an injection rate quotient Ͻ 1.2, implying that physiologic information and modulation were conserved in contrast-motion patterns. 14 In the Q i ϭ 3.0 mL/s group, 45% of the points had an injection rate quotient Ͼ 1.2. Therefore, the optical flow measures significantly matched the Doppler outcome for 2.0 mL/s injected cases when circulating blood flow ranged between 1 and 6 mL/s. With a 1.5 mL/s injection rate, flow measures matching was for circulating blood flow ranging from 1 to 3 mL/s. Figure  5 shows, for case 13, the contrast wave maps and flow curves corresponding to IR ϭ 1.5, 2.0, and 3.0 mL/s. As one can notice, the contrast wave maps corresponding to IR ϭ 1.5 and 2.0 mL/s show the best premises for velocity extraction (several pulses) compared with IR ϭ 3.0 mL/s (1 detectable pulse).

Quantification of the Flow through DSA
The DSA-based approach represents a logical evolution of the criterion standard method for intracranial vascular anatomy assessment. Although DSA is an imaging technique used for endovascular treatments, the lack of tools for periprocedural quantitative flow measurements did not ease real-time therapeutic decision-making. 1,2 The advantages for measuring blood velocity and flow changes during the procedure are numerous. For instance, one could evaluate the hemodynamic effect of an angioplasty in a stenotic segment for vasospasm or intracranial stenosis or the impact on the intra-aneurysmal hemodynamics of a flowdiverter stent to decide whether to place another one or coils. 1,2 Other possible applications could be in acquiring patient-specific boundary conditions (flow curves) for computational flow dynamics. Different methods have been proposed to quantify both flow changes and characteristics from DSA sequences. 3,7,[10][11][12]19 However, most could not provide a reasonable level of confidence for clinical applications because of their inability to fully cope with both the complex advective and diffusive dispersion of the dye in blood stream and the complex pulsatile nature of the flow in the artery. 3 Additionally, not all studies paid attention to CA injection protocols that might alter the physiologic modulation of CA motion waveform. 13,[20][21][22] This physiologic modulation is a key factor in dealing with the mixing problems related to CA transport in quantifying blood flow. The extraction of the modulated component with proper filtering is a necessary step that will condition the performance of such a functional algorithm. Con-trary to the technique of Waechter et al, 10 which also extracted the modulated component, our OF algorithm is model-free and can be applied to quantify flow in both the vessels and aneurysms. 1 In our study, we used 3 different IRs: 1.5, 2.0, and 3.0 mL/s. The Craya-Curtet number was 0.5-0.9 with a 5F catheter of 0.035-inch inner diameter, the tip of which was positioned at least a 10-artery diameter below the measurement zone. 13 Measurable flow pulsatility was obtained for IR ϭ 1.5 mL/s and 2 mL/s (injection rate quotient Յ 1.2). 13,14 Consequently, the OF technique performance was reasonably high, with slightly better results for 2 mL/s and for volume flows of Ͻ6 mL/s. The performance with IR 1.5 mL/s was lower than expected because of the larger dilution of the dye with a deleterious effect on signal-to-noise ratio. The injection at 3 mL/s presented excessive reflux at the injection site, which altered the matching between dye and blood-flow dynamics in the region of interest.

General Trend of the OF Method
The general trend of our current OF algorithm was to measure flow with an underestimated slope bias that depends on the injection rate: 0.48, 0.75, and 0.35 for Q i ϭ 1.5 mL/s, 2.0 mL/s, and 3.0 mL/s, respectively, with RMSE ϭ 0.06. Those values do not take into account sampling errors. The maximum slope value variations (95% CI) over the whole Q D show that the slope bias for Q i ϭ 1.5 and 2.0 mL/s was not statistically significant (the maximum slope is, respectively, 0.84 and 1.21, with RMSE ϭ 0.12, P ϭ .05), unlike that for IR ϭ 3.0 mL/s group. This general trend is consistent with the fact that the performance of the OF method depends first on the presence of the physiologic information on contrast wave maps and, consequently, on the injection rate quotient (Յ1.2).
Thorough examinations of slope dependency with Q D show that the slope bias for group IR ϭ 2.0 mL/s was not statistically significant (P Ͼ .05) up to 6 mL/s, setting the limit on reliable velocity measurement to 60 Ϯ 6 cm/s, taking into account the average inlet diameter of 4.8 mm and the relative RMSE ϭ 0.12. For the IR ϭ 1.5 mL/s group, the slope bias was small, within Q D Յ 3.0 mL/s, and increased with Q D . As Q D increases, the corresponding CA signal-to-noise ratio is reduced, possibly because poor mixing occurs at the injection site due to high Craya-Curtet number. 7,13 Consequently, intensity gradients in 2D images would have been overestimated in the denominator of the OF equation. More improvements to the OF technique are expected in the near future to cope uniformly with different injection rates and flow rates.

Study Limitations
Doppler Measurements. The OF trend and limitations presented in the above section are discussed under the assumption that USD measures are infinitely accurate, which is, in fact, not correct. Biases due to geometric spectral broadening and to angulations between the beam direction and the vessel axis may lead to an overestimate of the USD velocity of approximately 30%. 23 The Doppler measurement repeatability inherent to the experienced investigator was reasonable (Ͻ23%). However, the mean-flowrate repeatability of the measurements was Ͻ15% in 13 of 16 cases that had at least 2 USD measures, and for the other 3, it was Յ23%. The volumetric flow shapes were similar for repeated measures.
Acquisition Protocol. Sources of mismatch between OF and USD could have been due to the delay between the 2 acquisitions or to inaccurate estimates of the parent artery diameter from the 3D geometry. The diameter is influenced by the threshold selection and by the location of the measurement on the artery, which could deviate from the USD measurement zone. In this work, we focused on the most rectilinear part of the internal carotid artery to ensure the reliability of USD velocity measurements. In the submandibular part of ICA, the vessel narrowing is negligible within the OF measurement-segment length (2 cm). The range of flow rates in the ICA is quite exhaustive to validate the method. Nevertheless, it will be necessary to validate in other cerebral vessel locations where flow rates are smaller.

CONCLUSIONS
The DSA-based OF approach can be used to measure circulating flows in ICAs. The OF measure was validated with USD for flows up to 6 mL/s with CA injected at 2.0 mL/s and up to 3 mL/s with CA injected at 1.5 mL/s. The values obtained from the IR ϭ 3.0 mL/s group were significantly biased because only partial dynamic information was propagated to the CA motion pattern. Improvements to the OF technique are expected in the near future to improve the accuracy and uniformity of velocity measurements over ranges of flow rates and injection rates. The capability of using this technique during the procedures and assessing fast flows in cerebral vessels would open new insights in intraoperative therapeutic decision-making.

Statistical Analysis
The linear regression method we used applies generally to cases in which the independent variable (Q D ) is free from errors. In such situations, the robust fitting method consistently rejects the outliers and focuses on the bulk attenuation of data points. In our case, because 16 of 22 points were characterized by an estimated error on Q D , errors-in-variables regression methods would have been used (eg, Deming regression) to avoid underestimating the regression coefficients. Some OF and Doppler measures were unique, and using Deming-like regression would have led to speculations on relative variances. Requiring that all data points have more than a single measure simultaneously in either OF or Dopp-ler would have reduced the number of points by 50%. The robust regression fit is then a reasonable compromise to capture trends of the optical flow method.
The H test is generally performed among groups that are statistically independent. In our study, Q D cutoff-based subsampling created the 3 IR (Q i ) groups. However, with such a low number of points (5) in each group, the lack of that statistical independence should not have dramatic consequences, given the level of precision and accuracy that we were dealing with.

USD
For several decades, the potential of Doppler instruments to measure volume flow rates at the same time as velocity waveforms has been discussed. 24 One article has reviewed all of the factors of inaccuracy and provided the orders of magnitude of the errors that were derived from the Doppler acquisition constraints and that were evaluated during validation tests. 23 Validation experimentation estimates the performance of the technique with in vivo and in vitro tests. 25 The authors found errors in the range of 5%-54% for in vivo experiments and errors Ͻ18% for in vitro flow-phantom experiments.