Automatic 3D adaptive vessel segmentation based on linear relationship between intensity and complex-decorrelation in optical coherence tomography angiography

Introduction

Port wine stain (PWS) is a congenital cutaneous capillary malformations composed of ecstatic vessels, and the vascular targeted photodynamic therapy (PDT) has been widely used in the clinical treatment of PWS (1-4). Both the PWS diagnosis and PDT management could be improved by the motion-contrast optical coherence tomography angiography (OCTA) which is capable of visualizing the tissue microvasculature and perfusion down to capillary level in a non-invasive, label-free and high-resolution manner (5-7).

The motion-contrast OCTA records the depth-resolved scattering signals by measuring the echoes reflected from tissue, and then enhances the dynamic scattering of the red blood cells (RBCs) over the static background tissue by computing motion index (such as variance, difference or decorrelation) of the recorded backscattering signals (8,9). Using the motion index as the grayscale value of each voxel, a 3D OCTA image is generated. Furthermore, to reproducibly and objectively quantify the outcomes of OCTA in scientific research and clinical application, various quantification metrics of vascular morphology have been proposed, including vessel diameter index (VDI), vessel skeleton density (VSD), vessel area density (VAD), vessel perimeter index (VPI) and vessel complexity index (VCI) (10-15). The fundamental of calculating these metrics effectively is an accurate segmentation, which classifies each pixel as either a vessel or background. To obtain such a segmentation map, manual labeling by experts is the gold standard. However, since this is a time-consuming procedure and the outcome might vary from expert to expert, efforts have been focused on the automatic extraction of the vascular structures.

To achieve the automatic segmentation, a binarization threshold is usually applied to already generated OCTA images to assign all image pixels to be either black (background) or white (vessel), resulting in a binary image. Because the OCTA contrast mainly originates from the computed motion index, an empirical global motion-threshold has been used to extract the vessels from the background (i.e., binarization) (13). The fixed empirical threshold can be further improved with the automatic techniques such as histogram threshold using Otsu’s method (16-18), which adjusted the threshold for each image individually. Furthermore, local adaptive threshold approaches have been proposed for the processing of OCTA images (10,14,19), which produced better results than the global approaches. However, most of the local adaptive thresholds were applied on a 2D (x-y plane) enface maximum intensity projection (MIP) image regardless of the depth (z) dimension (12,13,17,20).

In addition, the histogram of the motion index is also affected by the local intensity or signal to noise ratio (SNR) of OCT signal (9,21). To accurately classify the SNR-dependent motion index, a simple approach is assigning an intensity-threshold to remove all the voxels without sufficient SNR (22-24), and then applying the motion-threshold on the MIP image for binarization. However, the OCTA motion index has a complicated dependence on the intensity of OCT signal, thus, a fixed intensity threshold will either include many strong noise-induced decorrelation artifacts or exclude lots of valid dynamic signals (9,25). Gao et al. (21) calculated the OCTA decorrelation values (SSADA signals) at different OCT intensity levels in the foveal avascular zone (FAZ), and concluded a numerical relation between decorrelation and intensity through linear regression. Accordingly, an intensity-adjusted decorrelation threshold was applied to the MIP images for a reliable OCTA quantification. However, their method is based on the MIP OCTA image and neglected the variation of SNR in the depth direction.

Recently, Huang et al. (9) derived an asymptotic linear relation of decorrelation to inverse SNR (iSNR) based on a multi-variate time series (MVTS) model, and accordingly they defined a linear 3σ boundary to distinguish the dynamic voxels from the static and noise in the iSNR and decorrelation (ID) space and to create SNR-adaptive OCT angiograms, termed as ID-OCTA. However, the conservative 3σ criterion cannot directly used as the binarization threshold for the purpose of minimizing binarization errors, because it was originally designed for generating OCTA image with high flow contrast and lots of low grayscale residual static tissues were reserved.

In this work, we proposed a binary image similarity (BISIM) method for automatically determining an optimal threshold in ID space entirely based on the acquired OCT and OCTA images. Firstly, we presented the ID-BISIM method. Then, human skin experiments were performed to validate benefits of the proposed ID-BISIM threshold both qualitatively and quantitatively compared with conventional 2-step intensity and binarization thresholding method in the decorrelation-contrast OCTA. Besides, vascular quantitative metrics were analyzed to demonstrate the potential of PWS evaluation with OCTA. Finally, the value of our work, advices for practical implementation and further improvement of the study were discussed.

Methods

Subject population and ID-OCTA

This study was approved by Ethics Committee of Chinese PLA General Hospital (NO.: 0404) and informed consent was taken all the patients. Ten subjects aged 5–15 years old with facial PWS were enrolled in this study. OCT scans were performed on two sites of cheek skin for each subject: (I) the lesion skin, and (II) the normal skin (the symmetrical site of the lesion on cheek). Each condition was scanned twice with their quantification results averaged to reduce the fluctuation of measurements, and a total of 4 datasets were acquired from each subject.

The skin imaging was performed with a lab-built OCT system, which was primarily based on a typical spectral domain configuration (26). A stepwise raster scanning protocol (x-y) was used for volumetric imaging over a field of view of 2.5 mm × 2.5 mm × 2.5 mm (z × x × y), with 300 A-lines per B-scan (fast-scan, x direction), 5 repeated B-scans per step (slow-scan, y direction) and 300 steps per volume, corresponding to a total acquisition time of ~8.3 s. Structural image was generated based on the intensity I(z,x,y) of OCT signal. As shown in Figure 1A, representative layered structures of human skin were clearly visualized in the structural cross-section, including epidermis and dermis. The OCT depth profile revealed that the signal intensity decreased exponentially with the increase of the penetration depth due to the light attenuation in the tissues, resulting in a decay of the SNR versus the penetration depth (insertion in Figure 1A).

Figure 1 ID-OCTA of human skin using a SNR-adaptive dynamic-static/noise classifier. (A) OCT structural cross section (i.e., intensity mapping). (B) Decorrelation mapping calculated with a temporal spatial kernel. The region with low SNR level (indicated by the red oval) also presents high decorrelation value, generating noise-induced dynamic artifacts. (C) Dynamic flow signal (red area) identified with the linear 3σ boundary (see the red line in D). (D) ID space mapping. Red line is the linear 3σ boundary. (E) Gaussian-filtered (sigma = 1, filtersize = 5) 3D ID-OCTA image. Decorrelation value is encoded with color. (F) Enface projection of ID-OCTA image. OCTA, optical coherence tomography angiography.

Angiographic images were created by calculating the change of OCT signal between successive tomograms taken at the same location. Here the inter-frame decorrelation D(z,x,y) was computed as a motion index to quantify the RBCs-induced dynamic changes (27), as shown in Figure 1B. However, it is impossible to extract the dynamic flow simply based on a threshold in decorrelation dimension, because the decorrelation value of noise region is also high, which results in pseudo-dynamic artifact (such as the region indicated by red oval in Figure 1B). To characterize the dependence of decorrelation on SNR, Huang et al. (9) proposed a MVTS model and found that the mathematical expectation of the decorrelation has a linear relation with iSNR Eq. [1] and the standard deviation σ is described as Eq. [2]:

$$D\to iSNR,a.s.,$$[1]

$${\sigma }_D=c_v\overline{D},$$[2]

where $iSNR=\frac{s^2}{I}$, s² refers to the intensity of system noise and I is the intensity of OCT signals. a.s. denotes convergence with probability one. c_v is the coefficient of variance (CoV) of decorrelation, $\overline{D}$ and is the mean value of the decorrelation. Accordingly, SNR-adaptive angiograms (see Figure 1C) were created by defining a linear 3σ boundary (see Figure 1D) to distinguish the dynamic voxels from the static and noise in the iSNR-decorrelation (ID) feature space, and the 3D ID-OCTA imaging and corresponding enface projection were obtained (Figure 1E,F). However, the conservative 3σ criterion cannot directly used as the binarization threshold for the purpose of minimizing binarization errors, because it was originally designed for generating OCTA image with high flow contrast and lots of low grayscale residual static tissues were reserved.

Binarization thresholding and quality metric

We proposed a SNR adaptive binarization method based on the linear boundary of static signals in ID space, which was illustrated in Figure 2. To quantitatively evaluate the linear SNR adaptive threshold, we defined an ID threshold D_T (see the red line in Figure 2F):

Figure 2 Proposed ID-binary image similarity (ID-BISIM) thresholding method for data-based 3D-adaptive binarization. (A) Binary image sequence created by changing the ID-threshold α_T from 1° to 90° with an interval of 1°. Typical binary OCTA cross-sections when (B1) α_T=5°, (B2) α_T=25° and (B3) α_T=40°. Cross-sectional maps of morphologic vector corresponding to (C1) B1, (C2) B2 and (C3) B3. Subtraction of the paired vector maps was performed to characterize the local structure differences between the two binary images: (D1) intra-class, (D2) inter-class, (D3) inter-class. (E) Plot of BISIM index versus ID- threshold α₁ and α₂. The BISIM index is minimal when α₁=7° and α₂=37°. (F) Voxels were projected into the ID space. ID-BISIM threshold α_T=7°, as indicated by the red solid line. (G) Binarized OCTA cross-section with the ID-BISIM method.

$$D_T=\cot \left({\alpha }_T\right)·iSNR$$[3]

where α_T is the ID-threshold evaluating the slope of the ID line D_T. The pixels were assigned to be white if its ID-threshold, $\alpha \in \left[0,{\alpha }_T\right]$ otherwise be black. Changing the threshold α_T from 1° to 90° with an interval of 1°, a binary image sequence can be generated (see Figure 2A) and can be roughly divided into three classes: (I) the binary image has only dynamic flow signal when ${\alpha }_T\in \left[1°,{\alpha }_1\right]$ (see Figure 2, B1); (II) the binary image has both flow signal and static tissue when ${\alpha }_T\in \left({\alpha }_1,{\alpha }_2\right]$ (see Figure 2, B2); (III) the binary image has flow signal, static tissue and noise when ${\alpha }_T\in \left({\alpha }_2,90°\right]$ (see Figure 2, B3). The binary images of the same class have a similar structure while the images of different classes exhibit different structures. To estimate the intra-class similarity and the inter-class difference, we defined a morphologic vector array $\overrightarrow{v}\left(\alpha ,z,x\right)$ to represent the structure of the binary image $B\left(\alpha ,z,x\right)$. As shown in Figure 2C, selecting a window from the binary image and setting the center of the window as the origin, the vector $\overrightarrow{v}\left(\alpha ,z,x\right)$ is the coordinate center of non-zero pixels in the window.

$$\overrightarrow{v}\left(\alpha ,z,x\right)=\stackrel{\left(k-1\right)/2}{\sum _{h=-\left(k-1\right)/2}}\stackrel{\left(k-1\right)/2}{\sum _{j=-\left(k-1\right)/2}}B\left(\alpha ,z+h,x+j\right)\cdot \left(h,j\right)$$[4]

where a square window with k pixels is used. In this paper, k is set to 35 pixels for an image with total 300 * 300 pixels. h and j are the pixel indexes inside the window. The cross-sectional vector maps $\overrightarrow{v}\left(\alpha ,z,x\right)$ reflect the local structure of the binary image at the position (z,x) (see Figure 2C). The vector cross-sections within each class were paired. The Euclidear distances of the difference between the paired vector maps were performed to characterize the local structure difference between two binary images (see Figure 2D), and they were further summed over the full z-x dimensions with a step of k as the global difference:

$$\Delta v\left(m,l\right)=\stackrel{Z}{\sum _{z=1}}\stackrel{X}{\sum _{x=1}}\left|\overrightarrow{v}\left(m,z,x\right)-\overrightarrow{v}\left(l,z,x\right)\right|$$[5]

where $|·|$ denotes Euclidear distance, Z and X are the number of windows in the depth and width directions of the image, respectively. m and l are the α values of the paired vector maps $\overrightarrow{v}\left(\alpha ,z,x\right)$. The binary images with similar structure (i.e., intra-class, see Figure 2, D1) would have a smaller ∆v score than the images with distinct structures (i.e., inter-class, see Figure 2, D2,D3). The BISIM of each class (V) was defined as the summation of all the paired images within each class:

$$V_i\text{=}{\displaystyle \sum }_{m,l}\Delta v\left(m,l\right),m,l\in classiandm\ne l$$[6]

where i is an index of the class,

$$BISIM=\frac{\stackrel{3}{\sum _{i=1}}{V}_i}{\stackrel{3}{\sum _{i=1}}{C}_{n_i}^2}$$[7]

where n_i is the number of vector maps in class i. The change of the BISIM value was plotted versus the angles α₁ and α₂ (see Figure 2E,F). When BISIM has the minimum value, which means obtained angles α₁ and α₂ divide signals in the ID space into three groups with minimum structural difference within each group, so the corresponding angle α₁ was assigned to the threshold α_T (see Figure 2F). In this work, α₁ and α₂ were initially set as 15˚ and 30˚, respectively. The gradient descent method (28) was used to quickly determine the minimum value of BISIM and the corresponding coordinates (${\alpha }_1^{*}$,${\alpha }_2^{*}$). Accordingly, the binarized cross-sectional angiogram was generated by applying the ID-threshold with ${\alpha }_T={\alpha }_1^{*}$ (see Figure 2G).

The quality of the automatic segmentation result was evaluated by its consistence with the ground truth. Each pixel of the binary map was classified as either true-positive Tp, true-negative Tn, false-positive Fp, or false-negative Fn. Based on the ratio of the amount of pixels in each set, the following quality metrics are calculated for a segmentation (20):

$$Sensitivity:Se=\frac{Tp}{Tp+Fn}$$[8]

$$Specificity:Sp=\frac{Tn}{Tn+Fp}$$[9]

$$Youden’sindex:J=Se+Sp-1$$[10]

Manual segmentations on Gaussian-filtered (sigma=0.5, kernel size=3) MIP image of OCTA were used to generate ground truth. Four manual segmentations were created by an expert for each data. Then, Simultaneous Truth and Performance Level Estimation (STAPLE) algorithm (29) was used to obtain the probabilistic estimate map of the true segmentation according to these four manual segmentations. Every pixel in the probabilistic estimate map was regarded as the foreground only if its probability was 75% and more in the probabilistic estimate map, while every other pixel was marked as background (20).

In order to facilitate comparison, conventional 2-step (intensity and binarization) thresholding method in the decorrelation-contrast OCTA was also performed as the control group, which was abbreviated as “2-step thresholding” in the following part. The similar binarization strategies were used to generate binary map for further quantitative evaluation (18,23,30) or serve as a control group (21). In control group, two kinds of threshold were used for binarization: intensity-threshold (removing the regions without sufficient SNR level and generating flow image) (22), decorrelation-threshold (i.e., binarization threshold). The intensity-threshold level was adjusted to two different levels: 6 standard deviations (2-step-6σ) or 3 standard deviations (2-step-3σ) above the mean of the estimated noise process (31). The estimation of noise was taken once per acquisition, meaning that it is possible for the threshold levels to change, though typically by a small amount (31). The decorrelation threshold was determined with the global Otsu method in the enface MIP OCT angiogram (16,18).

Vascular quantitative metrics

In clinical study, to quantitatively evaluate the differences between several groups of vasculatures, vascular quantitative metrics were calculated based on the binary maps, corresponding skeleton maps and perimeter maps (10,32,33). Conventional process to generate skeleton map and perimeter map was performed on the binary projection (10-12), however, the depth information was lost during the process of projection. Consequently, vessels in different layers might be identified as a single continuous vessel. To solve this problem, we obtain the skeleton map and perimeter map based on the 3D binary result and then acquire corresponding projection images. To generate the depth encoded projectional image, the skin surface was flattened based on the structural information and the depth was defined as the vertical distance between the vessel pixel and the skin surface boundary.

Here, several widely accepted statistical parameters were calculated, including VDI, VSD, VAD, VPI and VCI (10), the definitions are given as:

$$VDI=\frac{{\sum }_{x=1,y=1}^{n}A\left(x,y\right)}{{\sum }_{x=1,y=1}^{n}S\left(x,y\right)}$$[11]

$$VSD=\frac{{\sum }_{x=1,y=1}^{n}S\left(x,y\right)}{n^2}$$[12]

$$VAD=\frac{{\sum }_{x=1,y=1}^{n}A\left(x,y\right)}{n^2}$$[13]

$$VPI=\frac{{\sum }_{x=1,y=1}^{n}P\left(x,y\right)}{n^2}$$[14]

$$VCI=\frac{{[{\sum }_{x=1,y=1}^{n}P\left(x,y\right)]}^2}{4\pi {\sum }_{x=1,y=1}^{n}A\left(x,y\right)}$$[15]

where n represents the width and height of the enface OCTA image. A(x,y) represents the pixels registered as vessel area. S(x,y) represents the pixels registered as vessel skeleton. P(x,y) represents the pixels as vessel edge. VDI reflects the average diameter of blood vessels in the image and is sensitive to the dilation of blood vessels. VSD was calculated as the ratio of the length occupied by the blood vessels to the total area in the skeletonized vessel map. VAD was calculated as the ratio of vessel area to total binary image. Vascular length and diameter changes can cause significant changes in VAD. VPI was calculated as the ratio of the vessel perimeter to the total area of the OCTA image. VCI represents the visual complexity of graphics in an image. In tortuous and irregular patterns, the value of VCI is high.

In addition, vessel metrics map could be obtained by calculating quantitative descriptive metrics in local kernel. A down sampled block map could be generated with 30 × 30 (pixels) block size, the map was resized to original size and Gaussian filtered with a 3 × 3 (pixels) kernel. All vessel metrics maps were generated in the same fashion.

Statistical analysis

Vascular quantitative metrics mentioned above were firstly averaged for each condition and then presented as mean ± standard deviation. Paired sample t-test were implemented between normal skin and lesion site. Statistical significance is presented at two levels: *, P≤0.05; and **, P≤0.01.

Results

Binarization is a fundamental operation in OCTA quantification. The binarization performance of the proposed ID-BISIM threshold was quantitatively compared with the 2-step thresholding method (Figure 3). As shown in Figure 3A, based on the structural and angiographic images, the cross-sectional image can be decomposed into 4 components: air (green), superficial epidermis (blue, segmented from structure image), dermal blood flow (red, identified from OCTA image) and tissue around the dermis. The air (random noise) region above the epidermal surface is characterized with low intensity (large iSNR) and high decorrelation in the ID space (green dots in Figure 3B). In contrast, the superficial and avascular epidermis layer exhibits high intensity (small iSNR) and low decorrelation (no dynamic blood flow) in the ID space (the blue dots in Figure 3B). In the dermal layer, the dynamic flow (red dots in Figure 3B) and the surrounding tissue (black dots in Figure 3B) were above and below the ID-BISIM threshold (red line in Figure 3B), respectively. The ID-BISIM method enables a 3D-adaptive threshold, while the global thresholds used in the 2-step thresholding method cannot separate the dynamic signals from the static background completely. To be specific, in 2-step thresholding method, a high intensity threshold (6σ) lost amount of blood flow signals (compare Figure 3D,E), while a low intensity threshold (3σ) left numerous static background (compare Figure 3D,F), which were also indicated by the dashed lines in Figure 3B. To quantitatively demonstrate the advantages of ID-BISIM over 2-step thresholding method, a total of 40 scans from ten subjects were used for comparing difference methods with manual segmentations (Figure 3C). The result was summarized in Table 1. Generally speaking, the proposed ID-BISIM threshold enables a 3D adaptive binarization with an improved sensitivity and specificity over the conventional 2-step thresholding method, corresponding to a 37–65% improvement of Youden’s index. It takes about 1.15s to get α_T for one data, which was tested using MATLAB (R2019a, MathWorks, Natick, Massachusetts) running on a 64-bit PC with Core inter(R) (TM)i7-9700k CPU, and 16GB RAM.

Figure 3 ID-BISIM threshold enables a 3D adaptive binarization of OCTA images with improved sensitivity and specificity. (A) Representative cross-sectional OCTA image overlaid with the corresponding structural (intensity) image. Green region: the air above the skin surface, pure noise. Blue region: the superficial, avascular epidermal layer segmented based on OCT structural image, pure static tissue. Red region: the dermal blood flow identified based on ID-BISIM, pure dynamic flow. (B) ID space distribution of the pixels of the air (green), epidermis (blue), dermal blood flow (red) and dermal surrounding tissue (black). The red line denotes the demarcation line obtained by ID-BISIM method. The dashed purple horizontal lines represent the intensity threshold of the 2-step thresholding method (mean value add 3σ or 6σ). σ is the standard deviation of the estimated noise process. The dashed purple vertical lines represent the decorrelation threshold. (C) Manual ground-truth segmentation. Binarized OCT angiograms with (D) ID-BISIM, (E) 2-step thresholding-6σ, and (F) 2-step thresholding-3σ thresholds. The three rows are 1) Binarized cross-sectional angiograms, 2) binarized enface angiograms (red dashed line indicates the position of the cross-section in 1), and 3) binarization image overlaid on manual segmentation mask (white: Tp, black: Tn, red: Fn, green: Fp), respectively. OCTA, optical coherence tomography angiography; ID-BISIM, ID-binary image similarity.

Table 1 Comparison of ID-BISIM and 2-step thresholding methods using manual segmentations as ground truth
Full table

Skeleton operation is performed on the binarization images for the quantification of VSD and VPI. The proposed ID-BISIM threshold enabled a 3D binarized vessel image (see Figure 4A). Rather than the projectional 2D binarized image (see Figure 4B), the additional depth information facilitated the differentiation of the vessels which are close to each other in the enface plane but at different depths. For example, the two vessels of different depths (marked with the blue and pink arrows) were successfully differentiated (see Figure 4C) in 3D binarization, while in 2D binarized image (see Figure 4D) they were presented as a single large vessel.

Figure 4 Comparison of 3D and 2D binarization in vessel skeleton. Binarized (A) 3D and (B) 2D vascular images using ID-BISIM threshold, i.e., the vessel area maps. Depth is encoded with color. (C) 3D and (D) 2D vessel skeleton maps. The blue and pink arrows indicate two vessels that are close to each other in the enface plane but at different depths.

Quantitative OCTA image can be readily achieved with the proposed ID-BISIM threshold. Apparent changes of the vessel morphology were observed between the normal skin (Figure 5A) and the PWS lesion site (Figure 5B). Vascular dilation and hyperplasia were observed in PWS (see VDI and VAD in Figure 5), which corresponding with the clinical symptoms of bright red patches.

Figure 5 Quantitative OCTA of skin vascular morphology. 1) Original enface ID-OCTA images. 2) Vessel diameter map. 3) Vessel skeleton map. 4) Vessel area map. 5) Vessel perimeter map. 6) Vessel complexity map. All the maps are integrated with the OCTA images. (A) Normal skin (the symmetrical normal site of the lesion on cheek). (B) PWS lesion skin. OCTA, optical coherence tomography angiography.

Statistical results of the vascular quantitative metrics for normal skin and PWS are shown in Table 2. Using the proposed ID-BISIM thresholding method, the VDI and VAD increased significantly in PWS (P=0.0015 and 0.0485 respectively). Besides, significant decrease was observed in VCI with P=0.0094, while there are no statistical significances in VSD and VPI. Similar changes were also observed in the 2-step-3σ thresholding method, but presented less statistically significant (only VDI was statistically significant with P=0.024).

Table 2 Quantitative analysis of vascular metrics for normal skin and PWS
Full table

Discussion

Vascular quantitative metrics have been widely used in scientific research and clinical application, which require accurate binarization and effective extraction of vessel skeleton (13,34-36). In order to obtain precise binarization results, we proposed a 3D adaptive ID threshold using the linear relationship between the local intensity and complex-decorrelation, and the ID threshold was automatically determined based on the similarity of binary image structure. Furthermore, we extend the skeleton extraction operator to 3D space so as to distinguish vessels in different depths. Finally, the validity of the proposed ID-BISIM method was demonstrated by comparing with the conventional 2-step thresholding method. Thus, the superior performance of the proposed ID-BISIM method will be helpful to the quantitative diagnosis and treatment of skin diseases.

Distinct from the conservative 3σ criterion used in the grayscale ID-OCTA, the proposed ID-BISIM works for the purpose of determining an optimal ID threshold to assign all OCTA image pixels to be either black (background) or white (vessel) with minimal segmentation errors. Based on the derived linear relation of the decorrelation to the iSNR in reference (9), the BISIM evaluates the intra-class similarity of the three groups divided by ID angles α₁ and α₂, and ID-BISIM was used to automatic determine the ID-thresholds with minimal BISIM, or in other words, with minimal segmentation error. The proposed ID-BISIM algorithm offers several advantages over existing binarization methods. Rather than simply removed low SNR region in 2-step thresholding method, which means a balance must be struck between including valid signals and excluding noise (31), the ID-BISIM threshold enables a 3D adaptive binarization with an improved sensitivity and specificity over the conventional 2-step method, and presented improved statistical significance in the differentiation of PWS from the health skin.

The ID-threshold obtained with ID-BISIM was mainly influenced by the diameter and area density of vessels. To be specific, if the tomogram is dominant by vessels, the dynamic signals are densely distributed in ID space, while the static signals are distributed sparsely, under the minimum classification error criterion, the corresponding ID threshold obtained by ID-BISIM was larger. This adaptive threshold is essential for the accuracy of binarization and further quantitative analysis. Besides, different OCT systems hardly influence the ID threshold because the intensity of OCT signals has been normalized by the noise level and the distribution of OCT signals in ID space will only have a global shift along the asymptotic distribution. However, although the threshold is almost constant for different systems, the SNR of OCT systems will still influence the final angiograms due to the overall displacement of OCT signals in ID space.

In addition, the window size k of the morphologic vector array is of importance for the performance of ID-BISIM. To be specific, if k is set too large, the morphologic vector array $\overrightarrow{v}\left(\alpha ,z,x\right)$ will lose its sensitivity, on the other hand, if k is set too small, the stability of the vector array is degraded. In both cases, the vascular structures cannot be effectively extracted. Here, as a rule of thumb, we recommend setting k to be 2~3 times the capillary diameter.

In conclusion, the proposed ID-BISIM enables an automatic 3D adaptive vessel segmentation with enhanced performance in quantitative OCTA. The vascular quantitative metrics would be a useful tool for improving the diagnosis and the treatment of PWS.

Acknowledgments

Funding: This work is supported by National Natural Science Foundation of China (62075189), Zhejiang Provincial Natural Science Foundation of China (LR19F050002), Zhejiang Lab (2018EB0ZX01), Fundamental Research Funds for the Central Universities (2018FZA5003).

Footnote

Provenance and Peer Review: With the arrangement by the Guest Editors and the editorial office, this article has been reviewed by external peers.

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at http://dx.doi.org/10.21037/qims-20-868). The special issue “Advanced Optical Imaging in Biomedicine” was commissioned by the editorial office without any funding or sponsorship. Peng Li reports grants from National Natural Science Foundation of China, Zhejiang Provincial Natural Science Foundation of China, Zhejiang Lab, and Fundamental Research Funds for the Central Universities, during the conduct of the study. The authors have no other conflicts of interest to declare.

Ethical Statement: This study was approved by Ethics Committee of Chinese PLA General Hospital (NO.: 0404) and informed consent was taken all the patients.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Smith TK, Choi B, Ramirez-San-Juan JC, Nelson JS, Osann K, Kelly KM. Microvascular blood flow dynamics associated with photodynamic therapy, pulsed dye laser irradiation and combined regimens. Lasers Surg Med 2006;38:532-9. [Crossref] [PubMed]
2. Babilas P, Shafirstein G, Baumler W, Baier J, Landthaler M, Szeimies RM, Abels C. Selective photothermolysis of blood vessels following flashlamp-pumped pulsed dye laser irradiation: in vivo results and mathematical modelling are in agreement. J Invest Dermatol 2005;125:343-52. [Crossref] [PubMed]
3. Chen JK, Ghasri P, Aguilar G, van Drooge AM, Wolkerstorfer A, Kelly KM, Heger M. An overview of clinical and experimental treatment modalities for port wine stains. J Am Acad Dermatol 2012;67:289-304. [Crossref] [PubMed]
4. Wang Y, Gu Y, Liao X, Chen R, Ding H. Fluorescence monitoring of a photosensitizer and prediction of the therapeutic effect of photodynamic therapy for port wine stains. Exp Biol Med (Maywood) 2010;235:175-80. [Crossref] [PubMed]
5. Liu G, Jia W, Nelson JS, Chen Z. In vivo, high-resolution, three-dimensional imaging of port wine stain microvasculature in human skin. Lasers Surg Med 2013;45:628-32. [Crossref] [PubMed]
6. Zhao Y, Brecke KM, Ren H, Ding Z, Nelson JS, Chen Z. Three-dimensional reconstruction of in vivo blood vessels in human skin using phase-resolved optical Doppler tomography. IEEE J Sel Top Quantum Electron 2001;7:931-5. [Crossref]
7. Liu Y, Zhu D, Xu J, Wang Y, Feng W, Chen D, Li Y, Liu H, Guo X, Qiu H, Gu Y. Penetration-enhanced optical coherence tomography angiography with optical clearing agent for clinical evaluation of human skin. Photodiagnosis Photodyn Ther 2020;30:101734. [Crossref] [PubMed]
8. Cheng Y, Guo L, Pan C, Lu T, Hong T, Ding Z, Li P. Statistical analysis of motion contrast in optical coherence tomography angiography. J Biomed Opt 2015;20:116004. [Crossref] [PubMed]
9. Huang L, Fu Y, Chen R, Yang S, Qiu H, Wu X, Zhao S, Gu Y, Li P. SNR-Adaptive OCT Angiography Enabled by Statistical Characterization of Intensity and Decorrelation With Multi-Variate Time Series Model. IEEE Trans Med Imaging 2019;38:2695-704. [Crossref] [PubMed]
10. Chu Z, Lin J, Gao C, Xin C, Zhang Q, Chen CL, Roisman L, Gregori G, Rosenfeld PJ, Wang RK. Quantitative assessment of the retinal microvasculature using optical coherence tomography angiography. J Biomed Opt 2016;21:66008. [Crossref] [PubMed]
11. Alam M, Thapa D, Lim JI, Cao D, Yao X. Computer-aided classification of sickle cell retinopathy using quantitative features in optical coherence tomography angiography. Biomed Opt Express 2017;8:4206-16. [Crossref] [PubMed]
12. Reif R, Qin J, An L, Zhi Z, Dziennis S, Wang R. Quantifying optical microangiography images obtained from a spectral domain optical coherence tomography system. Int J Biomed Imaging 2012;2012:509783. [Crossref] [PubMed]
13. Liew YM, McLaughlin RA, Gong P, Wood FM, Sampson DD. In vivo assessment of human burn scars through automated quantification of vascularity using optical coherence tomography. J Biomed Opt 2013;18:061213. [Crossref] [PubMed]
14. Lozzi A, Agrawal A, Boretsky A, Welle CG, Hammer DX. Image quality metrics for optical coherence angiography. Biomed Opt Express 2015;6:2435-47. [Crossref] [PubMed]
15. Men SJ, Chen CL, Wei W, Lai TY, Song SZ, Wang RK. Repeatability of vessel density measurement in human skin by OCT-based microangiography. Skin Res Technol 2017;23:607-12. [Crossref] [PubMed]
16. Otsu N. A Threshold Selection Method from Gray-Level Histograms. IEEE Trans Syst Man Cybern Syst 1979;9:62-6. [Crossref]
17. Mehta N, Liu K, Alibhai AY, Gendelman I, Braun PX, Ishibazawa A, Sorour O, Duker JS, Waheed NK. Impact of Binarization Thresholding and Brightness/Contrast Adjustment Methodology on Optical Coherence Tomography Angiography Image Quantification. Am J Ophthalmol 2019;205:54-65. [Crossref] [PubMed]
18. Devarajan K, Di Lee W, Ong HS, Lwin NC, Chua J, Schmetterer L, Mehta JS, Ang M. Vessel density and En-face segmentation of optical coherence tomography angiography to analyse corneal vascularisation in an animal model. Eye Vis (Lond) 2019;6:2. [Crossref] [PubMed]
19. Schottenhamml J, Moult EM, Ploner S, Lee B, Novais EA, Cole E, Dang S, Lu CD, Husvogt L, Waheed NK, Duker JS, Hornegger J, Fujimoto JG. An Automatic, Intercapillary Area-Based Algorithm for Quantifying Diabetes-Related Capillary Dropout Using Optical Coherence Tomography Angiography. Retina 2016;36 Suppl 1:S93-S101. [Crossref] [PubMed]
20. Casper M, Schulz-Hildebrandt H, Evers M, Birngruber R, Manstein D, Huttmann G. Optimization-based vessel segmentation pipeline for robust quantification of capillary networks in skin with optical coherence tomography angiography. J Biomed Opt 2019;24:1-11. [Crossref] [PubMed]
21. Gao SS, Jia Y, Liu L, Zhang M, Takusagawa HL, Morrison JC, Huang D. Compensation for Reflectance Variation in Vessel Density Quantification by Optical Coherence Tomography Angiography. Invest Ophthalmol Vis Sci 2016;57:4485-92. [Crossref] [PubMed]
22. Enfield J, Jonathan E, Leahy M. In vivo imaging of the microcirculation of the volar forearm using correlation mapping optical coherence tomography (cmOCT). Biomed Opt Express 2011;2:1184-93. [Crossref] [PubMed]
23. Jia Y, Tan O, Tokayer J, Potsaid B, Wang Y, Liu JJ, Kraus MF, Subhash H, Fujimoto JG, Hornegger J, Huang D. Split-spectrum amplitude-decorrelation angiography with optical coherence tomography. Opt Express 2012;20:4710-25. [Crossref] [PubMed]
24. Carter HH, Gong P, Kirk RW, Es'haghian S, Atkinson CL, Sampson DD, Green DJ, McLaughlin RA. Optical coherence tomography in the assessment of acute changes in cutaneous vascular diameter induced by heat stress. J Appl Physiol (1985) 2016;121:965-72. [Crossref] [PubMed]
25. Choi WJ, Reif R, Yousefi S, Wang RK. Improved microcirculation imaging of human skin in vivo using optical microangiography with a correlation mapping mask. J Biomed Opt 2014;19:36010. [Crossref] [PubMed]
26. Yang S, Liu K, Ding H, Gao H, Zheng X, Ding Z, Xu K, Li P. Longitudinal in vivo intrinsic optical imaging of cortical blood perfusion and tissue damage in focal photothrombosis stroke model. J Cereb Blood Flow Metab 2019;39:1381-93. [Crossref] [PubMed]
27. Nam AS, Chico-Calero I, Vakoc BJ. Complex differential variance algorithm for optical coherence tomography angiography. Biomed Opt Express 2014;5:3822-32. [Crossref] [PubMed]
28. Li P, Huang Z, Yang S, Liu X, Ren Q, Li P. Adaptive classifier allows enhanced flow contrast in OCT angiography using a histogram-based motion threshold and 3D Hessian analysis-based shape filtering. Opt Lett 2017;42:4816-9. [Crossref] [PubMed]
29. Warfield SK, Zou KH, Wells WM. Simultaneous truth and performance level estimation (STAPLE): an algorithm for the validation of image segmentation. IEEE Trans Med Imaging 2004;23:903-21. [Crossref] [PubMed]
30. Gong P, McLaughlin RA, Liew YM, Munro PR, Wood FM, Sampson DD. Assessment of human burn scars with optical coherence tomography by imaging the attenuation coefficient of tissue after vascular masking. J Biomed Opt 2014;19:21111. [Crossref] [PubMed]
31. Cole ED, Moult EM, Dang S, Choi W, Ploner SB, Lee B, Louzada R, Novais E, Schottenhamml J, Husvogt L, Maier A, Fujimoto JG, Waheed NK, Duker JS. The Definition, Rationale, and Effects of Thresholding in OCT Angiography. Ophthalmol Retina 2017;1:435-47. [Crossref] [PubMed]
32. Jia Y, Bailey ST, Hwang TS, McClintic SM, Gao SS, Pennesi ME, Flaxel CJ, Lauer AK, Wilson DJ, Hornegger J, Fujimoto JG, Huang D. Quantitative optical coherence tomography angiography of vascular abnormalities in the living human eye. Proc Natl Acad Sci U S A 2015;112:E2395-402. [Crossref] [PubMed]
33. Dongye C, Zhang M, Hwang TS, Wang J, Gao SS, Liu L, Huang D, Wilson DJ, Jia Y. Automated detection of dilated capillaries on optical coherence tomography angiography. Biomed Opt Express 2017;8:1101-9. [Crossref] [PubMed]
34. Choi WJ, Wang H, Wang RK. Optical coherence tomography microangiography for monitoring the response of vascular perfusion to external pressure on human skin tissue. J Biomed Opt 2014;19:056003. [Crossref] [PubMed]
35. Deegan AJ, Talebi-Liasi F, Song S, Li Y, Xu J, Men S, Shinohara MM, Flowers ME, Lee SJ, Wang RK. Optical coherence tomography angiography of normal skin and inflammatory dermatologic conditions. Lasers Surg Med 2018;50:183-93. [Crossref] [PubMed]
36. Babalola O, Mamalis A, Lev-Tov H, Jagdeo J. Optical coherence tomography (OCT) of collagen in normal skin and skin fibrosis. Arch Dermatol Res 2014;306:1-9. [Crossref] [PubMed]