Automated segmentation and quantification of OCT angiography for tracking angiogenesis progression.

Angiogenesis is recognized as a crucial component of many neurovascular diseases such as stroke, carcinogenesis, and neurotoxicity of abused drug. The ability to track angiogenesis will facilitate a better understanding of disease progression and assessment of therapeutical effects. Optical coherence angiography (OCTA) is a promising tool to assess 3D microvascular networks due to its micron-level resolution, high sensitivity, and relatively large field of view. However, quantitative OCTA image analysis for characterization of microvascular network changes, including accurately tracking the progression of angiogenesis, remains a challenge. In this paper, we proposed an angiogenesis tracking algorithm which combines improved vessel segmentation and brain boundary detection methods to significantly enhance time-lapse OCTA images for quantification of microvascular network changes. Specifically, top-hat enhancement and optimally oriented flux (OOF) algorithms facilitated accurate segmentation of cerebrovascular networks (including capillaries); graph-search based brain boundary detection enabled coregistration of 3D OCTA data sets from different time points for accurate vessel density assessment and analysis of their changes in various cortical layers. Results show that this algorithm significantly enhanced the accuracy of vessel segmentation compared to Hessian method. Application to chronic cocaine intoxication study shows effectively reduced errors in chronic tracking of microvasculature and more accurate assessment of vessel density changes induced by angiogenesis.


Introduction
Optical Coherence Angiography (OCTA) is a non-invasive and label-free method to obtain high contrast 3D image of blood vessel network in vivo utilizing the dynamic OCT speckle signal induced by moving red blood cells [1]. For example, B. J. Vakoc, etc, applied optical frequency domain imaging to monitor tumor vascular change in response to antiangiogenesis therapy [2]. Vivek J.Srinivasan etc used OCT to study the progress in experimental stroke by measuring capillary perfusion, vessel diameter and cerebral blood flow volume [3]. Zhongdi Chu etc detected and assessed retinal vascular abnormalities with multiple quantitative indexes obtained from OCTA [4]. J. You, et al., showed microvascular adaptation in mouse cortical brain resulting from chronic cocaine exposures [5]. Different image acquisition and processing methods have been introduced in these studies tools to identify changes in the microvascular network caused by neurovascular diseases; however, more accurate quantitative assessment of OCTA images is needed to better understand angiogenesisassociated disease mechanisms and to develop more efficient therapeutic strategies.
Two major issues remain for OCTA to track chronic angiogenesis. Firstly, the accuracy of widely used quantitative tools such as vessel area density and vessel skeleton density [2, 4, 6] rely on the quality of the segmented binary vessel image. Several approaches for automatic vessel segmentation on OCTA images have been reported, among which the simplest method is to apply a local adaptive threshold to obtain vasculature with a higher intensity than background tissue [7]. However, this method suffers from its tolerance to background noise and sensitivity to minute vessels. Another method is to use a matched filter for vessel enhancement [8]. Although this method presents better performance on noise suppression, it is limited by the assumption that cross-sectional intensity profile of a vessel follows a Gaussian shape, which is not the case in OCTA images. Hessian-based methods were frequently applied to segmenting OCTA images [5, 9,10], which utilizes the ratio of eigenvalues of Hessian matrix to distinguish tubular structures (e.g., vessels) from other structures. However, visual inspection and quantitative evaluation show that Hessian method tends to generate vessel discontinuity due to its blurring effect on small vessels. Overall, these methods are able to effectively segment vessels with high SNR, but suffer from compromised sensitivity for small vessels with a low image contrast. Considering the fact that angiogenesis leads primarily to capillary density increase [11], the errors for small vessel segmentation may likely result in a biased assessment of angiogenesis. Secondly, the comparative study must be made on identical tissue region to accurately assess microvascular network changes over time. However, common rendering 3D software fails to select identical volume from each data set because of irregular cortex surface and variable mounting angle. This limitation can lead to an inconsistent pattern of angiogenesis progression.
To address these two issues, we report an angiogenesis progression tracking algorithm, which combines a new vessel segmentation method and an automatic brain boundary detection method. The former enhances the performance on segmenting microvasculature in OCTA images based on top-hat enhancement [12] and optimally oriented flux (OOF) [13].The latter detects brain boundary by graph search after applying dual filtering to remove reflection artifacts, by which changes in microvasculature and flows in different cortical layers can be analyzed and characterized. We compare the method with Hessian for improving the accuracy to quantify microvascular density and demonstrate the efficacy for tracking chronic cocaine induced angiogenesis in different cortical layers of mouse brain.

Animal preparation
Mice (C57/B6, male, 12-14 weeks of age, Jackson Lab) were anesthetized using 2.0-2.5% inhalational isoflurane. Then, a rectangular cranial window (roughly 2.5 × 2.5mm 2 ) was created above the frontal cortex. The exposed cortical surface was covered by 2% agarose gel and affixed with a 100μm-thick glass coverslip using biocompatible cyanocrylic glue. Dental cement was used on edges of the window to secure its attachment with skull. For chronic cocaine study, mice received intraperitoneal (i.p.) injections of cocaine (30mg/kg/day/each) for consecutive 28 days during which the animals underwent periodical µOCTA scans. All of the animal experiments were approved by the Institutional Animal Care and Use Committee of Stony Brook University.

System setup
A ultrahigh-resolution optical coherence tomography setup (μOCT) was used to acquire invivo structural (μOCT) and μOCTA images of the cerebral microcirculatory networks [14], in which an ultra-broadband light source(λ = 1310nm; Δλ FWHM ≈220nm) was used for illumination to achieve 2.5µm axial resolution. The beam exiting from sample arm is collimated to ф5mm, steered transversely by a galvo scanner (VM500, General Scanning) and then focused by an achromatic doublet (f16mm/NA0.25) onto mouse cortex. The lateral resolution is determined by the effective NA to be 5.17μm. The backscattered light from mouse cortex was recombined with reference light in the detection fiber connected to a custom spectrometer in which the interference fringes spectrally encoding the depth profile (A-scan) were detected by a fast linescan InGaAs CMOS camera (2048-pixels; GL2048, Sensors Unlimited) synchronized with sequential transverse scans for image acquisition. Meanwhile, GPU with custom GUI programming was implemented to enable real-time FFT for rendering and display of 2D/3D µOCT images. Specifically, 4 consecutive B-scans (x-z plane) at each cross section (y-axis) were acquired at 10 fps to derive a cross-sectional μOCA (x,z) image as its relative standard deviation [15].

Algorithm overview
The overview of our algorithm is shown in Fig. 1, which outlines the two-step architecture of our algorithm design and development. First, the brain boundary was detected on OCT intensity B-scan by a graph search method after cover glass reflection artifacts removed. The detected brain boundary was applied on OCTA B-scan to segment brain volume. The volume was split into 4 layers and OCTA maximum intensity projection image (MIP) was generated from each layer. Then the vessels including capillaries from the OCTA MIP images were identified by the proposed vessel segmentation method combining top-hat enhancement and OOF. Enhanced OCTA MIP images and density map were obtained by using binary vessel mask. Our algorithm was implemented with custom software written in Mablab 2016a (Mathworks, Natick, MA).

Vessel segmentation
To accurately segment microvasculature, we developed an algorithm based on top-hat filter and optimally oriented flux (OOF) for small vessel detection. Figure 2 shows the schematic diagram for vessel segmentation. Top-hat filtering is for preprocessing to reduce noise and enhance linear structure, in which an opening of reconstruction operation is first applied to the original image using linear structure element of length l at different angles. Then, multiple top-hat operations with the same structure element are utilized until a maximum response is yielded, i.e., contrast enhanced image. Here, length l is selected according to capillary diameters (e.g., 5-8µm) to improve the contrast of microvasculature ( Fig. 2(d)). After preprocessing with top-hat for contrast enhancement, OOF is applied for vessel segmentation (Fig. 2(e)). OOF detects vessel structure by measuring oriented gradient flux which is the amount of gradient projected along directionρ flowing in or out a local 2D circle or 3D sphere. If the projected directionρ matches the cross-sectional direction, the oriented gradient flux achieves the maximum value. The oriented gradient flux for a 2D image is given by [13]: where G is Gaussian function with a scale factor of 1 used to obtain image gradient of I , ρ is the projected direction of outward gradient flux, r is the radius of circle, ĥ nr = is the relative position vector, r S ∂ is the boundary of the circle. ( , ; ) f x r ρ can be written in quadric form of matrix , , is the entry of matrix at the ith row and jth column. OOF detects vessel structure by finding the optimal projection direction that maximizes ( , ; ) f x r ρ , a procedure to solve eigen decomposition problem on matrix , r x Q from which two eigenvalues can be computed and ordered as 1 | | λ > 2 λ . To eliminate false response on edge of OOF, a gradient antisymmetric function, oriented flux antisymmetry vector, was used to measure edge response [16]: Then, the antisymmetry gradient along vessel cross-sectional direction can be calculated as: where 1 v is the eigenvector corresponding to 1 λ . Using this function to compensate the false positive response on edge, the final vesselness is computed as: is the maximum value of 1 λ over scales r . However, it was found that OOF's response to large-scale objects was attenuated by vessel's irregular boundary. To solve this problem, we utilize OOF just for detection of microvasculature by selecting a small maximum scale value. Since large (e.g., branch) vessels have high contrast, they can be extracted by binarizing the image with a high global threshold. The threshold method does not take vessel shape into account, thus an anisotropic diffusion filtering is performed to enhance boundaries of large branch vessels before image binarization ( Fig. 2(b)-2(c)) [17].
With both microvasculature and large vessel masks derived, the final segmented image is generated by merging (i.e., 'or' operation) these two masks ( Fig. 2(g)).

Automatic segmentation of brain surface
To accurately evaluate the progression of angiogenesis, there is a need to compare the vascular density change of identical brain tissue over time. However, manual segmentation using commercial 3D image rendering software may fail to extract identical brain tissue from time-lapse image cubes (e.g., over days and even weeks) due to potential variations in head mounting angle and uneven brain geometrical shapes. To enable extraction of identical brain tissue from the acquired images at each time point, we implemented an automatic brain boundary detection algorithm. Direct boundary detection may fail due to strong reflection from the cover glass of cranial window as highlighted by the dash box in Fig. 3(a). Thus we first use Laplacian of Gaussian (LoG) filter in combination with first order of Gaussian filter to remove the artifacts. The cross-section of glass artifacts often has a symmetrical profile, so LoG filtering gives high response to a glass surface. However, as shown in Fig. 3(c), it yields strong response not only to a glass interface but also to a sharp edge of tissue, which leads to many false detections if the output image is directly binarized. Thus first order of Gaussian filtering (FoG) is used to suppress the edge response of LoG. Figure 3(b) shows that despite the high turnouts at the edge of tissue and glass, that at glass center is close to zero or eliminated. This allows us to identify glass by selecting pixels with high LoG output but low FoG output. The thresholding matrix used to binarize LoG response can be constructed as: where Fr and Lr are responses of 1st-order Gaussian and LoG filters, respectively, and t is a constant. The final mask of cover glass is shown in Fig. 3(d), which can be applied to the original OCT image to remove glass artifacts. After artifacts from glass surfaces are removed, the brain boundary is delineated by a graph search method, which treats pixels in an image as graph nodes and assigns weights to edges between them [18][19][20]. We first define the image gradient using Sobel operator: where Sv and Sh are the response of vertical and horizontal Sobel operator and w is a weight to enhance the horizontal structures. Then, a normalized gradient graph is established as: where C is a normalized value, which is low for brain boundary. Here boundary detection is equivalent to search for the path with a minimal overall graph cost ( , ) t i j defined as: This shortest path problem can be easily solved by dynamic programming algorithm [21]. Figure 3(e) shows boundary detection result before artifacts removed. Due to attraction of reflection artifacts, graph search method fails to localize the brain boundary. In comparison, boundary is successfully detected after removing artifacts in the preprocess (Fig. 3(f)).

OCTA image enhancement
To better visualize the microvascular network changes in angiogenesis progression, we applied the binarized vessel mask to enhance the original OCTA image. The contrast of original OCTA images is improved by contrast limited adaptive histogram equalization (CLAHE) [22]. Amplified background noise and boundary blurring induced by CLAHE are minimized by multiplying the enhanced image with the binary vessel mask. Figure 4(a-a') shows that the raw OCTA image (MIP) suffers from nonuniform signal and background noise, thus low-contrast microvasculatures are not clearly resolved. After enhancement with CLAHE and binary vessel mask, the OCTA image shows overall significantly improved signal-to-noise ratio and the boundaries of microvasculature network are sharp and well defined.

Vessel density quantification
To quantify brain microvascular density, the microvascular network is skeletonized from the segmented vessel mask using a skeleton algorithm in ImageJ [23]. In the evaluation of the segmentation algorithm, the vessel skeleton of the ground-truth mask was manually corrected to avoid the bias induced by skeletonization. An adaptive rolling window is applied on the skeleton map to allow for spatially resolved density measurement. The size of window is 27 × 27 pixels in our study. The vessel density within a local window is quantified by a fill factor (FF) as:

Total skeleton pixel Total pixel of window
Then, a spatial resolved density map can be obtained by moving the window in both vertical and horizontal directions.

Verification of results
Dice similarity coefficient (DSC) [24] is used to quantitatively evaluate the results from our vessel segmentation method, which is defined as: where TP, FP, and FN are true positive, false positive, and false negative rates. DSC has a restricted range of [0,1], where a larger value indicates higher similarity between two segmentation results. The ground-truth vessel masks used for quantitative evaluation were performed by two different specialists. We defined vascular density error to compare vascular density obtained from binary vessel masks of different segmentation methods, which is given as: where VD is vessel density obtained from masks of automatical methods and GT is the ground truth of vessel density obtained from manual segmentation mask. We compared our OOF method with Hessian method that has been frequently used to delineate vessel structures in OCTA images. Figure 5 compares the vessel segmentation results in which the upper panels show a raw MIP of 3D OCTA (a), manual segmentation as ground truth for comparison (b), automatic vessel segmentation using the hybrid method [25] combining intensity-based method and Hessian matrix with scales of 2≤σ≤10 (c), and the segmentation using our proposed method with the radius of local circle in OOF set to 1≤r≤5 (d). The lower panels (a'-d') are the corresponding zoom-in images to show pixel-wise comparison among these methods. False-positive (FP) and false-negative (FN) pixels were highlighted in orange and blue colors, respectively. The vessel skeletons generated from each segmentation were overlaid with vessel masks for a better evaluation. As described above, Hessian method has poor performance on segmenting microvascular of low contrast and microvascular adjacent to other large structure. As OOF measures the gradients of a local circle which only achieves the maximum value when the circle touches vessel boundary, it has less chance than Hessian matrix method to extend beyond the boundary of target vessel to the structure in the vicinity. As a result, OOF has better performance on detecting small vessels especially those adjacent to large structures (e.g., veins or arteries). In addition to avoiding inclusion of nearby objects, OOF can retain the sharpness of vessel edges because it does not rely on Gaussian kernel with different scales for multiscale detection. With additional contrast enhancement by top-hat preprocessing, it was found that OOF was able to accurately detect minute microvasculature (weak capillary flows) from raw images. As shown in Fig. 5(c'), Hessian method missed many small vessel segments in the microvascular network. In comparison, OOF was able to significantly improve the microvascular network continuity. For instance, 2 red arrows point to 2 capillary fragments which were missed by Hessian method but were fully restored by OOF in Fig.  5(d'). In addition, 3 green arrows point to 3 capillary fragments adjacent to large branch vessels missed by Hessian method but segmented by OOF (Fig. 5(d')), showing that OOF can effectively reduce false negative rate. Due to the close correlation between vessel mask and skeleton, OOF leads to improved accuracy for vessel skeletonization as shown in Fig. 5. In terms of capillary detection, a great number of skeletion fragments were observed as a result of errors of the Hessian mask (i.e., blue highlights in Fig. 5(c')); whereas the skeleton network derived from OOF vessel mask (Fig. 5(d')) matched well with that of the manual segmented mask (Fig. 5(b')).

Evaluation of vessel tracking
To quantitatively evaluate the segmentation performance, results from OOF and Hessian methods based 8 OCTA MIP images were compared using Dice similarity coefficient (DSC) as defined in Eq. (13). Figure 6(a) shows that the DSC of segmentation by OOF is 0.848 0.024 ± , which is close to ground truth (DSC = 1); whereas DSC by Hessian method (0.692 ± 0.063) is significantly lower (p<0.05, n = 8). This validation result is consistent with those shown in Fig. 5 that the new OOF method is able to segment more microvessels, resulting in a closer agreement with the ground truth. Additionally, Fig. 6(b) indicates that based on statistical analysis (n = 8 OCTA images), the error of vascular density of Hessian method (34.1% ± 6.1%) is significantly reduced to 3.6% ± 2.1% by OOF. This result shows that OOF is able to more accurately quantify vascular density for quantitative tracking of angiogenesis. It must be noted that both OOF and Hessian methods share a very similar computational frame-work. The major difference is the way to construct the matrix which encodes the gradient information of a vessel. The computational complexity of constructing Hessian matrix is O (Nlog(N)) if the image is convolved with a Gaussian kernel in Fourier space. The matrix of OOF in Eq. (3) can also be computed using fast Fourier Transform, and the overall computational complexity is the same O(Nlog(N)) [13].

Brain layer segmentation
After removing glass artifacts, the boundary was detected on an OCT intensity B-scan and then as a mask applied to the corresponding OCTA B-scan as shown in Fig. 7(a) and Fig.  7(b), respectively. With this process performed on all B-scans in 3D OCT data set, we divided 3D OCTA volume into 4 layers which roughly cover the range of cortical layers I to IV [26] (Fig. 7(d)). Compared with the original 3D volume (Fig. 7(c)), we observed that the artifacts above brain surface as indicated by yellow arrows are removed. To further demonstrate the effectiveness of the proposed method, we compared MIP images from tissue volume that were respectively acquired from a common rendering software and the proposed boundary detection method at different time points. MIP images by manual segmentation with commercial software show inconsistent shadows and patterns of vasculature at two time points (upper panel, Fig. 7(e)), likely resulting from irregular brain surface and inaccurate optical alignment (animal head mounting differences). However, with our boundary detection method, such inconsistency and artifacts are well compensated and the vascular patterns in 2 images are matched (lower panel, Fig. 7(e)), thus facilitating more solid quantitative comparisons. Fig. 7. Implementation of automatic brain segmentation. a)-b) Boundary was detected on OCT cross-sectional intensity image and then applied on OCTA image. c) Original 3D volume of OCTA images. Arrows are used to highlight noise before applying brain boundary detection. d) 3D volume of OCTA is segmented into four layers. e) MIP images at two time points from maunal segmentation with commercial software and the proposed boundary detection method.

Chronic cocaine elicited angiogenesis
Cocaine has been reported to induce acute decrease in cerebral blood flow (CBF) and cerebral blood volume [27,28]. In addition, chronic cocaine abuse can lead to cerebral hypoperfusion and thus hypoxic environment which may trigger the cascade of releasing angiogenetic agonists [29]. In this section, we utilized proposed segmentation method to track angiogenesis progression over cocaine treatment time on different layers of frontal cortex. The enhanced OCTA images introduced before were performed to better visualize the vascular network change before and after cocaine induced angiogenesis. To investigate chronic effect of cocaine, we longitudinally imaged mice frontal cortex every 7 days for consecutive 21 days. The Fig. 8(a) compares OCTA images of full thickness cortex (400µm) at day 1 and day 21. The zoom-in images corresponding to blue boxes clearly shows vessel density at day 21 has dramatic increase than that at day 1. With help of enhanced OCTA image, newly grown microvasculature can be easily identified as indicated by red rows. Figure 8(c) presents quantification of vessel density change of full thickness cortex during cocaine treatment. The vessel density from baseline (day 1) to 3 weeks after cocaine exposure increased by 12.5% ± 1.75%(m = 7, p<0.05). We further examined the vascular density change on different layers of frontal cortex. Vascular density was quantified on 4 different layers (100 m μ ) of frontal cortex over 3 weeks. The first column of density map array ( Fig. 8(b)) displays vascular density variation over depth at day 1 before cocaine-induced angiogenesis happened. We found vascular density(FF) increase from 0 to 400 m μ in which density of layer 4 (300-400 m μ ) increased by 23.1% compared to Layer 1 (0-100 m μ ) from 0.09 ± 0.009 to 0.1108 0.005 ± . This vessel density distribution over depth is consistent with Two-Photon Microscopy(TPM) and Optical Doppler Tomography(ODT) measurement reported previously [14,30]. The horizontal direction of density map array shows obvious density increase on all 4 layers, suggesting angiogenesis was ongoing both on surface and deep tissue due to cocaine treatment. It also can be observed from Fig. 8(d) that vascular density(FF) of layer 1 to layer 4 have similarly increased by ~0.014-0.02 after 3-week cocaine exposure. All layers present a linear vascular density increase over time, indicating that cocaine-induced angiogenesis is dependent on the time (e.g., accumulative dose) of chronic cocaine exposure.

Discussion and conclusion
Accurate quantification and detailed visualization of vascular network change in brain cortex are crucial for understanding neurovascular disease and progression. OCTA is a 3D imaging technique that can provide label-free, capillary-resolution angiograms based on decorrelation or variance of interferometric signals caused by moving red blood cells. Recent studies show that OCTA can be uniquely useful for identifying angiogenesis triggered by various diseases including cocaine toxicity and carcinogenesis. However, longitudinal tracking of angiogenesis progression remains technically challenging. Two problems are encountered that hinder unbiased tracking of angiogenesis: (1) time-to-time variation of selected data volume induced by brain geometric change or optical alignment, and (2) inaccurate vessel quantification caused by current vessel segmentation method. In this paper, we have proposed a new angiogenesis progression tracking algorithm to address the problems and validated the effectiveness of the method in studying cocaine induced angiogenesis. By employing a brain boundary detection method, our algorithm effectively eliminates the influence of brain geometric variation, which allows the data volume at different time point to be coregistered. Results in Fig. 7(e) demonstrate the efficacy of our brain boundary detection method for improving the consistency of vasculatural patterns over time. Moreover, the boundary detection method can be used for cortical layer segmentation, facilitating quantitative analysis of angiogenesis progression in different cortical layers. Our new vessel segmentation method enables effective reduction of errors in vessel quantification by combining top-hat enhancement and OOF. Comparative results clearly show that this new method significantly improves small vessel segmentation in comparison to Hessian matrix method that has been widely used for OCTA studies. The error of vessel density was decreased from ~34% to ~3.6% with the new method ( Fig. 6(b)). For invivo validation studies, we showed the effect of chronic cocaine treatment on the vascular networks, including a drastic increase of microvascular density in various cortical layers.
As a limitation, our angiogenesis tracking method does not address the potential errors induced by the automatic skeletonization. The common drawbacks of automatic skeletonization such as spurious branch and missing skeleton may potentially bias the vessel density assessment. However, this type of error mainly occurs on large vessels and thus would not severely affect the tracking of angiogenesis which primarily leads to changes in capillary density. Additionally, the application of anisotropic diffusion filtering algorithm in our method can minimize the error by smoothing large vessel boundaries. In future studies, a more advanced skeletonization algorithm may be applied to further remove this type of errors [31]. Another drawback is that our vessel segmentation method identifies microvessels and large branch vessels separately due to attenuated responses of OOF at very large scales. This increases the number of parameters needed for manual sets dependent on the quality of the original OCTA data set. The introduction of anisotropic diffusion filtering might reduce the computational efficiency of our method. Future work may include simplifying the algorithm and improving computational speed. It should also be noted that our analysis of vessel density in different cortical layers may neglect the possible false positive errors induced by artifacts such as shadowgraphic projection from superficial large vessels on deeper layers. This problem can be mitigated by using a simple morphology operator to remove isolated pieces on vessel mask of deep layers. Another potential solution is to integrate advanced algorithms to reduce shadowing artifacts such as step-down exponential filter [2] and projection-resolved method [32]. Other improvements in the algorithm that can be made in future studies may include correction of refraction distortion [33], extension of limited depth of focus (DOF) [34] and application of a more reliable brain boundary detection algorithm [35].
In conclusion, we report an angiogenesis tracking approach based on top-hat enhancement and optimally oriented flux (OOF) algorithms to facilitate accurate segmentation of microvascular networks (including capillaries). Further incorporation of graph-search based brain boundary detection enables coregistration of 3D OCTA data sets from different time points for accurate vessel density assessment and quantification of their changes in various cortical layers. Application of this algorithm to in vivo longitudinal OCTA data sets of mouse cortical brain during chronic cocaine exposures demonstrates the utility of this method for accurate tracking of microvascular changes, which is clinically relevant for enhancing our understanding of disease progression and thus noninvasive assessment of therapeutic effects.