OCTAVA: An open-source toolbox for quantitative analysis of optical coherence tomography angiography images

Optical coherence tomography angiography (OCTA) performs non-invasive visualization and characterization of microvasculature in research and clinical applications mainly in ophthalmology and dermatology. A wide variety of instruments, imaging protocols, processing methods and metrics have been used to describe the microvasculature, such that comparing different study outcomes is currently not feasible. With the goal of contributing to standardization of OCTA data analysis, we report a user-friendly, open-source toolbox, OCTAVA (OCTA Vascular Analyzer), to automate the pre-processing, segmentation, and quantitative analysis of en face OCTA maximum intensity projection images in a standardized workflow. We present each analysis step, including optimization of filtering and choice of segmentation algorithm, and definition of metrics. We perform quantitative analysis of OCTA images from different commercial and non-commercial instruments and samples and show OCTAVA can accurately and reproducibly determine metrics for characterization of microvasculature. Wide adoption could enable studies and aggregation of data on a scale sufficient to develop reliable microvascular biomarkers for early detection, and to guide treatment, of microvascular disease.


Introduction
The microcirculation comprises the smallest elements of the circulatory system, a dense network of arterioles, capillaries, venules, and lymphatic vessels with a diameter of less than 150 μm [1]. These small vessels, the microvasculature, account for about 99% of blood vessels in adults and play a key role in oxygen transport and nutrient delivery to the tissue, as well as in waste and carbon dioxide removal. Many studies have suggested that dysfunction in the microcirculation may be as important as dysfunction in larger blood vessels in the context of cardiovascular pathophysiology [2]. However, due to the lack of accessible/available diagnostic modalities tailored to the microvasculature, knowledge of its precise role in disease pathogenesis and progression is lacking and no specifically microvascular therapies are currently available [3]. With the introduction of optical coherence tomography angiography (OCTA), the study and visualization of the microvasculature is well underway in the human retina and, to a lesser extent, in human skin and animal models of human disease [4][5][6]. Standard "structural" OCT uses interference to measure the echo arrival time and intensity of backscattered light to generate cross-sectional and volumetric images of optical scattering in tissue [7]. By using differences between sequential OCT images (usually cross-sectional B-scans) at the same location, caused by motion in the sample, OCT volumetric images can be processed to generate volumetric or two-dimensional representations of blood flow, angiograms, which distinguish structures containing flowing blood from the surrounding static tissue [8]. OCTA angiograms typically comprise en face maximum intensity projection (MIP) images, where MIPs have been shown to be superior to other projections [9]. OCTA is an ideal technique for clinical imaging because it does not require exogenous contrast agents; contrast is generated from blood flow, largely from red blood cells. Based on OCTA images, it is possible to describe the microvascular network architecture and vessel morphology and introduce, describe, and quantify various metrics to characterize them. Some commonly used metrics include: vessel area density, total vessel length, and mean vessel diameter, all within a region of interest [10]. It is also possible to undertake a temporal analysis of pulsatility towards characterizing pathophysiology or arterial stiffening [11], or measure microvascular response to external stimuli, such as heating or pressure, to enhance diagnostic accuracy [12,13].
There are a number of commercial and non-commercial, clinical and laboratory-based OCTA instruments and systems that enable visualization and characterization of microvasculature. These systems invariably use different imaging engines, imaging protocols, image processing methods, and metrics to describe the microvasculature. As a result, comparing study outcomes obtained with different instruments and generating large databases across instruments and institutions is not possible. Indeed, it has been proven that quantifiable metrics obtained using different OCTA instruments are not readily comparable [14,15]. As a more specific example, differences in binarization thresholding methodologies have been proven to significantly influence the quantification of OCTA metrics in healthy eyes [16,17]. Standardized approaches to process and analyze OCTA images are needed. It is crucial that metrics generated by the research and clinical community are accurate, consistent and reliable if these metrics are to be used to generate agreed-upon biomarkers associated with diagnosis, monitoring, and treatment guidance of disease [18,19].
Fortunately, there is a growing interest within the research community in delivering universal tools that enable visualization and quantification of microvasculature. A number of tools originally developed for microscopy can be applied to OCTA images [20][21][22][23][24], even if none have yet been optimized for OCTA. Some toolboxes limited in scope, in terms of application or instrument, have been made available for evaluation of OCTA images [25][26][27], with differences in processing methods and quantitative metrics. The lack of consistency between such toolboxes limits their ability to compare results acquired with different instruments, in different studies, or in different institutions, and, thus, the ability to draw conclusions about the biological relevance of the results, pointing to the need for a more versatile and universal alternative.
We have developed and validated a new, easy to use, open-source software toolbox, OCTA Vascular Analyzer (OCTAVA), towards standardization of OCTA image analysis and characterization. Standardization implies choices of the free parameters in such a toolbox. We have identified from the extensive existing literature the best algorithms for OCTA image analysis and characterization. We compare in detail the performance of the five best performing segmentation algorithms and select eight metrics by assessing quantitatively and qualitatively their capability to visualize and characterize microvasculature. Finally, we implemented the optimal workflow in OCTAVA: i) pre-processing and segmentation of en face OCTA images; ii) identification of the vascular network and nodes using skeletonization; iii) automated measurements of the length, diameter, and tortuosity of the each identified vessel; and iv) generation of the microvascular network architecture and vessel morphology metrics: vessel area density, vessel length density, mean and median vessel diameter, total vessel length, branchpoint density, mean tortuosity and fractal dimension. Using images from various sources, we demonstrate that OCTAVA can accurately and reproducibly determine metrics for characterization of microvasculature independent of instrument and application.

OCTA imaging and human subjects
OCTA images used for toolbox development, validation and optimization were acquired using the multi-beam VivoSight Dx (Michelson Diagnostics Ltd, Maidstone, Kent, UK). This sweptsource OCT instrument allows imaging of microvasculature of the skin with 20 kHz line-scan imaging speed and provides an imaging resolution of 5.5 μm and 7.5 μm in the axial and transverse directions, respectively. Eight healthy individuals aged between 25 and 62 years old were enrolled in this study. This study was approved to be undertaken through the issue of a favorable ethical opinion by the University of Surrey Ethics Committee (FHMS 19-20 060). Written, informed consent was obtained from all participants in adherence to the Declaration of Helsinki. For all participants, the skin on the dorsum of the right hand, between the thumb and forefinger, was imaged. The handheld OCT probe was positioned for imaging on the skin through a plastic cap to reduce motion artifacts and maintain a constant distance between the imaging probe and the skin. Volumetric images were acquired over a 5 mm by 5 mm area with spatial sampling of 4.4 μm along the fast axis and 41 μm along the slow axis. Three images were acquired sequentially at the same location without adjusting the placement of the handheld probe. The built-in software "VivoTools" was used to generate volumetric representations of blood flow based on pairs of cross-sectional B-scans. A surface detection algorithm was then applied in MATLAB 2020a (The MathWorks, Inc., Natick, Massachusetts, USA) to identify and flatten the curved skin surface. The top 175-μm-deep section of the image was removed to exclude the avascular epidermis and dermal papillae. Dermal papillae were excluded since vessels in the dermal papillae are oriented along the imaging axis; thus, they appear as dots and their interconnectivity with the network is difficult to analyze. The depth range of our images corresponds with the superficial plexus, a dense network of vessels oriented parallel to the skin surface. Two-dimensional (2D) en face OCTA angiograms were generated from the volumetric angiograms in MATLAB using a maximum intensity projection (MIP) over a physical thickness of 500 μm in depth (unless otherwise specified), assuming an average group refractive index of 1.4.

Comparison of segmentation algorithms
Previous work has noted significant variability in resulting metrics using different segmentation algorithms [28][29][30][31]. As such, we have evaluated the performance of five different image segmentation algorithms to determine the optimal approach for OCTA MIP images in skin: a global thresholding approach using the built-in "convert to mask" function in ImageJ; k-means [32]; iterative self-organizing data analysis technique (ISODATA) [33]; adaptive thresholding [34]; and fuzzy thresholding [35]. These algorithms were selected by reviewing viable approaches in the retinal OCTA literature based on the number of papers that use them and the reported accuracy of the results. Out of tens of algorithms, we selected five from a preliminary assessment as demonstrating the best performance [28][29][30].
The five segmentation algorithms were first optimized empirically; any free parameters were adjusted to maximize recognition of vessel structures and minimize misidentification of background signals as vessels, assessed manually. Then, the performance of the five segmentation algorithms was compared using selected OCTA MIP images with different densities of vascular structures to assess how accurately the binarized image represented the OCTA data. Even if it is currently used as a gold standard in many other works, we opted not to use manual segmentation of the vascular network as the ground truth because it generally has poor accuracy and precision compared with automatic segmentation [29], and the results vary depending on the experience of the person performing the segmentation [36]. Instead, the algorithms were compared empirically and quantitatively using the metrics of vessel area density (VAD), network connectivity factor (CF), and repeatability. VAD is described in detail in Section 2.5. The network connectivity factor, CF, quantifies the percentage of identified vessels that are connected to the main network as defined by the formula CF ¼ 1 À S i S t , where S i is the number of isolated elements in the image (segments with neither endpoint connected to the network) and S t is the total number of identified vessel segments in the image. Repeatability of VAD was assessed using the within-subject standard deviation (described in Section 2.6) [37].

OCTAVA development environment
Software development was mainly undertaken in MATLAB for ease of optimization, graphical user interface (GUI) development, and customization. We utilize the ImageJ-MATLAB package to access ImageJ (US National Institutes of Health, Bethesda, Maryland, USA) libraries from MATLAB, which allows us to fully integrate the image processing capabilities in ImageJ and the automation and data management tools in MATLAB [38]. This approach allows us to build on existing functionality within ImageJ while giving us the flexibility to investigate different segmentation algorithms and collation of metrics, which is more straightforward to implement in MATLAB. Our OCTAVA software adopts and integrates the Angiogenesis Analyzer [22] developed within ImageJ [39] with several important modifications. First, we have developed and optimized image pre-processing and segmentation in MATLAB, which ensures good quality recognition of vessels in OCTA MIP images. Second, we have expanded the available quantitative metrics describing vasculature beyond those generated by Angiogenesis Analyzer to include additional information about the length, diameter, and tortuosity of individually identified vessels and the interconnectivity of the network, which enables us to perform more sophisticated statistical analyses on vessel data and to generate additional metrics beyond those available with other software. Finally, we have developed a dedicated GUI, shown in Fig 1, that allows automated segmentation, identification of vessel components, and compilation of results from within a single interface. Because our software is open source and developed in MATLAB, it can easily be modified to adapt to the needs of the research and clinical community while still being easy to use without modification. A MATLAB license is required to modify the software; however, a standalone version of OCTAVA is available that can be used without MATLAB.

OCTAVA software
2.4.1. OCTAVA GUI. The OCTAVA GUI provides a single interface for users to process OCTA MIP images and to view results. The user has the option of selecting a region of interest (a subregion of the image) for processing (rectangular or circular) and rescaling the number of pixels in the image-either to increase processing speed or to improve the precision of the metrics (red box in Fig 1). OCTAVA has two different operating modes: one which allows the user to process individual images with the aim of optimizing the analysis protocol for specific image types, and a batch-processing mode which allows the user to analyze multiple images without additional user input (purple box in Fig 1). In both cases, the same image processing workflow is utilized; each of these steps is described in detail in the following sections and visualized in Fig 2. Similarly to REAVER [23], we give the user an option for manual curation to allow for correction in the case of improper automated segmentation or vessel identification (yellow box in Fig 1). Manual curation is performed by manually adding or removing items from the list of generated regions of interest (ROIs) within ImageJ. Many metrics have been implemented to enable the identification of those most meaningful for use in diagnosis and monitoring of disease and response to treatment. These will be discussed further in Section 2.5.
2.4.2. Image processing workflow. OCTA MIPs are pre-processed using MATLAB to generate a binary mask using the following procedure. OCTA MIPS are first assessed by the user as having sufficient image quality to proceed. A 2D Frangi "vesselness" filter is applied since the objects of interest are blood vessels [40] (Fig 2C). The Frangi filter is commonly used in analysis of angiography images since it reduces the impact of intensity variations along a vessel and suppresses background noise, thereby improving image segmentation [41]. The 2D "vesselness" filter is used rather than a 3D filter for accessibility, since many OCTA instruments limit access to the volumetric data (volumetric data is acquired, but often not made readily available to the user). A multi-scale approach is used to minimize artificial vessel dilation (optimization of the filter is discussed further in S1 File). The OCTA MIPs are segmented into "vessel" and "not vessel" regions ( Fig 2D) using the optimal segmentation algorithm as determined by our analysis in Section 3.1. The binarized image is then skeletonized using a 3D thinning algorithm in ImageJ ( Fig 2E) [41,42] and a heatmap of vessel thickness is generated ( Fig 2F) [43]. The algorithm for measuring local thickness is described in Section 2.5. Nodes are identified based on their diameter relative to other nearby structures and removed so they do not impact diameter data. The dimensions of individual vessel segments are recorded and written to a spreadsheet for post-processing and compilation of results using MATLAB.

Identification of network elements.
OCTAVA identifies different elements that make up the network and displays this information in a color-coded overlay of the skeletonized image on top of the original image ( Fig 2G). Linear regions are separated into three categories based on their interconnectivity with the network: segments, colored yellow, are connected to other vessels at both endpoints (via nodes); branches (green) are connected to . Colored borders in (b-g) correspond to the step with the matching-colored arrow or box in (a). OCTA MIP image (b) is loaded into OCTAVA. Images are first pre-processed in MATLAB to remove noise and enhance the signal intensity of vessel-like structures using a Frangi filter (c). Images are then segmented using the chosen segmentation algorithm to generate a binary mask (d). The image is then sent to ImageJ, where the network is skeletonized (e). The thickness and interconnectivity are measured (f), and the ROIs and network elements are identified and classified based on their thickness and interconnectivity. An overlay image is generated (g) which helps the user confirm that an accurate map of the network architecture was measured. Colors in (f) correspond to the local vessel thickness. Colors in (g) correspond to different architectural components of the network including segments (yellow lines), branches (green lines), mesh regions (light blue lines), isolated elements (dark blue lines), and nodes (red and blue circles). The insets in (g) allow a closer examination of the identified network elements. Finally, the outputs of the ImageJ analysis are sent to MATLAB, which generates and compiles metrics of the network. The large images in (b-g) are 5 mm × 5 mm. The insets in (g) are 1.6 mm × 1.6 mm.
https://doi.org/10.1371/journal.pone.0261052.g002 another vessel on one end but have one free endpoint; isolated elements (dark blue) are not connected to other vessels on either end. Isolated elements below a certain diameter (specified using the "twig size" control in the OCTAVA GUI) are likely noise and are excluded from the analysis. Twig size, specified in pixels, may require adjustment based on the presence of noise in individual datasets. Nodes, marked in the image with blue and red circles, are counted but are removed from the diameter measurements to avoid artificially overestimating vessel diameters. Mesh regions, which are areas completely enclosed by the vessel network, are identified in cyan. Examples of each of these structures can be seen in the insets in Fig 2G.

Quantitative metrics for characterization of microvascular network
The microvascular network forms a complex architecture of interconnected vessels. The complex spatial network undergoes remodeling in embryonic and adult tissues, as well as in quiescent and pathological states. Such complexity signifies that full characterization of changes in vessel morphology will require multiple metrics. We chose eight metrics to comprehensively characterize microvascular network architecture based on an extensive qualitative review of the literature, balancing coverage and overlap of characteristics [13,[44][45][46][47][48]. These metrics allow users to obtain comprehensive information about the microvasculature [49]. Table 1 and Fig 3 give an overview of our selected metrics, how they are calculated, and indications of biological relevance. We note that the potential biological relevance is not exhaustive-this is an open area of research which will be facilitated by our software-and is included here as a brief justification of the breadth of metrics we chose to implement. All metrics are based on the 2D OCTA MIP image.
VAD is calculated from the binary mask by counting the total area of the identified vessels (white pixels) divided by the total number of pixels. The vessel length density (VLD) is calculated by counting the total number of black pixels in the skeletonized image, which represents the sum of the length of all the vessels (total vessel length) and dividing by the total number of pixels in the image. For the purposes of the VLD calculation, a vessel diameter of 1 pixel is assumed in order to have units of %. For other purposes, vessel diameters are calculated using

PLOS ONE
the local thickness algorithm in ImageJ applied to the binarized OCTA MIP image [43]. Overall mean and median diameter values are calculated as the average diameter of the vessel segments. The fractal dimension is calculated using the box counting method [50,51], and is reported as a number between 0 and 2 with a standard deviation; a higher value and lower standard deviation corresponds to more fractal-like structures. Branchpoint density is calculated by counting the number of nodes in the image and dividing by total vessel length. Tortuosity is calculated from the skeletonized image using the arc length-over-chord ratio, where the arc length is the total length of the identified segment between two nodes, L s , and the chord length is the length of a straight line between the two endpoints, L C , so that Tortuosity ¼ L s L c À 1 [52]. Additional methods for quantifying tortuosity for different applications will be incorporated into future versions of OCTAVA. In addition to the average length, diameter, and tortuosity measurements, we present histograms of the distribution of these values within the image. Previous work has shown that the distribution of these values can provide further information, for example, for distinguishing changes in total number of perfused vessels versus vessel dilation for increased VAD [53].

Evaluation of OCTAVA
The OCTA images of human skin described above were used to evaluate various pre-processing filters and segmentation algorithms to optimize the OCTAVA workflow. These images were also used to evaluate the repeatability of the OCTAVA metrics. The repeatability of each metric was calculated based on the within-subject standard deviation (S w ) method introduced by Bland and Altman [37]. The standard deviation of repeated measures for each subject was calculated, squared to obtain the variance for each subject, and then the square root of the average variance for all subjects gives the measurement error, S w . The coefficient of repeatability (CR) is defined as 2.77S w and 95% confidence intervals (CI) as: CR � 1:96 � ðS w = ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 2nðm À 1Þ p Þ, where n is the number of subjects and m the number of measures for each subject [54]. Additionally, we have used OCT images of intralipid flow in microfluidic devices obtained with a spectral domain OCT system described elsewhere [55] and a simulated OCTA MIP image to evaluate the accuracy of algorithms for metrics calculation. Manual characterization of the microfluidic device was performed in ImageJ. Exact metric values for the simulated OCTA MIP image were generated using a priori knowledge of the image structure. As well, to demonstrate the utility of the software in various types of OCTA images, we have used OCTA (MIP) images from several previously published works. These include: OCTA images of human skin collected with a laboratory-based polarization-sensitive OCT system [56], OCTA images of mouse ear skin acquired with the Telesto spectral-domain OCT system (Thorlabs, Newton, NJ, USA) [57], and images of the human retina obtained using the RTVue XR Avanti (Optovue, Inc., Fremont, CA, USA) [58].

Validation of OCTA segmentation algorithms
A comparison of the performance of each of the optimized segmentation algorithms for three images with different vascular network densities and SNRs can be seen in the top part of Fig 4. In general, the built-in global thresholding algorithm within ImageJ does not accurately represent the complexity of the vascular network for in vivo OCTA MIP images of skin. It can be seen by visual inspection that the ImageJ global thresholding approach tends to underestimate the density of the vascular network and the resulting segmentation is poor. For images with high vessel density, such as the bottom row, all four of the remaining algorithms performed reasonably well; however, for images with lower VAD or poor SNR, ISODATA and adaptive thresholding tend to add additional structures within noisy but largely uniform regions of the image regardless of the kernel size specified. This means that these algorithms may have limited utility for images of skin but may still be useful for retinal images, which tend to have fewer uniform regions. The fuzzy thresholding and k-means algorithms replicate the vascular network in the binary mask over a larger variety of images. The five segmentation algorithms were also compared quantitatively using two metrics: VAD compared with manual calculation from the OCTA MIP image and network connectivity factor. The results of the quantitative comparison can be seen in the bottom of Fig 4. The ISODATA and adaptive thresholding algorithms led to a high connectivity factor for all three images, which is not representative of the apparent differences in the images (higher CF does not necessarily correspond to better segmentation). The k-means and fuzzy thresholding algorithms performed similarly on both VAD and connectivity factor and led to a VAD closer to the manual calculation, indicative of good segmentation. Based upon this analysis performed over a large range of images, we have determined that the fuzzy thresholding algorithm provides the most accurate segmentation for our OCTA MIP images of skin and is applicable to a wider variety of images than the other algorithms. As a result, we have only implemented the fuzzy thresholding and adaptive thresholding algorithms in the final version of OCTAVA. While adaptive thresholding did not perform well with our images, it is used commonly in retinal OCTA; we included it only as a basis for comparison for users who may be more familiar with this algorithm.

Validation of vascular architecture metrics
To validate the accuracy of our generated metrics, we used two starting images: a conventional OCT en face image of a microfluidic device with known channel diameters; and a simulated OCT MIP image with hand-drawn structures mimicking a vascular network (Fig 5). The image of the microfluidic device was collected with a spectral domain OCT system described in [55]. The original OCT image depicted six small channels merging into one large channel. Although it conveys the advantage of known sizes, it was difficult to analyze this image since the pattern of the microfluidic device is not representative of a vascular network. Therefore, we synthesized a new image by copying and tiling the original image in a grid pattern to better mimic a network of vessels. For both images, the binary mask (Fig 5B and 5G) was generated using the fuzzy thresholding algorithm. Due to the bimodal distribution of vessels in the microfluidic image, the Frangi filter was not used; the performance of the Frangi filter is highly dependent on the range of vessel sizes within an image, so it is not possible to optimize the filter parameters for two distinct diameter ranges without distorting the image [59]. For the simulated OCTA MIP image, intensity variations were introduced along the simulated vessels so that the performance could be evaluated including the image enhancement step (without intensity variations, the Frangi filter would have no effect). The Frangi filter was used with a σ max value of 7 and was optimized using a procedure similar to that described in S1 File. The pixel size in each image (4 μm for the microfluidic image and 9.3 μm for the simulated OCTA MIP) determines the expected error bounds for diameter and length measurements.
The segmentation was evaluated by plotting the intensity of the OCT image and the binary mask along a vertical or horizontal line. The results demonstrate that the segmentation step does not impact the apparent channel diameters in most cases, although the Frangi filter may increase the diameter of small, low-intensity vessels depending on the chosen value of σ max . This is further confirmed by the histogram of channel diameters for both images (Fig 5E-5J). As expected, the histogram of diameters for the microfluidic device indicates a bimodal distribution with peaks at the correct values of 50 μm and 300 μm representing the small and large channels. Measured diameters are typically within the error bounds determined from the pixel size with a few exceptions. For example, some of the channels in the microfluidic device are measured to have larger diameters than expected due to the functioning of the algorithm at junctions between the small and large channels, as can be seen from the overlay image ( Fig  5C). Someone manually characterizing the image might identify a node at the endpoint of each of the 50-μm channels before any broadening or change in direction indicating the junction with the larger channel; the automated characterization often includes part of this transition region as part of the smaller channels and places the node closer to the 300-μm channel. This can be seen most clearly in the top and bottom channels of each row. As a result, the measured average diameter of the 50-μm channel is larger than it would be if the transition region were identified as an independent segment. For comparison, the nodes identified by OCTAVA were taken to be the segment endpoints for the manual characterization of the network. The thickness maps (Fig 5D-5I) demonstrate accurate diameter measurements throughout both images. The histograms of the lengths also match closely with the expected values (not shown).
The metrics generated by OCTAVA were compared with values calculated manually for both the microfluidic image and simulated OCTA MIP image ( Table 2). The VAD was calculated manually using the mean intensity of the image as a global threshold. This is a valid approach for these images since there is a clear difference in intensity level between the network of channels and the background. Diameter and length measurements for the microfluidic images were measured manually in ImageJ in order to manually calculate VLD, mean and median diameter, and branchpoint density; for the simulated OCTA MIP image, exact diameter and length measurements were calculated based on a priori knowledge of the image structure. In all but one measurement, the relative difference between the automatic and manual metric calculations did not exceed 10%, although OCTAVA tends to overestimate the vessel diameter when compared with manual characterization based on the way it identifies nodes and vessel endpoints as described above. The automated method we use for measuring the vessel diameter (the local thickness algorithm in ImageJ) is commonly used for measuring the diameter of structures within an image (irrespective of the image content), however, it is difficult to replicate it exactly in manual analysis. This difficulty, in conjunction with the placement of nodes described above, and the limitation based on pixel size, accounts for the differences in the diameter measurements.

Repeatability of OCTAVA metrics
Three repeated scans of the same skin location of eight study participants were used to validate the repeatability of OCTAVA metrics. Table 3 and Fig 6 summarize the outcome of the analysis. In Table 3, mean is the mean of all 24 measurements and measurement error (ME) is the mean of S w for all the participants. The factors that determine CR are related to differences in the raw OCTA images and the analysis. The repeatability of the software to obtain consistent

Application of OCTAVA to other OCTA instruments and imaging applications
The performance of OCTAVA was further evaluated using several other published OCTA MIP images of skin and retina [56][57][58]. The results are presented in Fig 7. In all cases, the segmentation (Fig 7B, 7F and 7J) and network identification (Fig 7C, 7G and 7K) showed highquality reproduction of vessel morphology. A heatmap was generated of the locally identified

PLOS ONE
vessel diameters throughout the image and the range of measured diameters is within the expected range for each case. This outcome is notable given that the same software was used to process and analyze all three images. The selected images came from several published works demonstrating the applicability of the software on different instruments, imaging targets, and applications. Given the range of metrics implemented within OCTAVA, we were able to describe additional metrics of these images beyond those presented in the original publications.

Discussion
While we do not propose any new algorithms here, we argue that standardization is the main attribute required to move the field forward. Many algorithms and processing workflows have already been proposed, usually proprietary to a research group or instrument vendor, to the detriment of open science. Overall, there is no generally accepted best workflow for segmenting and analyzing OCTA data [31,59], either overall, independent of the imaging target, or for a particular application. Our innovation comes from making informed choices through comparison of some of the best performing algorithms with the ambition of establishing a widely accepted and utilized standardized processing workflow suitable for all OCTA image types. To maximize its ease of uptake, we have packaged our proposed solution as an open-source, easyto-use toolbox which is accessible to users at all levels of technical competency. This is important because it can provide a starting point for new entrants to the field of OCTA and is accessible to clinicians and end users while still being customizable for users seeking to adapt the software to their own needs. All these features together represent the power of OCTAVA.

Workflow optimization and comparison with other software packages
The vessel enhancement and segmentation steps are key for accurate identification of vessels, and different segmentation algorithms produce significant variation in measured metrics.
With this in mind, we will briefly compare OCTAVA with other reported software used for analysis of OCTA images. This comparison will contextualize OCTAVA as a step towards a standardized workflow that could be used for all OCTA images.

Vessel enhancement.
In angiography more generally, vessel enhancement (filtering) is often used to improve the visibility of vessels. Here, we briefly justify our choice of the Frangi filter for vessel enhancement. Hessian-based filtering (including the Frangi filter) has been commonly used in magnetic resonance angiography, ultrasound angiography, and photoacoustic angiography [64][65][66]. Specifically for OCTA, there are many possible approaches, each with advantages and disadvantages, including Hessian-based filtering, rodfiltering, top-hat filtering, optimally oriented flux enhancement, Gabor filtering, weighted symmetry filtering, and active shape models [25,59]. Among these approaches, Hessian-based filtering has been shown to result in the best quality of final image [25]. Other works have employed a Frangi filter to improve visibility of vessels in OCTA images of the retina [49], choroid [67], and skin [41]. Despite its widespread use, the Frangi filter has also been criticized in some studies for introducing errors in the vessel architecture, either by missing vessels rendered in the image with low SNR or generating spurious vessels depending on the structure of the background noise [68,69]. Even using the multiscale Frangi filter approach, it has been shown that the results are highly dependent on the range of vessel sizes within the image [70]. Our study of Frangi filter behavior with OCTA images of the skin demonstrates that for σ max values in the range of 1-8, vessel dilation is minimal. However, parameter optimization must be adapted for each imaging application, which is a limitation on standardization. We finally note that modifications to this algorithm to further optimize it for OCTA images represents a possible further improvement [70].

Filling in gaps in vessels:
Manual curation vs automation. Improving network connectivity by extrapolation is another issue that must be considered as a potential image pre-processing step. It is expected in a real microvascular architecture, that all vessels are fully integrated in a network. Therefore, breaks in vessels and disconnected endpoints far from the image boundaries do not accurately represent vessel morphology. AngioTool (an open source tool optimized for confocal fluorescence images), among other reported tools, has incorporated the ability to extrapolate to fill in gaps in vessels to improve connectivity [21]. A required trade-off in such approaches is the acceptable degree of image manipulation balanced against reasonable expectations about the vessel network. For example, the discontinuity may be due to a change in direction or depth of the vessel, thus, adding the connection may be erroneous. Furthermore, since OCTA can only visualize perfused vessels with flow rate above a minimum threshold, non-patent microvessels and obstructions appear as discontinuities and could be overlooked if the gaps are filled in automatically [71]. At this stage, we have chosen not to incorporate automated gap-filling to avoid ambiguous manipulation of images. In some cases, where added connection is justified, the user may use the manual curation option to add connections. Although, to date we have found that (depending on the vessel density) these cases cause minimal change in the values of the overall metrics (see Section 4.2).

Image segmentation: Computer-learning vs traditional algorithms.
Some recent studies have investigated the use of deep-learning algorithms to improve vessel segmentation [25,27]. For example, Stefan and Lee have compared their deep learning segmentation algorithm (optimized for images of the mouse brain) with global thresholding and adaptive thresholding for volumetric OCTA image processing [25]. Although they demonstrated good performance, deep-learning and other computer-learning-based approaches for image segmentation are incompatible with the development of a standardized workflow since computer learning algorithms must be trained using a particular set of images.
Several recent studies in retinal OCTA have evaluated the performance of various traditional segmentation algorithms and found significant variability between them [17,[28][29][30]. Of these, the fuzzy means algorithm has shown promise for imaging in the choroid and retina [67,72]. Combined with the results presented here in skin, the fuzzy means algorithm is beginning to stand out as a good option for standardizing the processing workflow.

Choice and definition of metrics.
Many metrics have been used for characterizing microvasculature; an overview of common metrics can be found in [44]. One of the challenges impeding widespread adoption of clinical OCTA is the wide variability in how these metrics are defined and calculated. This is compounded by the fact that most software available for analyzing vascular networks does not allow access to the raw data used to generate those metrics. For example, AngioTool is commonly adapted for OCTA images but only provides a summary of metrics and does not provide information about individually identified vessels [21]. By contrast, OCTAVA makes the measurements from all identified vessel segments available to enable further statistical analysis and generation of additional metrics. Similarly to [25], OCTAVA uses a graphical approach to present the outcome of the quantitative analysis through histograms of diameter, length and tortuosity, which allows for further analysis of the mechanisms driving changes in metrics which characterize the image as a whole, such as vessel density and mean diameter.

OCTA image artifacts and the impact on segmentation and metrics
One of the main limiting factors in the application of OCTAVA is the presence of motion artifacts. Motion artifacts are a main limiting factor in any segmentation-based analysis because the algorithm does not differentiate artifacts from vessels. Casper et al. have recently proposed a segmentation technique which incorporates a refinement step to help correct for defocus and blurring caused by motion, which may prove helpful [69]. In any case, if motion artifact exceeds a certain level (such as the example in Fig 8A and 8B), there is no choice but to exclude the images from analysis. Thus, due attention should be paid to stabilize the patient/instrument interface to minimize motion artifact in the first place. For images with infrequent motion artifacts, OCTAVA allows for improperly identified features to be manually excluded from the analysis using manual curation. One example of this and the impact it has on the image metrics is shown in Fig 8C-8H. The impact of the motion artifact correction on the metrics depends on the overall information content of the image. However, we found that in cases where it was feasible to manually remove motion artifacts, the impact on the value of the metrics was minimal. Manual curation can also be used to remove artifacts from other sources, such as floaters in the eye or edema. Future versions of OCTAVA will include a quality index to indicate the degree of motion present in the image and determine automatically whether the image is a candidate for manual curation or should be removed from the statistical analysis.

Future directions
We are considering several improvements to the segmentation and analysis workflow in future versions of OCTAVA. One consideration is the ability to further segment images based on vessel diameter, which some studies suggest could bring extra diagnostic power [27,73,74]. Further segmentation based on vessel diameter may be used to differentiate between different vessel groups (e.g., venules, arterioles, and capillaries). It is known that different eye diseases and severity stages can affect the arterioles and venules differently [75,76], so it would be useful to be able to distinguish between them.
Another consideration is 3D-based analysis of the OCTA data. One inherent limitation of analyzing MIPS is that it is impossible to discern between two vessels that intersect or cross at different depths [44]. As shown by Sarabi et al., 3D analysis enables better visualization and captures details of vascular geometry not seen in 2D analysis [77]. However, challenges persist in 3D analysis of OCTA images due to projection artifacts inherent in OCTA imaging [78]. Taking into consideration how common and clinically applicable 2D analysis is, we think a standardized image processing workflow in 2D is still a significant step towards facilitating open data sharing and cross-comparison and pooling of data from different studies.
While our initial optimization has focused on skin, we recognize the wide application of OCTA imaging to retinal analysis, which requires the ability to measure several additional metrics including the area and perimeter of the foveal avascular zone, and the vessel density within several subregions of the image [79]. Retinal analysis will be delivered in the next version of the software. Motion artifacts, poor-quality images, and manual curation. Images with significant motion artifact (a) must be removed from the analysis since the segmentation algorithm cannot distinguish between the vessels and motion artifacts; the resultant segmented OCTA MIP (b) does not accurately represent the microvascular architecture. For images with minimal motion artifacts, such as the images in (c) and (d), the individual motion artifacts can be removed manually. (e) and (f) show insets of (c) and (d), respectively, demonstrating the removal of the motion artifact using manual curation. Depending on the VAD, individual motion artifacts likely have little impact on the generated metrics, as shown in the histogram of diameters (g) and the table of generated metrics with and without manual curation (h). Note that in (g), the columns representing the results with manual curation are always equal in height or shorter than the automatically generated values since no vessels were added by manual curation in this instance. Square OCTA MIP images in (a-d) are 5 mm × 5 mm and comprise an axial range of 500 μm. The images in (e) and (f) are 0.65 mm × 5 mm. https://doi.org/10.1371/journal.pone.0261052.g008 Finally, one of the goals of our research is to create multi-source image databases available to be mined to find better biomarkers of health and disease. For example, with the creation of large databases, we will provide enhanced opportunity for rapidly emerging AI-based image classification schemes to identify biomarkers.

Conclusion
We have developed a fully integrated and easy-to-use software tool, OCTAVA, for quantitative analysis of OCTA MIP images based on a range of metrics. We have demonstrated that quantification of the metrics with OCTAVA is accurate and repeatable. By showing that the software can be applied to OCTA images acquired from a range of instruments on different imaging targets, we have demonstrated OCTAVA's versatility. We believe OCTAVA is an important step towards standardization of the processing workflow for OCTA data. Such standardization will lead to wide-ranging benefits in its application.