Enhancing Vessel Segment Extraction in Retinal Fundus Images Using Retinal Image Analysis and Six Sigma Process Capability Index

Abstract

Retinal vessel segmentation, skeletonization, and the generation of vessel segments are considered significant steps in any automated system for measuring the vessel biomarkers of several disease diagnoses. Most of the current tortuosity quantification methods rely on precise vascular segmentation and skeletonization of the retinal vessels. Additionally, the existence of a reference dataset for accurate vessel segment images is an essential need for implementing deep learning solutions and an automated system for measuring the vessel biomarkers of several disease diagnoses, especially for optimized quantification of vessel tortuosity or accurate measurement of AV-nicking. This study aimed to present an improved method for skeletonizing and extracting the retinal vessel segments from the 504 images in the AV classification dataset. The study utilized the Six Sigma process capability index, sigma level, and yield to measure the vessels’ tortuosity calculation improvement before and after optimizing the extracted vessels. As a result, the study showed that the sigma level for the vessel segment optimization improved from 2.7 to 4.39, the confirming yield improved from 88 percent to 99.77 percent, and the optimized vessel segments of the AV classification dataset retinal images are available in monochrome and colored formats.

1. Introduction

The retina, a thin and transparent tissue located at the back of the eye, is essential for detecting light and transmitting visual data to the brain. Recent technological advancements have developed numerous diagnostic instruments and computer-aided diagnostic systems (CAD) for ophthalmologists, allowing them to detect and monitor vascular morphological changes in the retina [1]. Fundus photography, a widely used technique, aids in identifying retinal alterations and monitoring the progression of eye diseases, resulting in earlier and more accurate diagnoses [2].

Retinal fundus images provide valuable information regarding patient health from the status of the vascular anatomy, including veins, arteries, macula, fovea, optic discs, and numerous aberrant lesions, including cotton wool spots, exudates, and hemorrhages. Fundus retinal image vascular skeletonization is an essential topic in medical image processing that has received much attention in recent years. The extraction of retinal vessel skeleton segments is essential for retinal image analysis and for diagnosing ocular diseases. It involves identifying and extracting the skeleton structure of retinal blood vessels, allowing for the quantification of vessel biomarkers and the evaluation of vessel tortuosity. Nevertheless, manual segmentation of the retinal vasculature from the colored fundus images is labor-intensive, time-consuming, and prone to inconsistencies between specialists [3,4]. Therefore, the demand for human-free automated retinal vascular segmentation has become crucial. Diagnosing retinal diseases such as diabetic retinopathy (DR) and hypertensive retinopathy (HR) is highly dependent on the study and analysis of blood vessel morphology changes [5,6].

Vessel segmentation, skeletonization, and the generation of vessel segments are significant steps in any automated system for measuring the vessel biomarkers of several disease diagnoses. For example, most current tortuosity quantification methods rely on precise vascular segmentation and skeletonization of the retinal vessels. Additionally, having a reference dataset that provides accurate vessel segmented images is crucial in developing a deep-learning solution and an automated system for measuring vessel biomarkers in various disease diagnoses, especially for optimized quantification of vessel tortuosity or accurate measurement of AV-nicking. A phenomenon noticed on the generated skeletons is that it segments the vessels at the cross-overs or branching points and within the vessel itself, thus creating many vessel segment cuts. Therefore, it implies that the generated skeletons need some enhancement or optimization to have the real vessel segments available only at the crossover or branching points.

Accurate segmentation of retinal blood vessels plays a crucial role in computer-assisted diagnosis and staging techniques, facilitating the detection of morphological changes resulting from various disorders. Unsupervised techniques, such as matching filters, vasculature tracing-based segmentation, and model-based segmentation, have been extensively studied for retinal vascular segmentation [7,8,9,10]. However, the absence of manually labeled ground truths frequently hinders the performance of unsupervised methods, making supervised approaches generally preferable. Unlike unsupervised techniques, supervised models can improve with access to standard gold labels and learn from annotated data. In supervised models, retinal vascular segmentation involves feature extraction and pixel classification. There is a dichotomy between autonomously learned and manually constructed characteristics. In machine learning, manual feature extraction is used for fundus images and common classifiers such as SVM [11], KNN [12], radial projection [13], and ridge-based schemes [14]. Manual feature selection is application-specific and may utilize domain knowledge, but it lacks generalization capability because it cannot autonomously learn new features [15]. Several other methods in [16] achieved vessel segmentation using shallow neural nets and deep learning approaches that showed promising results in multiple areas of fundus retinal image analysis and diagnostics problems. Deep learning algorithms for retinal vascular extraction were developed and tested on publicly accessible datasets for diabetic retinopathy and other retina diseases. These approaches perform better for retinal vessel segmentation because retinal image datasets contain more precise and apparent blood vessels than HR-specific datasets.

A seven-layer CNN was developed by Khalaf et al. [17]. To decrease intra-class variation, they classified pixels into big, background, and tiny vessels. In pre-processing, they isolated the green channels of pictures and used adaptive histogram equalization (AHE) and mathematical top-hat filtering. The green channel, AHE, and top-hat filtering strengthened training picture vasculature and reduced noise. FCNs generate dense and accurate picture patch pixel predictions with high performance for each of the seven FCNs.

Multi-model networks utilized by others, like Oliveira et al. [18], presented an FCN with skip connections in order to convey characteristics from shallow to deeper layers. SWT-transformation added additional channels as input to study the enabled multiscale structure of the vascular system, and it was concluded that the domain knowledge benefited the deep learning results.

Liu et al. [19] presented a skeletonization method using a multiscale vessel filter and adaptive thresholding. The method effectively captures vessel structures with varying widths. It exhibits robustness against noise but may encounter difficulties in accurately capturing vessel crossings. Wang et al. [20] provided a retinal image vessel crossing technique and a bifurcation detection technique. The authors used a multi-attention network for vessel segmentation and a directed graph search approach for vessel branching points. Although the suggested technique accurately detects vessel crossings and bifurcations, a critical feature of retinal image analysis, multi-attention networks improve vessel segmentation by concentrating on informative areas, improving detection performance. Exploring vascular-related issues using the directed graph search technique helps find branching spots. The work needs a detailed assessment of the suggested approach on big datasets, and a comparison with state-of-the-art methodologies is required to determine its efficacy and usefulness in real-world clinical contexts.

Despite the advancements in retinal vessel skeleton segment extraction, several challenges remain: Accurate extraction of vessel branches at crossovers and branching points remains a significant hurdle for many methods. Handling vessels with varying widths and irregular shapes poses difficulties for some approaches. Certain methods’ computational complexity and processing time may limit their practical application in real-time scenarios. The reliance on annotated datasets for training deep learning models may be a constraint in some contexts.

The main concerns addressed by this study are:

• Why does skeletonization of the vascular tree break the vessel segments in points, not at a bifurcation or cross-section locations?
• Why are vessel segments cut into several pieces between the connection points?
• How can other researchers be provided with a correctly segmented skeletonized large dataset by expanding the researchers’ previously released AV classification dataset with images that comprise properly segmented vascular segments?

Moreover, this study proposed an improved method for skeletonizing and extracting pertinent vessel segments. Our method enhances the tortuosity calculation process by optimizing the vessel segment extraction process. We used process capability analysis, derived from the six sigma methodology, to quantify the improvements in tortuosity calculation before and after implementing our optimized method. In addition, we expanded the existing AV classification dataset by incorporating image-level segment skeletons’ severity level.

The main contribution and novelty of this work are:

• Proposing an enhanced method to extract each vessel segment in the vessel’s skeleton tree from each intersection point to the next intersection point.
• The improvement of the vessel-segment extraction leads to enhanced vessel-tortuosity calculation.
• The tortuosity calculation improvements were quantified using the six-sigma process capability index, where the tortuosity was calculated for all the vessel segments in all the 504 retinal images twice, once before applying the enhanced method of vessel segments and the second time after enhancing the extracted vessel segments.
• For the first time in the field, this work used the process capability index to measure tortuosity improvement due to improving the process of vessel segment extraction.
• The approach was implemented and evaluated on a robust dataset of 504 retinal images.
• The introduction of a new extension dataset containing vessel segment fragments for the 504 images made available for researchers’ future work.

The remainder of this article is structured as follows: Section 2 reviews the relevant literature on skeletonization research. Section 3 describes our optimized vessel extraction technique. Results of the proposed optimized method are presented in Section 4. Section 5 concentrates on the influence of enhanced vessel segments on the calculation of tortuosity. Finally, Section 6 presents the conclusions.

2. Literature Review

Vessel segment extraction plays a vital role in various applications of retinal image analysis, including vessel tortuosity measurement, disease diagnosis, and monitoring. Numerous techniques have been proposed in recent years to address the challenges associated with accurate and reliable segment extraction. At the same time, significant progress has been made in this field in recent years due to the increasing availability of high-resolution retinal imaging and the evolution of sophisticated image processing techniques. This literature review provides an overview of state-of-the-art methods and their strengths and weaknesses. At the same time, it provides an overview of contemporary methodologies for extracting segments of the retinal vessel skeleton, emphasizing their essential characteristics and contributions.

Fundus retinal image vascular skeletonization is an active research subject, with numerous potential ways to increase the process’s accuracy and efficiency. The study of the retinal vessel-segments extraction needs multiple steps to reach the retinal vessels’ skeleton fragments. Hence, each step definition and related work are summarized below; therefore, to quantify the enhanced extraction of the vessel segments, the following vessel-researched steps were performed to extract the vessel segments: vessel segmentation, vessel skeletonization, identification of bifurcation/intersection points, and extraction of optimized vessel segments. Thus, those steps will be discussed in the method section as they are necessary steps to extract the vessel fragments.

Extracting the vascular tree skeleton from fundus retinal images is challenging due to the complicated and variable design of retinal arteries and veins, noise and poor illuminations, and lesions. Graph-based techniques, morphological operations, and machine learning-based approaches to vessel skeletonization have been proposed in recent years. For instance, Mahapatra et al. [21] proposed a vessel segmentation method based on an optimal enhanced Frangi filter and ranked spatial fuzzy C-means, then used a graph-cut technique for skeletonization, achieving the latest performance in terms of efficiency and accuracy.

Another recent work in this direction was based on an optimized BCOSFIRE filter for vessel segmentation [22], FCN deep learning for the segmented retinal vessels classification to arteries and veins [23], the segmentation of ring-cut, and calculating the iterative-AVR as a new method for measuring AVR that was the first in the field to be implemented and published for arteriovenous ratio [6]. Using a publicly available dataset, ’AV classification’ and retinal morphometry datasets form the stages in a decision support system for diagnosing eye diseases from retinal images; however, the skeletonization phase was only briefly explained and described.

Ouyang et al. [24] described a U-Net-based deep learning framework for retinal vessel segmentation and suggested modification to U-net by adding vessel-related local feature augmentation and focused techniques to increase the accuracy of vascular segmentation in the LEA U-Net model. U-net is a prominent convolutional neural network (CNN) model. During segmentation, local feature improvement modules gather comprehensive vessel information and attention mechanisms concentrated on the key vessel sections. LEA U-Net performs well in retinal vascular segmentation. It detects vessel features, resists noise, and bypasses picture fluctuations. The local feature improvement and attention techniques increase segmentation accuracy and vessel feature extraction.

Lyu et al. [25] introduced a benchmark dataset and assessment system. Image improvement and noise reduction precede vessel segmentation utilizing vessel-ness filters and post-processing. Vascular tree analysis extracts morphological traits and quantifies vascular tortuosity and branching patterns. The research provides a benchmark dataset, assessment methodology, and a high-performing retinal vascular tree analysis approach. The method’s shortcomings are to be tested with patients with severe retinal images that require more study. Yan et al. [26] proposed a vessel skeleton extraction approach based on multiscale line tracking. The method combines vessel enhancement and line-tracking algorithms to extract vessel skeletons. The main advantage is its ability to handle various vessel widths. However, it may need more accuracy in complex vessel structures. Chen et al. [27] introduced a hybrid approach combining vessel segmentation and skeletonization using a deep learning model. The method achieves high accuracy in vessel segment extraction but requires a large annotated dataset for training, which can be time-consuming and costly. Wang et al. [20] proposed a graph-based skeleton extraction technique that utilizes vessel segmentation and graph representation. The method achieves accurate vessel skeleton extraction and is robust against noise. However, it may encounter challenges in handling vessel crossovers and branching points. Hawas et al. [28] presented a vessel skeleton extraction method using a scale-space analysis and graph-cut optimization. The method achieves accurate vessel segment extraction and effectively handles vessel crossings. However, it may be computationally intensive and time-consuming.

In conclusion, retinal vessel segment extraction methods have seen significant progress in recent years. Each method discussed above has strengths and weaknesses, providing valuable contributions to the field. However, challenges related to vessel branch extraction, varying vessel widths, and computational efficiency persist. Future research should address these challenges and develop more robust and efficient vessel segment extraction methods.

In addition to the fact that blood vessels are visible on the retina’s surface, vessel morphology changes can also be identified using retinal images. Its variations are apparent indicators of the severity of various eye diseases.

Another example of morphological alteration is vessel tortuosity, defined as the presence of twists or turns in the morphology of the vessel [29]. As the tortuosity increases, the severity of eye diseases like center retinal vessel occlusion (CRVO) [30,31], diabetic retinopathy, hypertensive retinopathy [32,33,34], systemic hypertension [35], and retinopathy of prematurity (ROP) [36,37,38,39,40] increase. Ref. [41] is associated with female sex, elevated blood pressure, elderly age, and other cardiovascular risk factors. Several studies have surveyed and analyzed tortuosity metrics. For example, Ref. [42] categorizes the surveyed metrics as distance-based and curvature-based. Ref. [43] provides a comprehensive review of tortuosity formulas and their practical use, including an in-depth discussion of every particular method and its formula and a review and classification of the tortuosity databases used in those studies. Ref. [44] comprehensively analyzed the correlation between diabetic retinopathy and vessel tortuosity. The detected patterns of tortuosity, as described in [45], manifest themselves as:

• Tortuosity when the vessel appears in C- or S-shaped elongation;
• Looping when the S- or C-shaped with Multivessel symmetry sign;
• Coiling when the vessel is shaped with 360-degree turns in the vessel itself; and
• Kinking: when it manifests arterial angulation with an acute level.

Developing a CAD system for eye disease diagnosis is a difficult task. However, numerous metrics have been proposed to quantify tortuosity. Such systems have both advantages and disadvantages. Additional research must be conducted to reach a universal consensus on the most precise metrics, standardizing tortuosity grades and connecting such severity grades with every eye disease [44]. The techniques for measuring retinal vascular tortuosity can be divided into three categories:

• Arc to chord ratio methods;
• Curvature-based methods;
• Hybrid item methods.

The vessel tortuosity is error-prone to the noisy center line or skeletonization procedures. However, the fractal-wavelet-based algorithm may be used as it can be directly applied to segmented images. Furthermore, it analyzes the whole image and not only particular sample points on the vessel (curve). As tortuosity is a level 1 result of hypertensive retinopathy (HR), fractal analysis is a method that is more sensitive to small changes in vessel tortuosity. However, it is robust as accuracy remains unchanged regardless of the dataset size. Moreover, it extracts the full vasculature tree [46].

This study proposed a better approach for skeletonizing and extracting the right vessel segments from the five hundred and four (504) pictures in the AV classification dataset and measuring the improvement in the tortuosity calculation procedure before and after optimization. This paper used the distance metric from distance-based approaches and the inflection count metric from hybrid methods to examine the impact of optimized extracted vessel segments on enhancing the tortuosity measurement process.

3. Materials and Method

3.1. Material

The AV classification dataset described in [23] includes labels for two research problems, vessel segmentation, and artery-vein classification, which was developed specifically to satisfy the requirements of deep learning experiments and is also used in their evaluation. The proposed method was applied to the AV classification dataset retinal images. The AV classification dataset was created to satisfy the ‘ labels. These fundus images and their labels are 2002 by 2000 pixels. Each original retinal image has four labels, one that is colored to help with vascular segmentation research problems, and the other is monochromatic. In addition to these two labels, AV classification research challenges may be solved with two additional terms.

3.2. Method

The process of optimizing the vessel skeletonization and extracting the exact proper vessel segments, either root segments or between bifurcation points, cross-sections, or leaf segments, are summarized in Figure 1.

3.2.1. Vessel Segmentation

Although blood vessels’ separation from fundus images of the retina is difficult due to factors such as uneven illumination, the center light reflex, choroidal vascularization, poor contrast, impulse noises, and background homogeneity, numerous supervised and unsupervised methods exist for retinal vessel segmentation; the authors’ previous work in [22] was used for the retina’s effective vessel segmentation. The suggested vessel segmentation technique combines the optimization for B-COSFIRE with the proposed optimized filter parameters. The method was verified to be useful in vessel edge detection and overcoming the central light reflex problem. Furthermore, the suggested technique can appropriately segment all the vasculature, including the vessels, with the central light reflex. Therefore, it is included in the researchers’ method for automatic background removal and returning the monochrome image of the segmented vessels.

3.2.2. Vessel Skeletonization and Segments Extraction

Skeletonization is a mathematical operation performed in the images to enable the detection of an object’s center lines inside the image. In this research, it is the detection of each vessel’s centerline in the vasculature tree. In addition, it is a way of representing the significant topological features of the vessels to handle the operations related to the vasculature tree. Blum initially proposed a skeletonization transform in 1967 using a transformation called Grassfire [47]. Blum’s transformation algorithm starts from a random point and uses evolution propagation to transform and create the median access skeleton [47]. Most of the skeletonization literature approaches are inspired by Blum’s transformation method. However, skeletonization is widespread in several methods and procedures, and this can be grouped into three main approaches (distance-based approaches, iterative thinning approaches, and Voronoi-skeleton). See

Table 1. Topological and geometrical features of each skeletonization method type.

Skeletonization Methodology and Approach	Topological	Geometrical
(1) Distance-based transform	$\chi$	✓
(2) Voronoi-skeleton	✓	✓
(3) Thinning alliteratively	✓	$\chi$

Note: ✓ indicates the approach is classified under this skeletonization type, while $\chi$ indicates not.

:

In iterative thinning approaches, the researchers used morphological processing via iteratively repeated erosions of the object boundary to detect the image skeleton [48,49] as in Figure 2. While in distance-based approaches, other researchers used geometric distance-based methods by generating the distance map of the symmetrical shape parts, using different distance transformation methods (geodesic, Manhattan, and Euclidian) to directly compute the symmetric shape parts to identify the skeleton within, as in Figure 3. In addition, others recognized the skeleton points as centers of circles in 2D or the center of balls in 3D [50], as in Figure 4. While others employed non-iterative processes that generate a specific mid-line or centerline of the model, to be dilated immediately in one cycle, without looking at all the images in Figure 5. Others used the shape constraints to transform the problem into grassfire and then curve evolution.

Another way for vessel skeletonization is to consider the loci of centers of two-tangent circle rings with the shape boundary, which completely fits in the front [51]. As illustrated in the triangular form in Figure 4, the idea of including the contextual information about the complete branching segments with color-based features enhances the skeletonization and segment extraction Multiple methods were proposed for vessel skeleton extraction in [52]. First, the vascular skeleton is removed using a thinning process [53] that iteratively erodes one layer of border pixels from each connected component within the black-and-white image, considerably preserving its connectivity until each connected part is the skeleton itself.

3.2.3. Spur Pixels Cleaning

At the top of the MATLAB function bwmorph(), which has the below mathematical Formulas (1) to (4): When used with the ’thin’ option, bwmorph applies the following MATLAB logic [54]:

• In the first sub-iteration, delete pixel p if and only if the conditions A1, A2, and A3 are all true;
• In the second sub-iteration, delete pixel p if and only if the conditions A1, A2, and A3’ are all true,

where x1, x2,…, x8 are the values of p’s eight neighbors in counter-clockwise order, beginning with the neighbor to the east.

   Condition $A_1$:

$$X_H\left(p\right)=1where{X}_H\left(p\right)=\sum _a^b{v}_{i}and{v}_i=\left\{\begin{array}{c}1,if{x}_{2i-1}=0and(x_{2i-1}=1or{x}_{2i+1}=1)\\ 0,otherwise.\end{array}\right.$$

(1)

Condition $A_2$:

$$2\le min(n{1}_p,n{2}_p)\le 3where{n}_1\left(p\right)=\sum _{k=1}^4{x}_{2k-1}\wedge x_{2k}and{n}_1\left(p\right)=\sum _{k=1}^4{x}_{2k+1}\wedge x_{2k}$$

(2)

Condition $A_3$:

$$(x2\vee x3\vee \overline{x8})\wedge x1=0$$

(3)

Condition $A_{3^{\prime }}$:

$$(x6\vee x7\vee \overline{x4})\wedge x5=0.$$

(4)

In this work, the authors introduced an enhanced skeletonization method at the top of the MATLAB skeletonization step by eliminating the noisy spur pixels and smoothing the generated skeleton following the analysis.

The initial skeletonization results are illustrated in Figure 6. It has extra noise, with connected vessel segments on each vessel’s left and right sides. This was enhanced and removed by iteratively detecting the endpoints and eliminating them until the main vessel was reached. Then, the extra small segments were deleted from the vessel, and its pixels were no longer found as endpoints until the last of such noisy segments appeared.

3.2.4. Detection of the Vessel Tree Intersection/Bifurcation Points

Vessel branch point detection is required to measure, analyze, and quantify the tortuosity level. Firstly, a morphological operation was performed to detect the branch points, then dilate them. The researchers used the identified branch points pixels to separate the different segments, hence locating the vessel segment as connected components from the vasculature tree that do not include crossover or bifurcation points. The segment is either a starting root segment or a vessel part between bifurcation points, cross-sections, or ending leaf-points [55]. In contrast, the edges of the vessel tree vasculature are attained through a morphological operation that cleans the inner pixels and reserves the vessel’s pixels [55]. To detect the branch points, the researchers calculate the cross-point-number (CPN) using Equation (5) in [56], where P1–P8 are the 8-neighborhood pixel values counter clock-wise as illustrated in Figure 6 around the point $C_i$ that is under detection.

$$\begin{array}{c}\hfill CPN=\frac{1}{2}\sum _{i=1}^8\left(P_i-P_{i+1}\right).\end{array}$$

(5)

3.2.5. Eliminating the False Intersection Points and Segment Extraction

The MATLAB-implemented function for image skeletonization performs the morphological thinning of the image boundaries. This step achieves the image’s skeleton However, it has a drawback: generating a spur pixel when making the vessels’ skeleton creates confusion about whether it is the daughter branch of the vessel or just a spur extra pixels (Figure 7).

To solve this problem, the researchers have improved skeletonization results to perform the following steps in the Algorithm 1 pseudo-code.

A phenomenon has been noticed in the above MATLAB function output, where one pixel remains in the vessel skeleton from the root of each deleted spur piece, causing L-shaped angles in the vessel curve. Furthermore, this extra pixel has an impact on the subsequent steps as it causes adding two extra branch points in the MATLAB branch point extraction function, BWmorph(), for ’branch points’ detection, and it also causes breaking the vessel segment into multiple subsegments.

Algorithm 1: Optimization pseudo code for fragment extraction

3.3. Measuring Retinal Vessel Tortuosity to Quantify the Impact of the Enhanced Vessel Fragments

Extracting the vessel fragments from the fundus retinal images is a prerequisite to calculating vessel tortuosity. It is the level of twistedness in the retinal vessels. It is an important sign of the severity of diabetic nephropathy and hypertensive retinopathy. Many formulas are proposed and used to measure the tortuosity in retinal images. In this work, the researchers have selected the inflection count metric tortuosity measures to calculate the tortuosity of all the vessel segments in the AV classification dataset (See Figure 8).

It is essential to extract the vessel segments from the retinal image to proceed in the calculation of the tortuosity inflection count metric (ICMn) using the below formulas in Equations (6) to (8):

The inflection counts metric uses the number of angles by which the vessel twists in different directions as it gives a strong indication of the twists that happen in the vessel multiplied by the ratio of the geodesic vessel length devised by the straight line distance length between the start and the end of the vessel segment, see Figure 8 and Equation (6) for the calculation formula.

$$ICM=(Inflectio\mathit{n\_p}oints+1)*\frac{Arc}{Chord}$$

(6)

$$IC{M}_{Before}=\left\{A\right|R^3:A\Rightarrow (Ar{c}_{i,j},Chor{d}_{i,j},IC{M}_{i,j})\}$$

(7)

$$IC{M}_{After}=\left\{A\right|R^3:A\Rightarrow (Ar{c}_{i,j},Chor{d}_{i,j},IC{M}_{i,j})\}.$$

(8)

The above equations will be applied twice on the 504 images of the vessel segmentation dataset. As a result, the researchers generated a feature set for quantifying the impact of improving the vessel segment extraction on the tortuosity calculation results and to proceed in measuring the process capability before and after the improvement.

3.4. Performance Metrics

The process capability index ($CpK$) from the prominent six-sigma standard was employed in this study to quantify the improvement gained from the improved vessel extraction steps. $CpK$ is a real number that measures a process’s or a product’s feature’s performance per the criteria stated. Generally, the distance between the process specification limits is divided by the distance between the natural tolerance limits, as defined by six sigma process units. These indices let us determine how well the process meets the criteria [57]. A process is competent if the actual values’ quality characteristic of the generated data come within the lower and upper specification limits with a high degree of probability. Various statistics, such as Cp and $CpK$, assess a manufacturing process’s capability. As the $ICM$ tortuosity calculation becomes evident, the researchers must add the optimization steps to obtain the $ICM$ of each vessel segment. Adding the optimization steps to the optimized vessel extraction procedure will generate a new, enhanced method for calculating tortuosity. Therefore, the researchers have two processes: process1 is the tortuosity measurement before enhancement application, and process2 is the tortuosity measurement after enhancement. The experiment was performed to extract all the vessel segments using the old approach and to calculate the $ICM$ tortuosity for all the vessel segments extracted before applying the optimization. The second scenario is to generate all the vessel segments using the new optimized approach and to calculate $ICM$ tortuosity for all the vessel segments.

The quantified tortuosity values create a feature set containing ’ICM before’, where the tortuosity is measured for all the retinal images in the selected ’AV classification’ dataset.

Based on the $CpK$ results, the sigma level was identified, and finally, the researchers compared the two sigma levels for the two processes to identify the process improvement level. The above steps are illustrated in the flowchart diagram Figure 9 to represent the quantification process visually. It is worth noting that this quantification is just one example of the benefits of vessel extraction optimization. Furthermore, it was used here to have another view of the importance of the reached results in addition to being presented visually in Section 4.3 in the results section. Alternatively, Process or machine capability indices are quantitative assessments of how well machines and processes can operate compared to specifications. The capability index assumes that the process output parameter follows a normal distribution. The capability index combines the normal distribution parameters (the mean $\overline{X}$ and standard deviation) into a single parameter by comparing the observed process characteristics to a theoretical three-sigma process. The upper and lower specification limits are the theoretical process limits ($USL$ and $LSL$) [58].

$CpK$ considers the process mean and analyzes the process distribution about the actual state of the process. The ratio of $CpK$ to $Cp$ directly indicates how far off-center the process is performed. It presupposes that the process output follows a normal distribution. The calculated $CpK$ and $CP$ values are identical if the process variance lies between the defined limits. However, when the variance deviates from the center of the specification, it is penalized proportionally to its degree of deviation. $CpK$ is extremely beneficial and extensively employed. A $CpK$ greater than 1.33 indicates that the process is almost capable. Less than 1.33 indicates that the variation is too large relative to the specification or that its location is offset from the specification’s center. It could be a combination of both location and width. $CpK$ measures the three deviations between the process mean and the closer specification limit. To accurately measure $CpK$, it is necessary to use data that follow a bell-shaped “normal” (Gaussian) distribution. Others consider it an approximation. Only when the process is perfectly centered does $CpK$ equal $Cp$. $Cp$ represents the maximum possible value for $CpK$[59]. For processes where the mean of the measured input variable matches the value of the targeted mean, the process capability Cp index indices may be derived using Equation (9). It may be used to determine the process-capability $Cp$ index for processes where the mean of observed parameters matches the value of a target:

$$CP=\frac{USL-ICM}{2*k*\sigma },$$

(9)

where $USL$ and $LSL$ are the specification limits, k is 3 (the process capability), and $\sigma$ is the observed standard deviation of the input variable data in Equation (10) of the input variable data.

$$\sigma =\sqrt{\frac{{\sum }_{i=1}^n{(x_i-\mu )}^2}{N}}.$$

(10)

Suppose the measured parameter’s mean is not equal to the value of a target. If the mean is not centered within the specification limits, the corrected capability index is $CpK$. The process capability index $CpK$ (see Equation (11)) is used to measure the process capability before optimization and after optimization using the below formulas borrowed from the six sigma framework [57].

$$CpK=min(\frac{\left|USL-\overline{X}\right|}{3\sigma },\frac{\left|\overline{X}-LSL\right|}{3\sigma }).$$

(11)

$CpK$ will be measured before applying the optimization steps for all the 504 images (12) and similarly will be repeated in (13) after applying the optimization steps

$$Cp{K}_{Before}=\left\{R\right|R:IC{M}_{i,j}\Rightarrow Cp{K}_{i,j}\}$$

(12)

$$Cp{K}_{After}=\left\{R\right|R:IC{M}_{i,j}\Rightarrow Cp{K}_{i,j}\}.$$

(13)

In order to use $CpK$, data should be normally distributed; this can be achieved by applying a Box–Cox transform on the data of $ICM$ in (14) and (15)

$$ICM{n}_{Before}=\{R|R:IC{M}_{i,j}\stackrel{BoxCoxTransform}{\to }ICM{n}_{i,j}\}$$

(14)

$$ICM{n}_{After}=\{R|R:IC{M}_{i,j}\stackrel{BoxCoxTransform}{\to }ICM{n}_{i,j}\}.$$

(15)

If the tortuosity $ICM$ values for the 504 retinal images are not normally distributed, then transformed will convert it to become normally distributed. The Box-Cox transformation normalizes the variable optimally, bypassing the need to randomly test several transformations to find the optimum alternative. Moreover, it converts non-normal data on the required positive response variable X into normal data using Equation (16).

$$x^{\lambda }=\left\{\begin{array}{cc}\frac{x^{\lambda }-1}{\lambda }& \lambda \ne 0\\ \ln \left(x\right)& \lambda =0.\end{array}\right.$$

(16)

It operates by adjusting the shape and dispersion of the data via a power transformation. When the Box-Cox transformation is used to fit the data to a ’normal’ distribution, this may be used to test for normality. If the modified data have a normal distribution, the original data are probably also normal. The Box-Cox transformation may enhance data linearity and convert skewed data into a more symmetrical distribution. Furthermore, increase the precision of statistical tests like linear regression and $CpK$. The following assumptions need to be fulfilled before using the $CpK$ [60] that are:

• Comprehending the foundations of process capability analysis and its corresponding measurements, accumulating the process data, and computing the essential statistics;
• If the data is not normally distributed, apply Box-Cox Transform, as it is a prerequisite to $CpK$ calculation;
• If the data are not normally distributed, apply the Box-Cox Transform, as it is a prerequisite to $CpK$ calculation capability.

After validating the above assumptions, $CpK$ was calculated for ICM_fore and ICM_after; the results were compared, and the $CpK$ outcomes of the process were interpreted.

4. Results

This research experiment’s results are presented as follows; the visualization of the process steps results is illustrated in Section 4.1. Followed by depicting the enhanced extracted vessel segments in Section 4.2 and Section 4.3. In addition, presenting how the newly optimized vessel segments helped to gain more accurate results of the tortuosity metrics in Section 4.4. Moreover, Section 4.5 about the updated AV classification dataset with the new set of enhanced vessel fragments dataset.

4.1. Visualization of the Results of the Process Steps

Figure 10 illustrates a high-level overview of each step’s visual results in extracting the vessel fragments and calculating the tortuosity.

The proposed enhancement steps were added after the step in Figure 10d,e. Figure 11a shows the existence of a spur problem, while Figure 11b illustrates the spur removal improvement.

After performing the standard segmentation and skeletonization, the small spur segments were cleaned from the skeleton. This result leads to the segments containing ‘L’-shaped consecutive sets of pixels in the skeleton (See Figure 12).

Such pixels cause the wrong breakdown of the vessels to the segments at each pixel in the “L”-shaped corner. Hence, the researchers have detected and eliminated those pixels from the skeleton (see Figure 12c) to achieve smooth, proper, complete vessel segments (see Figure 13).

After applying the proposed method of enhancing the skeletonization to obtain a proper segment for each vessel branching or crossover, the researchers noticed the optimization difference by detecting the right full segments before and after optimization of the optimized skeleton in Figure 14, as these pixels cause the fragmentation of the corresponding segments at false branch points.

4.2. Vessel segments identification

As described in Section 3.2.2, the list of the proper segments is achieved by extracting the branch points from the optimized skeleton generated in step 6 in the process in Figure 1 above and performing the logical operation between the optimized skeletonization and the branch points image complement; Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15 and Figure 16 illustrate the enhancement of the segment’s identification.

4.3. Visualization of the Extracted Optimized Segments

Figure 15, Figure 16 and Figure 17 show examples of the proposed method implementation on the 504 images of the AV classification dataset and the improvement of the skeletonized vessels that led to proper vessel segments instead of having the initial broken vessel segments from the MATLAB bwmorph() function results of datasets.

4.4. Results of Calculating the Tortuosity before and after the Improvement

The inflection count metric $ICM$ and its Box–Cox-transformed reading is applied before and after applying the optimized steps to all vessel segments in each retinal image included in the AV classification dataset to measure the impact of optimizing the vessel segment extraction. A sample of this work is presented in

Table 2. A sample from the tortuosity results before and after enhancing the vessel segments.

	Before Enhancing Vessel Segments Extraction			After Enhancing Vessel Segments Extraction
Image_No	Segments Count	ICM	ICMn	Segment Count	ICM	ICMn
102	86	85.0	3.4	79	305.9	2.2
110	47	69.5	4.1	39	134.5	2.1
114	63	93.6	3.9	52	217.5	2.1
139	54	73.2	4.0	43	227.4	2.1
140	76	44.3	3.7	59	110.8	2.1
271	68	127.5	3.8	59	226.3	2.2
309	46	30.7	4.2	36	112.9	2.2

below:

4.5. AV Classification Dataset

The AV Classification dataset is enhanced by having another set of skeleton images for its original retinal images (See Figure 18). This set of new images can be available for researchers of such research problems by email.

5. Discussion

This section measures and discusses the advantages of optimal vessel segment extraction. The process capability index compares the tortuosity calculation’s improvement before and after the vessel extraction optimization. The improvement of the quantification of the tortuosity computation is found in Section 5.1.

Its worth indicating that the steps of image processing and machine learning needed to achieve the above results were performed using MATLAB 17b. In contrast, the below analysis of the process capability index was performed using MINITAB version 22.

5.1. Quantifying the Tortuosity Calculation Improvement

$CpK$ was used in this work to quantify the effect of vessel extraction on the values calculated for $ICM$ tortuosity. Other process-related metrics (sigma level, yield, and DPMO) were calculated and compared in addition to $CpK$. The two columns of $ICMn$ before and after, shown in Table 2, represent the results of calculating the $ICMn$ tortuosity before and after optimization for all 504 images and were used as input to the $CpK$ index calculations. In order to quantify the impact of the enhanced vessel segments, the stages above will be described and illustrated in Section 5.1.1, Section 5.1.2 and Section 5.2.

5.1.1. The Impact on the Vessel-Segment Length

After optimizing the vessel segment extraction, the extracted vessel segment becomes longer since each segment represents a portion of the vessel from the vessel start/intersection/bifurcation point to the subsequent intersection/bifurcation point/vessel end. Therefore, if the researchers take the mean length of all the vessel segments in the retinal image before and after optimization and plot those averages in a linear regression line, they obtain the following graph, demonstrating that the average vessel segment length has increased for all 504 retinal images. These results are promising to enhance the computation results for tortuosity (see Figure 19).

5.1.2. The Impact on the Tortuosity Inflection Count Metric

Figure 20 shows two different linear box plots for 504 retinal image’s inflection count metric values before and after optimization. The optimized vessel fragments process connected multiple broken pieces of the same vessel fragment into one connected piece from its start to its end; this decreased the total number of segments, ending up with an image with a longer segment and, from the graph, it can be seen that the average and mean of the ICM_after value have greater ICM_before average and mean. This suggests that the tortuosity readings are generally greater after optimization.

5.2. The Impact on the Process Capability Index of Calculating the ICMn

A set of assumptions must be validated before calculating the $CpK$ of the two retinal vessel tortuosities, ’ICM before’ and ’ICM after’, for optimized extracted segments from the 504 images. Kotz and Montgomery [60] state that the following important assumptions are to be asserted and confirmed.

• The process output must be statistically controlled;
• The distribution of the quality attribute is normal;
• Observations must be random and unrelated to one another.

Section 5.2.1, Section 5.2.2 and Section 5.2.3 discuss the above assumptions’ validity to calculate the process capability index before and after vessel segments enhancement.

Table 3. Sample data collected from the ’ICMn before’ and ’ICMn After’, each sample is a group of five readings, where the mean (X-bar) and the range (R) are calculated for each sample.

ICMn_Before									ICMn_After
S#	C1	C2	C3	C4	C5	X-bar	R		S#	C1	C2	C3	C4	C5	X-bar	R
1	3.1	3.1	3.4	3.3	3.9	3.35	0.80		1	2.1	2.0	2.1	2.1	2.2	2.09	0.24
2	3.3	3.8	3.2	3.0	3.4	3.34	0.75		2	2.0	2.1	2.2	2.0	2.0	2.07	0.19
3	3.3	3.3	3.0	3.3	3.3	3.26	0.31		3	2.1	2.1	2.1	2.1	2.1	2.08	0.02
4	3.8	3.6	3.2	3.2	2.9	3.36	0.91		4	2.1	2.1	2.0	2.0	2.1	2.07	0.11
5	3.2	3.2	3.1	3.3	3.1	3.20	0.22		5	2.2	2.1	2.0	2.0	2.1	2.07	0.28
6	3.1	3.3	4.3	3.2	3.7	3.52	1.15		6	2.0	2.0	2.2	2.0	2.2	2.09	0.20
7	3.3	3.7	3.7	2.8	3.4	3.38	0.89		7	2.1	2.1	2.2	2.0	2.1	2.07	0.21
8	3.7	3.0	3.4	3.0	4.2	3.44	1.29		8	2.1	2.0	2.0	2.1	2.2	2.08	0.22
9	2.9	3.4	3.0	3.4	3.2	3.18	0.51		9	2.1	2.0	1.9	2.0	2.0	2.02	0.15
10	3.1	3.5	3.2	3.8	3.7	3.46	0.76		10	2.0	2.1	2.1	2.1	2.1	2.08	0.13
11	3.2	3.5	3.6	3.1	3.4	3.35	0.48		11	2.1	2.0	2.1	2.0	2.1	2.06	0.08
12	3.7	3.2	3.4	2.9	3.1	3.27	0.87		12	2.1	2.1	2.0	2.0	2.0	2.02	0.11
13	2.8	4.0	3.2	2.7	3.5	3.25	1.29		13	1.9	2.1	2.1	1.9	2.0	2.01	0.25
14	3.3	4.1	2.9	3.4	3.6	3.44	1.18		14	2.0	2.2	1.9	2.1	2.1	2.07	0.30
15	3.0	3.4	3.3	3.2	3.5	3.29	0.53		15	2.0	2.1	2.0	2.1	2.1	2.06	0.16
16	3.1	3.4	3.0	3.6	3.5	3.35	0.57		16	2.1	2.2	2.0	2.0	2.1	2.07	0.16
17	3.2	3.1	3.2	3.1	3.5	3.21	0.48		17	2.2	2.1	2.0	2.2	2.1	2.14	0.24
18	3.8	3.5	3.8	3.3	3.1	3.50	0.70		18	2.1	2.1	2.1	2.2	2.1	2.09	0.10
19	3.3	3.9	3.2	3.4	3.0	3.38	0.87		19	2.1	2.1	2.0	2.1	2.0	2.06	0.15
20	3.6	3.2	3.2	3.5	2.9	3.29	0.69		20	2.0	2.2	2.1	2.1	1.9	2.06	0.21

is a pre-processing sample from the ICMn_Before and ICMn_After data to draw the run charts and probability plots for validating the assumptions and calculating the process capability indices. The table shows a moving window of five readings taken to calculate their average (X-bar) and range R-bar and the total target averages of X-bar and R-bar according to the Equations (9)–(11).

5.2.1. Verifying the Statistical Stability of the Tortuosity Measurement Processes before and after Optimization

The control chart can be used to review the data and to ensure process stability. Figure 21a,b show that all displayed sample ranges and average scores are inside the control boundaries on the R-Chart and the X-Bar chart, with no sign of clustering, run, trend, or shift. Therefore, the process is statistically controlled and influenced only by rare variations. In other words, the process is steady over time.

5.2.2. Verifying the Normality Assumption

Histogram and probability plots and transformation were utilized to test the normality assumption of the calculated tortuosity $ICMn$ results. The histogram and normal probability plot were created via MINITAB version 22 Statistical software.

The tortuosity results in each column were checked to determine whether they were normally distributed. If not, it was transformed using the Box-Cox transform to become normally distributed to fulfill the requirement to measure the process capability index as per [57,60]. Figure 22a,b depict the histogram of the $ICM$ tortuosity results before and after vessel extraction optimization; both histograms show that the data are skewed to the right and, as the p-value is less than 0.05, the data are not normal. It is crucial to remember that, for a more accurate process capability index, the data should be more symmetrical in a normal distribution shape [60], and this can be achieved by applying the Box-Cox transform for the two columns’ $ICM$ before’ and ’ICM after’ to become a normal distribution shape (see Figure 23a,b) and the new columns ’ICMn before’ and ’ICMn after’ values become ready for the process capability calculation.

Figure 23a,b show that the two histograms of the $ICMn$ before and after data become normally distributed after transforming the data using the Box-Cox transform.

Moreover, the normal probability plots are depicted in Figure 24a,b, highlighting that the normal probability plot test results for the tortuosity $ICMn$ data before and after optimization reveal that, via the Anderson Darling test statistic, a mean: (3.045, 2.053), a standard deviation: (0.5094, 0.07516), and the significance threshold (0.05) is increased to (0.829, 0.257). Furthermore, p-values of (0.032, 0.721) for the ICMn_before and ICMn_after indicate that the transformed process data are now normally distributed about the mean. This indicates that the information has been successfully transformed. Consequently, it has been determined that the transformed data are generated through a normally distributed method.

5.2.3. Verifying the Randomness Assumption

The run chart in Figure 25 discovered the four main types of non-randomness: trend, mixture, oscillation, and cluster patterns, where:

• The trend: is a consistent upward or downward shift in data;
• A mixture: is distinguished by the lack of points along the center line;
• The oscillation: is the data swinging up and down;
• Finally, clusters are collections of connected points on a single side chart center line.

Run chart interpretation: the p-values for trends, clustering, mixtures, and oscillation are all larger than the alpha value of 0.05. The actual number of runs is similar to what was predicted. As a result, the data are judged to be independent and random.

The actual number of runs is similar to the anticipated number. Therefore, $ICMn$ tortuosity-reported results data observations are random before or after enhancing vessel segment extraction.

After validating the three assumptions, the researchers can calculate and compare the process capability indices of the $ICMn$ tortuosity computation process and analyze the $Cp$ and $CpK$ values for each process before and after adding the enhancement steps.

5.3. Process Capability Index-$CpK$

$Pp$ and $Ppk$ denote the process performance, while $Cp$ and $CpK$ denote the process capability. $Cp$ and $CpK$ are the most important capability indices. $Cp$ indicates how much the process data can fit within the criteria. In contrast, $CpK$ indicates whether the overall mean is centered. That is why $CpK$ is used when the process output distribution is skewed. If the aggregate average falls in the middle of the specification, $Cp$ and $CpK$ will be identical. If $Cp$ and $CpK$ scores are dissimilar, this indicates that the aggregate mean is not exactly centered. The greater the two score values’ disparity, the more the mean deviated from the overall average. This concept is depicted in Figure 26a,b. It is noticed that the overall and within bell curves are almost over each other in Figure 26a as $|Cp-\mathrm{CpK}|=1.03-0.89=0.14$. In contrast, the overall and within bell curves have almost deviated from each other in Figure 26b as $|Cp-CpK|=|1.46-1.67|=0.21$ since the between bell curve values deviate from the mean more than the overall bell curve.

Table 4. Process capability vs. sigma level and process yield before and after vessel segment enhancement.

The Process Capability ($\mathit{CpK}$) of Measuring the Tortuosity $\mathit{ICMn}$	$\mathit{CpK}$	Sigma Level	DPMO	Conforming (Yeild%)
Before vessel segments optimization	0.89	2.7	115,000	88%
After vessel segments optimization	1.46	4.39	1866	99.77%

compares $CpK$ values, sigma level, DPMO, and process yield for the data before and after the optimization steps. The $CpK$ value increased from 0.89 to 1.46, indicating that the process became less dispersed around the mean. Thus, the process can produce $ICMn$ tortuosity readings within expected specification limits (see Figure 27). Furthermore, the increase in the $CpK$ value from 0.89 to 1.46 moves the process capability from being almost capable to being a capable process as per the process capability index scale (See

Table 5. Process capability/sigma level/DPMO Scale [59,60].

Is Process Capable?	$\mathit{CpK}$	Sigma Level	DPMO	Standard Dev Compared to Specification Limits
Not Capable	0.33	1	691,462	Higher
	0.67	2	308,538
Almost Capable	1	3	66,807	Lower
	1.1	3.3	35,930
	1.2	3.6	17,864
	1.3	3.9	8198
Capable	1.33	4	6210	Lower
	1.47	4.4	1866
	1.6	4.8	483
	1.7	5	233
	1.8	5.4	48
	2	6	3.4

).

According to Table 4 and Table 5, the proposed improvement shifted the process capability from almost capable to a capable process.

Table 4 presents a big improvement in the DPMO level of the tortuosity process defects, as the defects are reduced from 115,000 per million before optimization to 1866 per million after optimization. The sigma level was calculated based on the number of defects that occur per one million opportunities (DPMO). There are six main categories, with the sixth being the most capable process.

• 1 Sigma: This sigma level permits 691,462 defects per million chances.
• 2 Sigma: This standard allows for 308,538 defects per million opportunities.
• 3 Sigma: 66,807 defects per million opportunities are permissible at this sigma level.
• 4 Sigma: 6210 defects per one million opportunities are acceptable.
• 5 Sigma: 233 defects are permitted at this sigma level.
• 6 Sigma: 3.4 defects per million opportunities are observed at this level.

Finally, the proposed improvement increased the sigma level of the tortuosity calculation process from 2.7 to 4.39 and improved the confirming readings (yield) from 88%.

Figure 26a,b depict the results of calculating the process capability indices $Cp$, $CpK$ for the tortuosity calculation process before and after enhancing the vessel segment extraction. After adding the vessel segment enhancement steps in the tortuosity calculation process, the difference between $CpK$ and PpK values converges to zero when the researchers compared their values before and after optimization, indicating that the process is under greater statistical control, as the sample standard deviation and $sigma$ become almost identical. $Cp$ and $CpK$ compare the performance output consistency to the process output average performance. The letter ’k’ represents a centralizing factor. The index considers the possibility that process output data are not centered. $CpK$ provides information regarding the future performance of a process when its results are under statistical control. By comparing the $CpK$ value after vessel segment enhancement (1.46) to its value prior to applying the proposed new steps (0.89), it can be seen that the process becomes more capable of calculating the tortuosity $ICMn$ formula after improving the vessel segment extraction and the vessel segments become one long connected segment from the start till the end of the segment, as opposed to being broken to multiple fragmented sub-segments before.

The new set of enhanced skeletonized 504 images was added as a subset of the researcher’s work on retinal vessel morphometry (RVM) generated datasets [6]. Figure 27 illustrates how the bell curve of the $ICMn$ data after the optimization becomes closely tight around its mean (2.05) compared to the ’ICMn_before’, which was spread away around its mean (3.05).

This effort aligns the extraction of vessel segments with the true definition of the vessel segment, allowing for a more accurate calculation of vessel tortuosity.

6. Conclusions

In this research, the authors have contributed to solving the problem of breaking the vessel segment into various parts while extracting the vessel segments from the vasculature skeleton, particularly when the vessel segment is between two bifurcation points and is broken into multiple pieces. The vessel segment extraction from the skeleton of the vasculature tree is a crucial step in quantifying vessel tortuosity for diagnosing several eye diseases. However, inaccurate extraction of the vessel segments can cause errors by increasing the count of segments, thereby impacting the accuracy of the tortuosity measurement. The author’s contribution was the proposal of an improvement to the MATLAB skeletonization bwmorph() function by identifying and eliminating an ’L-shaped’ set of pixels in the skeleton, which is the root cause of the incorrect identification of branch points, resulting in the incorrect extraction of segments at said points. Also, the Six Sigma process capability index method was utilized to quantify the improvement in calculating the level of vessels’ tortuosity before and after optimizing the extracted vessels. As a result, the vessel segment optimization improved the sigma level from 2.7 to 4.39. In contrast, the confirming yield improved from 88 percent to 99.77 percent.

Finally, the authors generated optimized vessel segments for each of the 504 images in the AV classification dataset, available in monochrome and color formats. The images are organized in distinct folders for future use in research projects. Researchers can request access to the dataset by contacting the authors and citing this work in their academic research.