Method for quantitative assessment of retinal vessel tortuosity in optical coherence tomography angiography applied to sickle cell retinopathy

: Tortuosity is an important geometric vessel parameter and among the first microvascular alterations observed in various retinopathies. In the current study, a quantitative vessel tortuosity index (VTI) based on a combination of local and global centerline features is presented. Performance of VTI and previously established tortuosity indices were compared against human observers’ evaluation of tortuosity. An image-processing pipeline was developed for application of VTI in retinal vessels imaged by optical coherence tomography angiography (OCTA) in perifoveal (6 mm × 6 mm) and parafoveal (3 mm × 3 mm) regions centered on the fovea. Forty-one subjects (12 healthy control (NC) and 29 sickle cell retinopathy (SCR)) and 10 subjects (5 NC and 5 SCR) were imaged in the perifoveal and parafoveal regions, respectively. The relationship between VTI and age was examined in the perifoveal regions in NC subjects. VTI was measured from the OCTA images and compared between NC and SCR subjects using generalized least square regression with and without adjusting for age and race. VTI was found to correlate better than the 4 previous indices with performance of human observers. In the perifoveal region, a significant correlation was observed between VTI and age (r = − 0.4, P<0.001, N = 12). VTI was higher in SCR than NC subjects in perifoveal and parafoveal regions (P ≤ 0.001). The results demonstrate that the proposed method shows promise for detection of increased tortuosity in


Introduction
Sickle cell disease (SCD) is a genetic life-long chronic illness [1], characterized by formation of sickle shape erythrocytes [2]. Sickle cell retinopathy (SCR) is a major complication of SCD which is known to affect peripheral and macular retinal vascular beds and can lead to vision loss in the progressed stages due to formation of new vessels [3,4]. Furthermore, frequent presence of retinal neovascularization, hemorrhages, schisis cavities, and black sunbursts have been reported to be associated with SCR [3,4].
Fluorescein angiography (FA) has been used as a gold standard in the past decades for clinical evaluation of microvasculopathy due to SCR [5]. However, FA has been mainly performed for advanced proliferative stages of SCR and is known to be non-quantitative and invasive with limited applicability for subjects with small pupil [6]. Previous studies using spectral domain optical coherence tomography (SD-OCT) in SCR showed temporal and central macular thinning [7,8], and indeed macular thinning was found to be associated with proliferative SCR [8]. Lately, optical coherence tomography angiography (OCTA) became available and offered an alternative to FA by providing high-resolution non-invasive imaging of capillary network in different retinal layers [9][10][11]. Recent studies in SCR using OCTA in parafoveal regions showed microvasculopathies such as decreased superficial and deep retinal capillary density [5,12,13], enlarged foveal avascular zone (FAZ) [5,13,14], capillary nonperfused areas [5,15] and disruption of the perifoveal anastomotic capillary arcade [5].
Tortuosity is an important geometric vessel parameter which has been considered as a risk factor of multiple retinal pathologies [16], and increased vessel tortuosity is known to be among the first microvascular alterations due to many retinopathies such as diabetic retinopathy [17]. In SCR, increased retinal vessels tortuosity has been observed, but not considered as a pathognomonic sign. Non-specific increased tortuosity has been reported in more than 32% of SCR subjects by qualitative analysis of large retinal vessels [18]. Furthermore, increased tortuosity of major retinal vessels has been commonly observed in a majority of SCR subjects with HbS/HbS (SS) genotype [3,12]. Finally, a recent analysis of OCTA images of parafoveal regions showed vessel tortuosity is more sensitive than thickness for SCR detection [13].
There is no widely accepted mathematical definition for tortuosity [19], and tortuosity evaluation has been mainly performed subjectively by clinicians. However, sensitivity and repeatability of visual tortuosity evaluation is limited due to high inter-observer variability [20], and indeed qualitative evaluation is inefficient for large studies. Quantitative assessment of tortuosity alterations has been reported in various retinal pathologies [17,[21][22][23], using one or combination of distance measure (DM) (i.e. the ratio of vessel length to its chord length) [13,17,21,23], curvature or integral of curvature [16,22,24,25], and number of inflection points [26,27]. However, these methods are limited since they cannot always accurately estimate vessel tortuosity [13, 16, 17, 21-24, 26, 28-30], or are scale-dependent [25,27]. More recently, a measure of tortuosity based on slope chain coding (SCC) was proposed which is invariant to rigid transformations [31]. However, as suggested in [32], SCC depends on length of linear elements which needs to be carefully selected for accurate tortuosity quantification.
In the current study, a method for assessment of retinal vessel touristy is proposed and shown to be invariant to rigid transformations and perform better than previous tortuosity metrics when compared with the evaluation of human observers. VTI does not require any parameter tuning, is computationally efficient, and more sensitive to changes in vessel curvature than methods that rely only on DM. Additionally, OCTA images in both parafoveal and perifoveal regions were analyzed for detection of alterations in retinal vessel tortuosity due to SCR.

Subjects
The study was approved by an institutional review board at University of Illinois at Chicago. Informed consents were obtained from subjects in accordance to the tenets of Declaration of Helsinki. The study was performed in 2 cohorts of subjects based on size of the imaged retinal region (i.e. 6 mm × 6 mm (perifoveal); 3 mm × 3 mm (parafoveal)). Perifoveal comprised 41 subjects (13 males and 28 females) with ages ranging from 15 to 62 years. Parafoveal comprised 10 subjects (5 males and 5 females) with ages ranging from 27 to 56 years. Based on a complete clinical history and ocular examination, subjects were categorized into 2 groups of healthy control (NC, N = 12; 7 OD and 5 OS) or sickle cell retinopathy (SCR, N = 29; 15 OD and 14 OS) in perifoveal region, and NC (N = 5; 1 OD and 4 OS) or SCR (N = 5; 3 OD and 2 OS) in parafoveal region. Three NC and 2 SCR subjects were in both cohorts. Subjects demographics and clinical data has been shown in Table 2. Exclusion criteria were inability to give informed consent or participate in the study, diabetes mellitus, glaucoma, or any other retinal disease. Each subject contributed data to the study by one eye with the best image quality.

Image acquisition
Imaging was performed by a commercially available OCTA instrument (Optovue Inc, Fremont, California, USA). The laser wavelength was 840 ± 45 nm with an axial scan rate and axial scan depth resolution of 70 KHz and 5 µm, respectively. B scans were acquired from identical retinal locations to generate blood flow map based on motion of red blood cells. Each B-scan was comprised of 304 A scans. Images of the superficial retinal vessels and capillary network were generated in perifoveal (6 mm × 6 mm) and parafoveal (3 mm × 3 mm) regions centered on the fovea. The superficial layer was defined by the Optovue software and displayed the retinal vasculature in the nerve fiber and ganglion cell layers with minimal flow projection and shading effects.

Image processing and analysis
Assessment of superficial retinal vessel tortuosity was performed using several imageprocessing steps as shown in Fig. 1. The algorithms were developed in MATLAB (Release 2015b, MathWorks, Inc., Natick, MA, USA) with image processing toolbox version 9.0.

Vessel segmentation
Segmentation of vessels in the perifoveal and the perifoveal regions in OCTA was performed using a k-means clustering algorithm [33], similar to a previous study in which human observers performed the classification based on Horton-Strahler schemes [34]. K-means clustering is an iterative classification algorithm that divides data into k clusters where k is a positive integer. Clusters were found to minimize the least square error as shown in Eq. (1).
where E is the error, x is intensity value of pixels, and µ is the centroid of clusters. The algorithm starts with K randomly selected initial centroids. The distance between intensity value of each pixel to the centroids was computed to assign the pixel to a cluster with the closest centroid. A new set of clusters were computed in the next iteration based on assignment of pixels to the clusters in the previous one. This process was repeated to the point that the centroids were not changed between 2 consecutive iterations. The k-means clustering algorithm was used in the current study with k equals to 2 for classifying pixels either as vessel or background pixels to generate a binary image. The binary image obtained from the k-means clustering was enhanced by a series of morphological operations. Filling was performed to eliminate any hole within the vessels, thickening was performed to bridge between separated vessel branches, objects smaller than 100 pixels, computed by counting the number of pixels in each binary object, were removed to eliminate noise and capillaries. Finally, dilation was performed using a disk shape structuring element to smooth vessel walls. Vessel centerlines between each pair of bifurcation points were extracted using distance transform through manual endpoint selection. The vessel endpoints were selected by simultaneous visualization of the grayscale and binary images in 2 separate windows. The grayscale image was used to identify bifurcations and the binary image was used to more accurately locate them with respect to the vessel boundaries and centerline. Cubic smoothing spline was utilized with a regularization parameter of 0.01 to obtain an adequate centerline and avoid aliasing. Vessel tortuosity index (VTI) was then computed for each of the centerlines as described in the following section.

Vessel tortuosity assessment
The most frequent tortuosity measures found in the literature were based on DM which is the ratio of arch length to chord length of a vessel segment [17,[21][22][23]. Despite the fact that DM is a simple and quick approach for vessel tortuosity quantification, there are circumstances where it fails to accurately determine the tortuosity due to its simplicity [29]. A tortuosity measure by combination of DM and number of inflection points was suggested by Bulliet et al. [26], which was shown not always match with visual perception of tortuosity [27]. Additionally, tortuosity density index (DT) defined as multiplication of number of inflection points with sum of DM determined between points of changes in centerline curvature was proposed by Grisan et al. [27]. The main drawback of this method was scale dependency.
Tortuosity measures based on curvature or integral of curvature have been proposed and applied in various pathologies [16,24,30]. Hart et al. showed that the sum of squared curvatures along retinal vessel centerline perfectly match with tortuosity perception of experts for classifying their data set into tortuous and non-tortuous blood vessels [16]. However, these methods can lead to misrepresentation of vessel tortuosity without considering changes in the sign of the curve, which is an important parameter used by clinicians for tortuosity evaluation [27].
Tortuosity assessment based on local angle changes was proposed previously [35][36][37]. Grisan et al. showed limitation of these methods which computed the same tortuosity for a semi-circumference and a curve obtained by juxtaposition of its two arcs. The two have the same mean angle changes but different tortuosity.
A method based on SD of distribution of vessel centerline incremental lateral displacement was proposed earlier by Wenn et al. [38]. Nevertheless, their method does not consider the magnitude of the curve, which plays an important role in clinical evaluation of vessel tortuosity.
Recently, Lisowaka et al. compared performance of 5 different retinal vessel tortuosity measures against sampling rates of vessel centerlines on a public data set [32]. Performance comparison of DM [16], DT [27], SCC [32], and 2 curvature based measures [16] showed that the overall performance of DT was good, but not always the best. They suggested that attention to numerical details and standardization is essential before choosing a tortuosity index. Further details on available tortuosity measures is beyond the scope of the current study and can be found elsewhere [27,29].
Comparatively little work has been conducted for quantitative assessment of retinal vessel tortuosity in OCTA. The VTI presented in the current study is sensitive to small changes in tortuosity, and hence is suitable for detecting tortuosity alterations in retinal vessels in OCTA. For each point along the centerline, the angle (θ) between a line tangent to the centerline and a reference axis was determined. SD of absolute θ values (SD θ ) along each centerline was computed, representing variation of local angle changes. The number of critical points at which the first derivative of centerlines vanishes (N) was quantified for each centerline based on frequency of changes in sign of the slope of the tangent lines. The N value was set to 1 for centerlines with no critical point to avoid zero tortuosity due to the lack of critical points. Inflection points at which the second derivative of centerlines vanishes (Ip) along the centerline were identified by detecting changes in the sign of curvature (k) where k was calculated using Eq.
The tortuosity index was multiplied by vessel length (L A ) since L A increases with tortuosity, and was normalized by vessel chord length (L c ) to allow comparison of vessel segments with variable chord length. The mathematical derivation and visual demonstration of VTI is given in Eq. VTI is unitless similar to previous tortuosity metrics [16,27,32]. The minimum value for VTI is equal to 0 for an ideal straight line with zero SD of local angle changes. In theory, there is no maximum value for VTI since it increases with higher variation in angles, number of critical points, and the magnitude of curve. However, VTI in OCTA was generally lower than 1. NC and SCR example of perifoveal OCTA, vessel segmentation, and extracted centerlines for tortuosity analysis are shown in Fig. 3.

VTI validation
VTI was validated by (i) using sinusoidal curves with variable magnitudes and angular frequencies, and (ii) by quantitative comparison against performance of human observers (MS, WO and MK).
Two sets of sinusoidal curves were generated, the first set had constant angular frequency (1.25 rad/sec) and variable magnitudes between 0 and 6. The second set had constant magnitude of 1 and variable angular frequencies between 0 rad/sec and 6.3 rad/sec. A dependable tortuosity index should rise with increasing magnitude and increasing angular frequency coequal with visual perception of tortuosity.
Quantitative comparison of VTI and human observers' grading was performed using a set of 25 sinusoidal curves generated from randomly selected magnitudes between 0.1 and 5, and randomly selected angular frequencies between 0.6 rad/sec and 4 rad/sec. The random magnitudes and frequencies were generated using Mersenne Twister pseudo-random number generator [39]. VTI and 4 previously established tortuosity indices, namely DI, T nl [40], DM, and integral of absolute curvature (T c ) were computed for the curves and sorted in an ascending order. Similarly, human observers visually evaluated the curves and arranged them in an ascending order.

Statistical analysis
Statistical analysis was performed using SPSS (version 22, SPSS, Chicago, IL, USA) and statistical analysis codes written in MATLAB applied separately to data from the perifoveal and parafoveal regions. The association between VTI and performance of the human observers was determined using Spearman's rank correlation. Subjects' demographics were compared using Chi-square or t-test. Correlation between VTI and age in NC subjects was assessed using weighted Pearson's correlation analysis. Finally, generalized least squares (GLS) was used to determine the effect of disease (NC and SCR) on VTI with and without adjusting for covariates (age (continuous) and race (categorical)). Statistical significance was accepted at P≤0.05.

VTI validation
The validation tests using sinusoidal curves showed that VTI increased exponentially with increasing magnitudes and angular frequencies. The Spearman's rank correlation statistics for comparison of VTI and the 4 previous methods with performance of human observers is shown in Table 1. VTI was correlated better with all human observers' evaluations than tortuosity indices derived from previously established methods.

Demographic data
Subjects' demographics are summarized in Table 2. In perifoveal region, age and eye examined were not different among NC and SCR subjects (P≥0.3), while sex and race were different (P<0.03). In parafoveal, age, sex and eye examined were not different among NC and SCR subjects (P≥0.2), while race was different (P = 0.008).

VTI in perifoveal and parafoveal region
In the perifoveal region, VTI was assessed in 1026 and 2444 vessel segments in NC and SCR subjects, respectively. Mean VTI per subject ranged from 0.25 to 0.9 in perifoveal region. There was a negative correlation between VTI and age in NC subjects (r = −0.4, P<0.001, N = 12). Finally, VTI was significantly higher in SCR (0.61 ± 0.11) than NC (0.31 ± 0.04) subjects with or without age and race adjustment (P<0.001).
In the parafoveal region, VTI was assessed in 181 and 154 vessel segments in NC and SCR subjects, respectively. Mean VTI per subject ranged from 0.34 to 0.88 in parafoveal region. Finally, VTI was significantly higher in SCR (0.69 ± 0.18) than NC (0.40 ± 0.04) subjects with or without age and race adjustment (P≤0.001).

Discussion and conclusion
In the current study, a quantitative vessel tortuosity index was formulated by extracting multiple parameters from vessel centerline. This method is similar to human observer's evaluation because variation of local angle changes, number of critical points, and magnitude of curve, each contribute to visual perception of tortuosity and were included in VTI formula. The method was applied to OCTA images obtained in perifoveal and parafoveal retinal regions and tortuosity alterations in superficial retinal vessels due to SCR were demonstrated. VTI increased exponentially with higher magnitude and angular frequency of sinusoidal curves, consistent with visual perception of tortuosity. More importantly, VTI was shown to match with visual perception of tortuosity and its performance was better than other commonly used methods such as DI and DM. Furthermore, VTI has several advantages for assessment of retinal microvasculature. First, it is invariant to rigid transformations such as translation, rotation, and scaling because these transformations have no influence on any of its components (SD θ , N, M, L A /L c ). This is important since rigid transformations are common in clinical applications in which different instruments with variable setting are used for image acquisition. Second, VTI was normalized with respect to chord length which is necessary for detection of tortuosity alterations in retinal tissue which is densely vascularized and contains vessels with variable lengths. Third, VTI incorporates vessel length which can contribute to tortuosity. In fact, retinal vessels can become tortuous due to longitudinal stretching or shortening of distance between their tethering points [16,41]. Although the pathophysiology of vascular stretching is not well established, the assumption is that the tone of smooth vascular muscles which determines vessel length and consequently tortuosity can be altered by mediators, blood gas, and metabolism due to pathology [42].
The proposed and many previously established methods do not account for the influence of vessel caliber on tortuosity. Therefore, the findings may misrepresent tortuosity of the entire network since small caliber and highly tortuous vessels can increase the overall tortuosity of the network. Attempts have been made to consider vessel width information in retinal tortuosity quantification [25], and to detect type of tortuosity alteration in 3D models of intracerebral vessels [26]. However, large clinical studies are required to determine the effect of vessel size on clinical perception of tortuosity to provide a deterministic weight factor to aggregate the impact of vessel size on tortuosity of the network. Although the effect of vessel caliber was not considered in the VTI formula, the k-mean clustering algorithm used in the current study automatically eliminated capillaries with diameter smaller than 20 µm, and hence tortuosity was assessed in relatively larger size vessels within the network which are mainly evaluated by clinicians for detection of tortuosity alterations.
In the current study, increased retinal vessel tortuosity due to SCR was reported in perifoveal and parafoveal regions. This finding is in agreement with previous qualitative [3,12,18] and quantitative [13] reports of increased retinal vessels tortuosity in SCR. Despite the report of a moderate increase (17%) in tortuosity of parafoveal vessels [13], we found a marked increase in tortuosity in both the perifoveal (98%) and parafoveal (70%) regions. The difference in findings can be attributable to differences in techniques.
Retinal vessels tortuosity in the perifoveal region was negatively correlated with age in NC subjects, consistent with a previous study [43] that quantified tortuosity in major retinal arterioles and venules. This finding suggests that age has an independent effect on retinal vessel tortuosity and needs to be matched in future studies for detection of retinal tortuosity alterations.
There were some limitations in the current study. Accurate junction point (i.e. bifurcation and crossover) detection in retinal images is crucial for tortuosity evaluation since they define vessel course and their misdetection can result in inaccurate tortuosity measure. In the current study, a human observer (MK) visually inspected images and selected junction points. Despite FA and FP where bifurcation points can be more easily detected as shown in [44,45], presence of numerous junction points and a dense capillary network in OCTA increases complexity and decreases reliability of automatic junction point detection. Future research is warranted for developing methods for automatic junction point detection in OCTA which can improve the efficiency and increase repeatability of tortuosity assessment. Retinal arteries and veins were not separated since there is no reliable method to distinguish between them in OCTA. In fact, previous studies have shown that there is a difference in tortuosity of arteries and veins [43]. Nevertheless, tortuosity measurements were averaged per subject in the current study to reduce variability due to vessel type. Additionally, we assumed the parameters used in the VTI formula had equal weight toward tortuosity and indeed the strong correlation between VTI and initiative perception of tortuosity further confirmed this hypothesis. Nevertheless, large clinical studies are useful to determine the effect of each parameter on visual perception of tortuosity and make adjustment if necessary. In the current study, the feasibility of application of the method was demonstrated in a relatively uniform population of subjects, predominately with SS genotype and stage II retinopathy. Future studies in larger cohorts of subjects are needed to investigate the effect of genotype and stage of SCR on tortuosity alterations. Additionally, future studies would be useful to investigate the relation of blood vessel tortuosity alterations with other retinal vascular and anatomical abnormalities detected by multimodal imaging [46]. Assessment of retinal vessel tortuosity alterations shows promise for clinical diagnostic evaluation and longitudinal monitoring of microvasculopathies due to SCR. Furthermore, the method presented in the current study can be potentially applied to microvascular images of other tissues for quantitative assessment of vessel tortuosity.

Funding
NIH research grants DK104393, EY001792, Senior Scientific Investigator Award (MS) and departmental award from Research to Prevent Blindness.

Disclosures
The authors declare that there is no conflict of interest related to this article.