Fully distributed absolute blood flow velocity measurement for middle cerebral arteries using Doppler optical coherence tomography

: Doppler optical coherence tomography (DOCT) is considered one of the most promising functional imaging modalities for neuro biology research and has demonstrated the ability to quantify cerebral blood flow velocity at a high accuracy. However, the measurement of total absolute blood flow velocity (BFV) of major cerebral arteries is still a difficult problem since it is related to vessel geometry. In this paper, we present a volumetric vessel reconstruction approach that is capable of measuring the absolute BFV distributed along the entire middle cerebral artery (MCA) within a large field-of-view. The Doppler angle at each point of the MCA, representing the vessel geometry, is derived analytically by localizing the artery from pure DOCT images through vessel segmentation and skeletonization. Our approach could achieve automatic quantification of the fully distributed absolute BFV across different vessel branches. Experiments on rodents using swept-source optical coherence tomography showed that our approach was able to reveal the consequences of permanent MCA occlusion with absolute BFV measurement. found that the BFV ranged from − 23.04% to + 33.97% across all four branches at baseline but was reduced to − 19.7% to + 30.88% after the occlusion. The proposed absolute BFV measurement approach shows good potential for further investigation of hemodynamic changes during permanent MCA occlusion.

Abstract: Doppler optical coherence tomography (DOCT) is considered one of the most promising functional imaging modalities for neuro biology research and has demonstrated the ability to quantify cerebral blood flow velocity at a high accuracy. However, the measurement of total absolute blood flow velocity (BFV) of major cerebral arteries is still a difficult problem since it is related to vessel geometry. In this paper, we present a volumetric vessel reconstruction approach that is capable of measuring the absolute BFV distributed along the entire middle cerebral artery (MCA) within a large field-of-view. The Doppler angle at each point of the MCA, representing the vessel geometry, is derived analytically by localizing the artery from pure DOCT images through vessel segmentation and skeletonization. Our approach could achieve automatic quantification of the fully distributed absolute BFV across different vessel branches. Experiments on rodents using swept-source optical coherence tomography showed that our approach was able to reveal the consequences of permanent MCA occlusion with absolute BFV measurement.

Introduction
Doppler Optical Coherence Tomography (DOCT) [1][2][3][4][5][6], or Optical Doppler Tomography, is a functional imaging modality based on Optical Coherence Tomography [7] which is used to image optical scattering media with high resolution based on light interference. First developed in the late '90s, DOCT takes advantage of Doppler effects happening at the intersection of light and moving substances (e.g., flowing particles in blood vessels) that result in shifts in the phase component of the back-scattering light. Leveraging this phase shift (Doppler shift), the velocity of such moving substances can be derived. DOCT has been used extensively for ophthalmology [8][9][10][11], cardiology [12][13][14], and dermatology [15][16][17][18][19] due to its flexibility to detect blood flow velocity at a high accuracy.
Recently, DOCT has been shown to be able to image cerebral blood vessel networks in the brain cortex and to reveal changes in blood flow velocity (BFV) within the brain cortex during different types of brain damage such as dura removal [20], cocaine administration [21], and ischemia induced by photodynamic therapy (PDT) [22]. Moreover, quantitative cerebral BFV measurement was investigated using DOCT to reveal better functional information from brain imaging [23][24][25][26][27][28][29][30][31]. Flow velocity is a unique parameter, distinct from the flow rate, because it is independent of the vessel diameter and area. Given the laser central wavelength λ and the refractive index n of the scattering media, the absolute velocity of the moving particles can be described by where f is the Doppler phase shift, ΔT is the time interval between two samplings, and θ is the Doppler angle. Equation (1) indicates that the velocity measurement not only relates to the properties of the laser and the scattering particles, but it also relates to the geometry of both directions of the laser beam and the flow. Therefore, one of the major difficulties with measuring absolute BFV is to determine the Doppler angle at each position along the vessel since this angle is continuously changing with respect to the vessel geometry, as is shown in Fig. 1 As one of the major arteries that supply blood to the cortex, the middle cerebral artery (MCA) plays an important role in brain activity and physiology. Investigation of blood flow velocity and its changes in MCA is of significant interest for revealing middle cerebral artery syndromes such as ischemic stroke [32,33]. Focusing on the analysis of cerebral hemodynamics, this paper presents a method to quantify the total absolute blood flow velocity in MCA based on volumetric vessel reconstruction from pure DOCT images. A modified region growing segmentation method is first used to localize the MCA on successive DOCT B-scan images. Vessel skeletonization, followed by an averaging gradient angle calculation method, is then carried out to obtain Doppler angles along the entire MCA. Once the Doppler angles are determined, the absolute blood flow velocity of each position on the MCA is easily found. Given a seed point position on the MCA, our approach could achieve automatic quantification of the fully distributed absolute BFV. Experiments on rodents based on sweptsource OCT demonstrated the feasibility of our approach and revealed the complex cerebrovascular physiology during permanent MCA occlusion.

Animal preparation
To conduct the absolute BFV measurement studies, three 295-400g male Sprague Dawley rats were used. The animals were injected intraperitoneally with a Nembutal bolus [55 mg/kg body weight (b.w.)] at the beginning. Supplemental injections of Nembutal (27.5 mg/kg b.w.) were given as necessary. After resection of soft tissue, a 6.5 × 8mm imaging area of the skull over the left primary somatosensory cortex (rostromedial corner positioned ~1 mm caudal and 2 mm lateral from the bregma) was thinned to ~150 µm using a dental drill. Body temperature was measured via a rectal probe and maintained at 37°C by a self-regulating thermal blanket.

The swept-source OCT system
A swept-source OCT system was built to demonstrate our proposed method. The schematic of the system setup is shown in Fig. 2. A swept source laser with a central wavelength of 1310 nm, a maximum A-line speed of 50 kHz, and a total average power of 16 mW (SSOCT-1310, Axsun Technologies Inc., Billerica, Massachusetts) was used in our system. The bandwidth of the laser source was around 80 nm and the axial resolution of the system was 9.3 μm in air (6.6 μm in tissue). A Michelson type interferometer was used, and 90% of the laser power was sent to the sample arm and 10% of the light to the reference arm. In the sample arm, a fiber collimator, a two-axial galvo mirror scanning system, and an achromatic doublet was used to deliver the beam to the rat brain. The fiber collimator gave a beam diameter of 2 mm, and with the 30 mm focusing length doublet, the system gave a lateral resolution of 14.6 μm. C + + platform based data acquisition software running on a 64 bit Windows 7 operation systems was used to control the galvo scanner and acquire the data.

Algorithm overview and pre-processing
As shown in Fig. 3, taking phase-resolved Doppler images captured with different time intervals as the input, the whole process of our absolute BFV measurement method is as follows: first, due to the limited dynamic range of the Doppler method, the input images are averaged to produce high-contrast vessel profiles; second, a restricted region growing algorithm is adopted to locate the MCA branches; third, the segmented MCA structure is skeletonized and smoothed after which the Doppler angle is calculated according to the geometry of the vessel skeleton. Finally, incorporating the acquired Doppler angles, the absolute BFVs at each location are calculated from the mean vessel signals within the segmented lumen through Eq. (1). As the first step, the raw input DOCT data are produced by the Phase-resolved Color Doppler (PRCD) method [4,15]. The phase-resolved Doppler images are obtained by differencing the phase components of adjacent A-lines such that they represent the average Doppler frequency: where φ j + 1,z and φ j,z are the phases for the signal at depth z of the j + 1th and jth A-line. Using the autocorrelation algorithm and letting Re(A j,z ) and Im(A j,z ) be the real and imaginary part of the complex data, respectively, Eq. (2) can be rewritten as [15] 1, , , In our implementation, after the calculation of all the Doppler phase shifts in a B-scan frame, the final DOCT image was filtered with an average filter (window size: 4-by-4) to smooth out the noise.
Using Eq. (1) and Eq. (3), the absolute value of the flow velocity could be determined given the Doppler angle. However, subject to the arc-tangent operation in Eq. (3), the acquired Doppler phase shift suffers from limited dynamic range (Δφ: [-π, π]), resulting in phase wrapping artifacts if the time interval ΔT is not accommodated to the flow velocity. Moreover, since the BFV and the Doppler angles at different positions along the MCA are different, not all the vessel profiles are prominent in different B-scans. As illustrated in Fig. 4, MCA branches at frame 123 can be seen clearly in ΔT = 2t (t is the minimum time interval between two adjacent A-lines), while the other branches at frame 180 can only be distinguished in ΔT = 20t. Therefore, in order to provide a solid foundation for vessel segmentation in the next step, images at different time intervals are acquired, normalized and averaged as a pre-processing procedure. Two advantages could be realized by doing so: all the vessels could be clearly seen in the averaged cross section since the contrast between the vessel and the background tissue is higher than the non-averaged profiles; and the speckle noise, along with the bulk motion artifacts, is suppressed from averaging the data. Moreover, before proceeding to the next step, each B-scan of the averaged OCT volume is further filtered by a 2D median filter. Figure 4 depicts the vessel profiles at the averaged frames from which can be seen the prominent vessels. Fig. 4. The acquired frame 123 and frame 180 of a total of 256 frames from a rat brain imaging experiment. The limited dynamic range of the resolved Doppler phase can be seen clearly: the cross-sectional profiles of different vessels are only prominent at particular time intervals while in the averaged profile, these vessels all stood out from the background.

Vessel segmentation
Due to its simplicity and fast speed in segmenting regions with well-defined edges, the region growing algorithm is widely employed in different image processing applications from object detection to CT/MRI structure segmentation [34][35][36][37]. Given one or multiple initial seed points in multi-dimensional images, this algorithm examines the neighboring pixels of these seed points and determines whether they should be added to the region. However, this algorithm commonly suffers from the "leaking" problem in images with blurred edges [38] resulting in odd regions with a huge discontinuity. Meanwhile, the vessel signals in the DOCT images are affected by the Doppler angles and the flow pulses, resulting in discontinuities in boundary pixel intensities. To overcome these problems and effectively segment the MCA from stacked DOCT images, a restricted region growing algorithm is proposed. The algorithm is named "restricted" because the growing process is constrained simultaneously by the adaptive threshold and a circular shape prior.
1) The adaptive threshold t is given by the compounding of the mean gray value t 1 of the boundary at the previous frame and a hard threshold t 2 derived from a Region-of-Interest (ROI) at a current frame: t = (t 1 + t 2 )/α, where α is the compounding ratio. In this study, the ROI at a current frame is given by the neighboring region of the previous vessel center, and threshold t 2 is defined as a specific portion of the maximum gray value in this ROI.
2) The circular shape prior is a circular boundary where the vessel region will stop growing when they encounter each other. The center c of the circular boundary is given by the previous vessel center coordinates while its radius r is derived from the mean vessel radius at the previous frame. Moreover, due to the fact that the diameter of the MCA decreases when moving from the main vessel to side branches, the diameter of the shape prior at the next frame will be reduced when the maximum diameter of the current vessel is less than the diameter of the shape prior of the previous frame, resulting in a smaller shape prior at all successive B-scans. By introducing the adaptive threshold which incorporates the gray scale information across sequential frames, the boundary pixel intensity fluctuation caused by changes in Doppler angle and blood flow pulse is overcome; by constraining the growing process with a shape prior, a 'leaking problem' can be prevented. Figure 5 shows how the segmentation of the vessel lumen is restricted by the aforementioned constraints. The growing process starts with a seed point on a branch of the middle cerebral artery and grows at B-scan-wise. The stopping criteria are (1) no pixel satisfies the constraints, (2) the variance of the current ROI is smaller than a pre-defined value and (3) the displacement between the vessel centers of adjacent frames is larger than a threshold. For clarity, Algorithm 1 summarizes the entire process of the proposed method. Initialized parameters: t 1 0 , t 2 1 , α, r 0 , and k; 2 while stopping criteria not met do 3 Get the k frame of DOCT B-scan images; 4 Update the adaptive threshold t k = (t 1 k-1 + t 2 k )/α; 5 Update the shape prior parameters: center c k [x k-1 ,y k-1 ], radius r k ; 6 Update the k frame region seed points by propagating the k-1 frame regions that satisfies the constraint t k , to the current frame; 7 Region grows from the seed points until all constraints are met; 8 Get the stopping criteria status; 9 k = k + 1; 10 end.

Doppler angle correction and absolute BFV calculation
The segmentation of the MCA finds the vessel centroids for each frame at the same time which is composed of a volumetric vessel skeleton that represents the geometry of the vessels. With this MCA skeleton, the Doppler angle at each position can be derived accordingly. Here, we present an averaged gradient method which takes the mean gradient value of the neighborhood ± w of each point as the flow direction at that position. Since the B-scan images are produced discretely, the gradient direction of a point P k = [x k , y k , z k ] T on the MCA skeleton can be represented in a discrete form as: where k = w + 1, w + 2, w + 3,…, then the Doppler angle at point P k is given by: where η = 1 for forward flow and η = −1 for reverse flow. Figure 6 depicts the discrete Doppler angle calculation principle at k position. A larger w will result in a smoother Doppler angle distribution (a smoother vessel skeleton), but also will neglect small angle changes along the vessel. In our implementation, w = 4 was chosen to achieve the best balance of such a tradeoff. The proposed averaged gradient method ensures the robustness and continuity of the calculation of the Doppler angles along the MCA skeleton which is crucial for determining the absolute BFV. After the determination of the Doppler angle, the images taken at the appropriate delay, which results in minimum phase wrapping within the vessel lumen, are selected for the calculation of the absolute BFV. The BFV in one B-scan image is considered uniform, and thus the BFV is measured as the average intensity within the vessel lumen at this frame. The absolute BFV is obtained by dividing this averaged intensity by the cosine value of the Doppler angle at the same location.

Phantom experiment
A flow phantom experiment was carried out to verify the accuracy of our method. 2% intralipid was used as the phantom. The pumping rate was set to 0.047mL/min, resulting in an average speed of 27mm/s using a microtube with an inner diameter of 254µm. The tube was intentionally bent in order to get a varying Doppler angle distribution. The cosine values of the Doppler angles along the flow phantom ranged from −0.354 to 0.307, equivalent to the angles of 72.12 to 110.75 degrees. During the calculation of the flow velocity, the phantom was segmented from the DOCT images using the proposed restricted region growing algorithm, and the Doppler angle distribution was calculated afterward (as shown in Fig. 7). However, because the Doppler angle was near 90 degrees in the center portion of the tube (resulting in a very small cosine value), the velocity at this portion could not be quantified accurately. Therefore, the velocity quantification was carried out at portions with absolute cosine values larger than 0.05 (Doppler angle outside the range of 87.134 to 92.866 degrees). As a result, the flow phantom was divided into two different segments, as shown in the flow velocity distributions in Fig. 7. The flow velocities of the two segments before Doppler angle correction were 4.718 ± 1.014mm/s and 6.484 ± 1.358mm/s. After Doppler correction, these velocity values go up to −26.424 ± 2.835mm/s and 26.272 ± 2.391mm/s, accordingly. The corrected absolute flow speed was very close to the theoretical velocity of the pumping rate, and the deviation was within 10.8% of the measured velocity.

Fully distributed absolute BFV measurement for MCA
Fully distributed absolute BFV quantification for MCA was demonstrated by the proposed volumetric reconstruction method. Three-dimensional OCT data were acquired at A-line rates of 50 kHz, 25 kHz, 10 kHz, 5 kHz, and 2.5 kHz resulting in time intervals of 0.02 ms, 0.04 ms, 0.1 ms, 0.2 ms and 0.4 ms, respectively. According to Eqs. (1)-(3), the axial velocity ranges for PRCD were ± 11.70 mm/s, ± 5.85 mm/s, ± 2.34 mm/s, ± 1.17 mm/s and ± 0.58 mm/s for ΔT = 0.02 ms, 0.04 ms, 0.1 ms, 0.2 ms, and 0.4 ms, respectively. Also, the tilted angle of the output laser beam was carefully tuned to avoid the phase wrapping of the blood flow. The imaging area for subject 1 was 8.4 mm by 8.4 mm while for subjects 2 and 3, it was 4.2 mm by 4.2 mm. These data were processed later to obtain the three-dimensional OCT structural and DOCT images, the dimensions of which were 512x1024x256. The DOCT images were then fed into our algorithm to verify the algorithm's applicability. For better reference of the vessel distribution, the intensity-based Doppler variance (IBDV) [ The IBDV method uses only the intensity information such that it is suitable for phase instable OCT systems. J = 4, N = 4 were used for the calculation of IBDV in our experiments.
The far left column of Fig. 8 shows the en face maximum intensity projections (MIP) of the IBDV images for all three subjects. The white arrows indicate the branches of the MCA.
After vessel segmentation, MCA structures were reconstructed (2nd column of Fig. 8) and further smoothed by averaging neighboring pixel coordinates to produce the refined vessel lumen structures (3rd column of Fig. 8). The MCA skeletons were obtained afterward as displayed in the far right column of Fig. 8 with artery branches coded in different colors. The reconstructed vessel structures coincide with the en face image and reveal the necessary vessel geometry for absolute BFV measurement, and the overlapping of different vessels and the MCAs was successfully avoided. Doppler angles on the entire reconstructed MCA skeletons were calculated using the proposed averaged gradient method. Figure 9(a) shows the angle measurement result for subject 1 with the cosine value of the Doppler angle at each position coded in red and black. It can be seen that the cosine value was in good accordance with the vessel geometry where the red color represents a steeper vessel, and the black color represents the Doppler angle of 90°. Figures 9(b) and 9(c) show the enlarged details of the dashed-line-frames in Fig. 9(a) with arrows showing the flow direction at each discrete point which represent the vessel center at each B-scan image. The flow directions are correctly identified using our method.  Once the Doppler angles were derived, the absolute BFV along all the MCA was found. Among all these vessel branches, the images captured at 50 kHz (ΔT = 0.02ms) were used to calculate the absolute BFV, resulting in an axial velocity range of ± 11.70mm/s for the PRCD method. Figure 10 shows the BFV calculation results for all three subjects obtained before and after Doppler angle correction. In the corrected BFV distributions (2nd column, Fig. 10), non-uniformity was clearly seen at the corrected absolute BFV distribution (indicated by the changing color of the vessels). This velocity fluctuation was believed to be caused by the natural cardiac cycle of the animal. Since the OCT was operating at a point-by-point scanning scheme, each sampling was acquired at slightly different time points. Therefore, when scanning a pulsatile flow, the changes in flow velocity were captured and resulted in the discontinuity of BFV distribution. Also, in the uncorrected BFV distributions (1st column, Fig. 10), the blood flow signals were too weak to measure at the bottom of the main branches (limited by the penetration depth of the OCT system) and at the tips of the side branches (due to Doppler angles near 90°). Thus, the absolute BFV shown in the second column of Fig. 10 was extended in these sections using the mean branch BFV values. Among all the artery branches we measured, the overall mean absolute BFVs of 15.27 mm/s for subject 1, 13.65 mm/s for subject 2, and 17.75 mm/s for subject 3 were found, and the overall ratio between corrected and uncorrected BFV was 1.7249-fold, 2.9148-fold and 3.0590-fold for these subjects, respectively. Also note that the ratios between the corrected and uncorrected BFVs varied across different branches which was caused by the varying vessel geometry (i.e., different cosine value at different vessels/positions). Moreover, a failure case of our algorithm can be found in the second bifurcation of subject 3 (the bottom graph in the 2nd column, Fig.  10) where very low BFV was obtained. This failure was due to the loss of vessel signal caused by the 90 degree Doppler angle. However, the flow signal in those segments could be revealed at a longer time interval. The total average fluxes of the MCA blood flow in different branches were further calculated for validation. We first measured the average diameters for all the branches of the MCAs in different subjects, and then calculated the total average flux for these branches with the diameters and absolute BFVs using the following equation: where v is the average absolute BFV, and D is the average diameter of the vessel branch.
The results are listed in Table 1. The main branches have larger diameters than the side branches, and the flux in the main branches coincides with the total flux in the side branches across different subjects. The maximum deviation of the flux in the main and side branches is 0.04-mm 3 /s. This error was believed to be introduced by the error of the angle calculation and the ambiguity of scanning a pulsatile flow.

Quantification of MCA absolute BFV changes reveals the consequences of permanent MCA occlusion
Stroke is one of the leading causes of death and long-term disability worldwide. Ischemia, with cell death resulting from decreased blood supply to a brain area, characterizes approximately 90% of all strokes [40,41]. Here, we propose to use DOCT to determine absolute blood flow velocity changes during permanent MCA occlusion. Ischemic conditions were achieved via surgical occlusion and transection of the M1 segment (just distal to lenticulostriate branching) of the left middle cerebral artery such that only MCA cortical branches were affected. The permanent MCA occlusion surgery follows the same protocol as in [31]. A double ligature was tied and tightened around the MCA, and the vessel was then transected (completely severed) between the two knots. Three-dimensional OCT data were obtained before (baseline) and right after the occlusion surgery, and the DOCT images were produced afterward. Figure 11 displays the en face standard deviation projection of the acquired IBDV images and the en face MIP projection of the Phase-Resolved Color Doppler (PRCD) images for baseline and post-occlusions at different time intervals. A reduction of the flow velocity was found in the MCA immediately after the occlusion. The velocity in all the MCA branches dropped to an undetectable level until acquired with an effective 5K A-line rate. Using the proposed Doppler angle correction algorithm, the absolute BFV was measured in different MCA branches. Figures 11(a)-11(f) display the absolute BFV changes in six different positions indicated in the PRCD MIP images, respectively. The velocity in the main MCA branches (a) and (b) suffered from sharp decreases of 93.11% and 82.37% while the side branches (c-f) witnessed larger BFV reductions at 94.17%, 91.49%, 95.32%, and 93.88% of the original flow velocity. Notably, in one of the MCA branches (c), blood existed with greatly reduced but opposite flow direction post-occlusion. The BFV fluctuations in the side branches (c-f) were further monitored by repetitively sampling the B-scans 100 times. It was found that the BFV ranged from −23.04% to + 33.97% across all four branches at baseline but was reduced to −19.7% to + 30.88% after the occlusion. The proposed absolute BFV measurement approach shows good potential for further investigation of hemodynamic changes during permanent MCA occlusion.

Discussion
The proposed absolute blood flow velocity measurement method has the following significance: 1) the distribution of the absolute blood flow velocity along the entire MCA (including main-branches and sub-branches) can be reconstructed. This provides a new way for in vivo monitoring of brain hemodynamics, especially the blood flow changes in stroke and other neural injuries. 2) The measurement principle is based purely on post-computation of the reconstructed DOCT signal. No additional hardware is involved for absolute BFV quantification, and the algorithm can be implemented to any other OCT system for brain imaging.
3) The one-click automated procedure enables faster blood flow data calculation and removes observer deviations in manual operations. In fact, our algorithm works not only for MCAs, but also other vessels with relatively large diameters.
Nevertheless, subject to the limited dynamic range of the Doppler phase shift measurement, the DOCT images at different time intervals need to be acquired, which leads to a relatively long acquisition time (~3 mins for 5 intervals used in our experiments). Also, the successful performance of the algorithm depends on the signal-to-noise ratio (SNR) of the original DOCT images. The segmentation of the vasculature network will be greatly enhanced when the system resolution and the SNR are improved.

Conclusion
We demonstrated the fully distributed measurement of absolute blood flow velocity for the middle cerebral artery by reconstructing the vessel geometry using pure DOCT stack images. Given a seed point location in the MCA images, the algorithm automatically localizes the vessel lumens, calculates the Doppler angle, and then quantifies the absolute BFV on each position of the whole MCA, including the main vessels and side branches. Our approach is verified by quantifying the hemodynamics of cerebral arteries during permanent MCA occlusion with custom-built, swept-source optical coherence tomography.