2.1.

Preprocessing

The preprocessing step contains two main processes, namely, first, increasing the contrast of the input image and second, computing the foreground and background marker images.

2.1.1.

Increasing the contrast of the input image

Noisy images have negative effects on the watershed algorithm and cause oversegmentation. In this study, we used the Optovue OCTA device for imaging. Using this device, enface images can be extracted in color and grayscale modalities. In this step, first, the red channel of the input image was extracted because, as is shown in Fig. 3, this channel has the highest contrast between the vessel pixels and background in comparison with the gray-scale image, red and blue channels provided by the Optovue system. Then, we applied two-scale ($\sigma =1$ and $\sigma =2\text{pixels}$) vessel enhancement filtering to remove background noise from the red channel of the input image and enhance the contrast between the vessel and background pixels.²⁸ The vesselness measure $V_0$ obtained by applying the vesselness filter on the red channel of a superficial image is shown in Fig. 4(b). As it is obvious, the contrast between the flow pixels and background has been increased.

Fig. 3

A visual comparison among the contrast of different variants of the input image (the provided grayscale image and three RGB channels). (a) RGB image, (b) the provided grayscale image, (c) red channel, (d) green channel, and (e) blue channel.

Fig. 4

Preprocessing step. (a) Red channel in Fig. 3(c), (b) Enhanced image (vesselness measure), (c) background marker image, (d) distance map image, and (e) foreground marker image.

2.2.

Marker-Controlled Watershed Segmentation

The watershed algorithm by immersion was used to segment abnormal CNP region candidates. In the following, we first present a brief description of the watershed algorithm by immersion and then we describe how it has been adjusted in the form of the marker-controlled watershed to segment CNP regions up to their true boundaries.

2.2.2.

Applying watershed algorithm to segment CNP region candidates

Generally, watershed algorithm is applied to the gradient image³⁷ because in the gradient image, objects appear as dark regions and their boundaries are located at pixels with high gradient values. The Sobel operator is one of the well-known and simple filters in terms of computation for obtaining gradient magnitude of the image. Using this filter, the gradient magnitude of the image $f(x,y)$ is computed as follows:⁴⁰

Eq. (4)

$$f_g(x,y)=\sqrt{{[G_x*f(x,y)]}^2+{[G_y*f(x,y)]}^2},$$

where $G_x$ and $G_y$ are Sobel edge masks in $X$ and $Y$ directions, respectively. Figure 5(a) shows the gradient magnitude of Fig. 4(b).

Fig. 5

The output of applying the marker-controlled watershed to the gradient magnitude of Fig. 3(b). (a) The gradient magnitude of Fig. 3(b), (b) applying the conventional watershed algorithm to (a), (c) imposing marker image [Figs. 4(c) and 4(e)] in (a), and (c) the marker-controlled watershed result.

Due to the existence of numerous regional minima in the gradient image, directly applying watershed algorithm to the gradient image caused oversegmentation. Figure 5(b) shows the result of applying the watershed algorithm to the gradient magnitude image in Fig. 5(a). To overcome this intrinsic drawback of the watershed algorithm, the marker-controlled watershed algorithm was used. Using this technique, the number of regional minima was controlled by defining marker image and thus flooding only occurred in basins corresponding to each marker.⁴¹ Therefore, in the marker-controlled watershed, first, the marker points should be computed and then these markers should be superimposed on the gradient image as global minimums.

Based on the foreground marker image $f_{fgm}(x,y)$ and the background marker image $f_{bgm}(x,y)$ provided by the preprocessing step, the marker image $f_m(x,y)$ is defined as follows:

Eq. (5)

$$f_m(x,y)=f_{fgm}(x,y)\cup f_{bgm}(x,y).$$

Then, minima imposition morphological operation was used to superimpose each marker on the gradient image.⁴⁰ In this operation, as it is defined by Eq. (6), first, the pointwise minimum between the gradient image $f_g$ and the marker image $f_m$ is computed ($f_g\wedge f_m$), then using morphological reconstruction by erosion, the minimum superimposed gradient image $f_{msg}(x,y)$ is obtained [Fig. 5(c)]:

Eq. (6)

$$f_{msg}(x,y)=R_{f_g\wedge f_m}^{ϵ}[f_m(x,y)],$$

where $R^{ϵ}$ is the morphological reconstruction by erosion of $f_g\wedge f_m$ from the marker image $f_m$. In fact, the pointwise minimum forces each region in the gradient image that is corresponded with each white region in the marker image to be a global minimum of the gradient image. Then, all gradient image information except the masked regions (global minimums) are recovered using the morphological reconstruction by erosion.

Finally, CNP region candidates were obtained by applying the watershed algorithm to the superimposed gradient image $f_{msg}(x,y)$. Figure 5(d) shows the result of applying the watershed algorithm to $f_{msg}$.

2.3.

Postprocessing

As shown in Fig. 5(d), some of the segmented CNP regions are fragmented and do not show an accurate range of abnormal region. In addition, some of the normal regions were segmented along with abnormalities. Therefore, the postprocessing step was considered to modify the marker-controlled watershed segmentation output and obtain correct abnormal CNP region with more accurate boundary. The postprocessing step was performed by merging fragmented candidates based on the proper similarity measure and then separating correct abnormal CNP regions from modified candidates.

To merge the fragmented regions, the abnormal CNP region candidates were modeled as an undirected weighted graph $G=(V,E,w)$ so that by defining $V={\{R_i\}}_{i=1}^n$, each vertex represents a candidate region $R_i$ and each edge $e_i=(R_i,R_j)$ in $E={\{e_i\}}_{i=1}^m$ connects two neighboring candidates. $n$ and $m$ are equal to the number of candidates and the number of neighbors, respectively. By considering the weight of each edge equal to the size of the common boundary between two neighboring regions, the weighting function $w:E\to \mathbb{N}$ for each two neighboring vertexes $R_i$ and $R_j$ is defined as follows:

Therefore, the merging process was accomplished for each subgraph $G_i^{\prime }$ independently. Two neighboring regions $R_i$ and $R_j$ are similar if their common boundary is larger than threshold value $T_3$. According to Sec. 2.4, this threshold value was computed equal to 6 pixels. Thus, the largest edge in each subgraph is extracted and the threshold condition is checked. If the condition is satisfied, the nodes associated with this edge merge together. It means that the same label is given to the corresponding regions of those two nodes and the boundary between them is removed. Then, the subgraphs are updated and rebuilt. Considering $R_i$ as a region, it should be removed and merged with $R_j$, the components $V^{\prime }$, $E^{\prime }$, and $w^{\prime }$ of each subgraph $G^{\prime }$ are updated as follows:

The merging process termination condition is defined as reaching unified regions or missing thresholding condition. Figures 6(a)–6(f) show the merging process on a sample of fragmented CNP candidates. As it is seen, first two yellow and dark green regions in Fig. 6(a), which have the largest common boundary, are merged in Fig. 6(b). By this merging, the dark green region is included in the yellow region so that the boundary of the yellow region changed and the corresponding weighted graph is updated. By continuing this procedure, the unified yellow region in Fig. 6(f) is obtained. Figure 7(a) shows the result of applying the graph-based merging algorithm to the marker-controlled watershed segmentation output in Fig. 5(d).

Fig. 6

(a–f) Representation of the graph-based merging process on a sample of fragmented abnormal CNP candidates.

Fig. 7

Postprocessing step. (a) Merging fragmented abnormal CNP candidates of Fig. 5(d), (b) separating correct abnormal CNP regions from (a), and (c) mapping abnormal CNP regions on the color enhanced input image.

Finally, to remove normal regions that are wrongly extracted along with the abnormalities, a circle-checking step with a circle of radius $R$ is considered. So, correct abnormal CNP regions are obtained by removing any CNP candidates where it was not possible to draw a circle of radius $R$ inside them. A circle of radius $R$ can be inscribed in each abnormal region because of using Euclidean distance function and the way of computing VD (the distance of each pixel from the nearest nonzero pixel). The radius $R$ was considered as 4 pixels. It is computed by Eq. (12) that will be described in Sec. 2.4.

Figures 7(b) and 7(c) show the final abnormal CNP regions and their mapping on the color enhanced input image, respectively. The color enhanced input image was obtained by applying vesselness filter to each red, blue, and green channels independently.

2.4.

Determination of the Threshold Values

According to the analysis of VD, values in the OCTA images of normal and diabetic subjects, the threshold values $T_2$ in the preprocessing step, $T_3$ in the merging algorithm, and the radius $R$ of the circle in the final step are computed as follows:

Table 1

Illustration of the parameter estimation in the dataset.

The radius $R$ equal to 4 pixels ($40\mu \mathrm{m}$) is a discriminator between normal and abnormal CNP regions. Due to the noise and functionality of the watershed algorithm, some of the abnormal CNP regions are not extracted up to their true boundary. Therefore, using threshold value $T_2$ less than $R$ ($T_2=3\le R$), the regions with less VD value are extracted. Finally, in the merging step, if these regions have a common border size larger than or equal to $T_3=2T_2=6\text{pixels}$ with abnormal CNP regions, they will be considered as a part of abnormal regions and combined with them. Otherwise, they are considered as a normal region and removed. In fact, T2 is an approximation of the CNP region radius and T3 is an approximation of the CNP region diameter.

2.5.

Dataset

A set of 36 images consisting of 18 pathologic and 18 normal images were used. Pathologic and normal images were acquired from 12 diabetic participants with DR at different stages and 18 healthy volunteers, respectively.

In the patients group, the mean age was $52.82\pm 11.31$ years and female/male ratio was 4/8. In the normal group, the mean age was $39.07\pm 8.27$ and female/male ratio was 9/9. The signal strength index ranged from 46 to 69 and 62 to 83 in patient and normal group, respectively.

Images were acquired using a commercial spectral-domain OCT imaging device (RTVue-XR, Optovue, USA). This instrument uses a light source centered on 840 nm with a full-width half-maximum bandwidth of 45 nm and has an A-scan rate of 70 kHz. The tissue axial resolution of this device is $5\mu \mathrm{m}$, and the transverse is $15\mu \mathrm{m}$.

Imaging was performed over $3\times 3\mathrm{mm}$ field-of-view centered on the fovea with a 1.6-mm depth to obtain OCTA volumes consisting of $304\times 304\times 512\text{voxels}$. Each data cube was acquired in $\sim 2.9\mathrm{s}$. Two repeated B-scans, each containing 304 A-scans, were captured at each fixed position before proceeding to the next sampling location. A three dimensional data cube was formed by orthogonal registration and merging two consecutive scan volumes including one $x$-fast scan and one $y$-fast scan.⁴²

Angiography information was extracted using the split-spectrum amplitude-decorrelation angiography algorithm by computing the signal amplitude-decorrelation between two consecutive B-scans of the same location.⁴³

Using the built-in software, the OCTA volume is segmented to different structural boundaries. Then, by computing maximum flow projection within the slabs defined by the segmented boundaries, en-face projection angiograms, including superficial [between $3\mu \mathrm{m}$ beneath the internal limiting membrane and $15\mu \mathrm{m}$ beneath the inner plexiform layer (IPL)], deep (between 15 and $70\mu \mathrm{m}$ beneath the IPL), and choriocapillaris (between 31 and $59\mu \mathrm{m}$ beneath the retinal pigment epithelium), are constructed.^43,44 Since DR influences superficial and deep capillary plexuses, we segmented CNP regions in these two angiograms.

In Optovue OCTA devices, en-face images can be directly extracted in color and gray-level modalities. In this work, we used the red channel of the color modality because it has higher contrast compared to the corresponding grayscale image.

3.1.

Evaluation Process

We evaluated our proposed method both qualitatively and quantitatively. The segmentation results were compared with manual delineation ground truth provided by two trained ophthalmologists (observer 1 and observer 2). Using painting tool of the Photoshop software, the experts manually drew the boundary of each abnormal CNP regions on the color en-face angiogram. Then, using a code we wrote in MATLAB, each color en-face image with manual delineated boundary was converted to the binary image (ground truth). The experts delineated only the boundary of abnormal nonperfusion regions, not all nonperfusion areas.

CNP regions were manually segmented twice by each observer. The variability in the manual segmentations delineated by different observers (interobserver agreement) and by the same observer (intraobserver agreement) were evaluated using Dice similarity coefficient (DSC), which is computed as follows:⁴⁵

In intraobserver variability evaluation, each observer segmented the CNP regions manually in two different sessions. For both observers, the time interval of segmentation between these two sessions was more than 1 month in the superficial plexuses and more than 3 weeks in the deep angiograms.

Table 2 shows the $\text{mean}\pm \text{standard deviation}$ of DSC for the evaluation of intra- and interobserver agreement in the manual segmentation of CNP regions.

Table 2

Mean±standard deviation of DSC for intra- and interobserver variability evaluation in the manual segmentation of abnormal CNP regions for three groups, including normal, diabetic, and all (normal + diabetic) subjects in the superficial and deep plexuses.

In this table, $\mathrm{Obs}{1}_i$ and $\mathrm{Obs}{2}_i$ show the manual ground truth provided by the first and second observer in the $i$’th session, respectively. $\mathrm{Obs}{1}_{1\&2}$ and $\mathrm{Obs}{2}_{1\&2}$ are the unions of the both ground truths provided by each of the first and second observers at different sessions, respectively. As can be seen, there is a good intra- and interobserver agreement in our study (a DSC value larger than 0.7 indicates good agreement⁴⁶).

We used precision and recall measures to evaluate our proposed method quantitatively. These measures are defined as follows:

The correlation of segmentation was evaluated using Jaccard index ($J$). Considering the segmented result (SR) obtained by the algorithm and ground truth given with the dataset, $J$ is defined as follows:

The result of the proposed method was compared with the proposed method by Zhang et al.,⁹ which has recently been published for automatic CNP region segmentation and the baseline approach used in the Schottenhamml et al.¹⁰ method. The comparison was done on our dataset because there is no public available dataset in this field.

In the implementation of Zhang et al.⁹ method, we used the preprocessing steps of the proposed method instead of a reflectance-adjusted thresholding step presented in Ref. 9, because we have not accessed the B-scans in the commercial version of the Optovue device. Therefore, in the preprocessing step, first, vesselness filter was applied to the input image and then the binary vessel mask was obtained using Otsu thresholding. The main and postprocessing steps were implemented based on their published paper. In the main processing, first, a VD map was acquired by applying Euclidean distance transform to the binary vessel mask and then the CNP map was obtained by applying the threshold value $DT=4$ to the VD map. In the postprocessing step, the extracted CNP map was modified using morphological operations. So that, first, the CNP regions were eroded by a five-pixel wide square structuring element and then regions smaller than eight pixel or regions with minor axis length smaller than two pixels were eliminated. Finally, remained regions were dilated by a seven-pixel wide square structuring element. Because of this alternation in the preprocessing step, in the following, this method is called modified Zhang et al.⁹

Another simple and valuable approach, which was used in the Schottenhamml et al.,¹⁰ is the baseline approach. As it was described in Sec. 1, in the baseline approach, first the vessel network is segmented, and then each nonvessel pixel in the logical vessel map is considered as a CNP region.

To have a valid comparison with the baseline approach and denoting the usefulness of watershed algorithm, we used the vessel segmentation map provided by our proposed method using Otsu thresholding rather than the vessel segmentation method used in the Schottenhamml et al.¹⁰ method. We also used the final circle-checking step to specify abnormal CNP regions for the baseline method. Therefore, the baseline method was implemented in four steps: (1) appling vesselness filter, (2) providing the vessel binary mask using Otsu thresholding, (3) finding the connected components of the inverted binary mask, and (4) obtaining final CNP map by removing any CNP candidates that it was not possible to draw a circle of radius $R=4$ pixels inside them (circle checking).

3.2.

Evaluation of Recall and Precision Rates for the Proposed Method and the Compared Methods

Table 3 shows the results of precision and recall on all images for our proposed method, the modified Zhang et al.⁹ and the baseline methods by observer 1 ($\mathrm{Obs}{1}_{1\&2}$), observer 2 ($\mathrm{Obs}{2}_{1\&2}$), and the intersection of observer 1 and observer 2 ($\mathrm{Obs}{1}_{1\&2}\cap \mathrm{Obs}{2}_{1\&2}$) in each of the superficial and deep capillary plexuses. As can be seen, in the superficial (deep) capillary plexus, the proposed method has better performance compared to other methods with recall rate 0.6864 (0.7068) for $\mathrm{Obs}{1}_{1\&2}$, 0.7108 (0.7196) for $\mathrm{Obs}{2}_{1\&2}$, and 0.7435 (0.8045) for ${\mathrm{Obs}1}_{1\&2}\cap \mathrm{Obs}{2}_{1\&2}$ while preserving the high precision rate 0.8522 (0.7854) for $\mathrm{Obs}{1}_{1\&2}$, 0.8471 (0.8102) for $\mathrm{Obs}{2}_{1\&2}$, and 0.8240 (0.7186) for $\mathrm{Obs}{1}_{1\&2}\cap \mathrm{Obs}{2}_{1\&2}$.

Table 3

The result of precision and recall on all images for modified Zhang et al.,9 baseline, and proposed method by Obs11&2, Obs21&2, and the intersection of them in each of the superficial and deep plexuses.

Considering $\mathrm{Obs}{1}_{1\&2}\cap \mathrm{Obs}{2}_{1\&2}$ as a ground truth, in the superficial (deep) plexus, the proposed method outperformed the compared methods by improving recall rate around 44% (34%) compared to the modified Zhang et al.⁹ and precision rate around 44% (37%) compared to the baseline method.

Tables 4 and 5 show the precision and recall rates, respectively, in normal and diabetic subjects for the proposed method and the compared methods in each of the superficial and deep plexuses, individually. As can be seen, in the normal subjects, the proposed method outperforms both the modified Zhang et al.⁹ and the baseline methods by decreasing the false positive rate (high precision rate) while preserving recall rate on a high value. In the superficial (deep) plexus of the normal subjects, precision rate was improved around 49% (21%) compared to the baseline method and recall rate was improved around 30% (27%) compared to the modified Zhang et al.⁹ In the superficial (deep) plexus of the diabetic subjects, compared to the baseline method, precision rate was improved around 42% (36%) and compared to the modified Zhang et al. method⁹ recall rate was improved around 47% (36%).

Table 4

The result of precision and recall rates in normal subjects for modified Zhang et al.,9 baseline, and proposed method by Obs11&2, Obs21&2, and the intersection of them in each of the superficial and deep plexuses.

Table 5

The result of precision and recall rates in diabetic subjects for modified Zhang et al.,9 baseline, and proposed method by Obs11&2, Obs21&2 and the intersection of them in each of the superficial and deep plexuses.

Although it can be observed in Table 5 for diabetic subjects that the precision of the proposed method by the modified Zhang et al.⁹ is higher (around 14% in the superficial and 18% in the deep plexus) compared to our proposed method, the significant difference around 47% in recall rate shows the drawback of this method in correctly segmenting the regions up to the boundary. And also, it can be observed in the baseline method that the recall rate is higher for both normal (around 10% in the superficial and 12% in the deep plexus) and diabetic (around 22% in the superficial and 16% in the deep plexus) subjects, although there is a significant difference in precision rate for both normal and diabetic groups compared to the proposed method.

Figures 8 and 9 show the segmentation results of the proposed method and the compared methods using the ground truth of $\mathrm{Obs}{1}_{1\&2}\cap \mathrm{Obs}{2}_{1\&2}$ on two diabetic (Fig. 8) and two normal samples (Fig. 9) of the superficial and deep plexuses. In these figures, the green and dark blue regions show the automatic segmentation results and the ground truth, respectively. The light blue regions show the overlapping of the automatic segmentation and the ground truth.

Fig. 8

A visual comparison between our proposed method and other methods on two samples of diabetic images in each of the (a,b) superficial and (c,d) deep plexus. (a–d) Original diabetic images, (e–h) ground truths, (i–l) baseline method, (m–p) modified Zhang et al.,⁹ and (q–t) proposed method.

Fig. 9

A visual comparison between our proposed method and other methods on two samples of normal images in each of the (a,b) superficial and (c,d) deep plexus. (a–d) Original diabetic images, (e–h) ground truths, (i–l) baseline method, (m–p) modified Zhang et al.,⁹ and (q–t) proposed method.

Figures 8(i)–8(l) show the result of applying the baseline method on two diabetic samples in the superficial and deep plexuses. As can be seen, this method has the highest false positive compared to the other methods due to its simple assumption that considered each nonvessel pixels as an abnormality.

The thresholding-based method proposed by Zhang et al.⁹ failed to extract most of the abnormal regions and correctly segment them up to their boundary. However, it had the least false positive regions compared to the other methods [Figs. 8(m)–8(p)].

The proposed method could segment the abnormal CNP regions up to their true boundary more accurately than the compared methods, as shown in Figs. 8(q)–8(t). However, our method segmented some false positive regions due to the extraction of physiologic CNP regions along the large vessels, low-quality OCTA images, and oversegmentation of some abnormal CNP regions.

Although the FAZ is an anatomical part of the retina, in pathological cases, destruction and enlargement of the FAZ happen. Therefore, we did not remove the FAZ area from the segmented abnormalities.

In the normal subjects (Fig. 9), the FAZ is the only avascular region that should be segmented. As it is clear in Figs. 9(i)–9(l), the baseline method segmented some physiologic regions in addition to the FAZ due to its segmentation strategy, which considered each nonvessel pixel as an abnormality. The method proposed by Zhang et al.⁹ failed to segment the FAZ up to its true boundary due to the use of thresholding segmentation technique [Figs. 9(m)–9(p)].

According to Figs. 9(q)–9(t), our proposed method could segment the FAZ more accurately than the compared methods, although it has one false positive region in Fig. 9(l) because of the low-quality of the input image.

3.3.

Evaluation of the Segmentation Correlation for the Proposed Method and the Compared Methods

To evaluate segmentation correlation, Jaccard index between the result of the proposed method and ground truths in each of the superficial and deep capillary plexuses was computed.

Table 6 shows the result of Jaccard similarity metric in each of the superficial and deep plexuses of normal, diabetic, and all subjects, for the modified Zhang et al.,⁹ the baseline and the proposed method by $\mathrm{Obs}{1}_{1\&2}\cap \mathrm{Obs}{2}_{1\&2}$ as a ground truth. As can be seen, in comparison with the other methods, the proposed method has higher Jaccard rate on all images with a mean of 0.72 and maximum correlation up to 0.86 in the superficial plexus and a mean of 0.70 with a maximum rate of 0.84 in the deep plexus.

Table 6

Mean±standard deviation of J for modified Zhang et al.,9 baseline, and proposed method by Obs11&2∩Obs21&2 in the superficial and deep plexus of three groups, including normal, diabetic, and all (normal + diabetic) subjects.

3.4.

Evaluation of the Proposed Method and the Compared Methods in Automatic Quantification of Abnormal CNPs

As it was mentioned previously, automatic quantification of abnormal CNP regions can be helpful in monitoring patient’s treatment process. In Refs. 9 and 10, the total avascular area (TAA) and the extrafoveal avascular area (EAA) in normal and pathologic groups were computed to quantify abnormalities. Since the output of the proposed method is a binary map of abnormal CNP regions, TAA was computed by counting the number of nonzero pixels in the binary map. And, same as Zhang et al.,⁹ EAA was considered as the avascular area outside the 1-mm central circle.

Table 7 shows the TAA and the EAA values in the superficial and deep plexus of the normal and diabetic subjects for the proposed method, the compared methods, and ground truth ($\mathrm{Obs}{1}_{1\&2}\cap \mathrm{Obs}{2}_{1\&2}$). As can be seen, in both capillary plexuses, for all compared methods, the TAA and EAA are significantly larger in the diabetic group compared to the normal group.

Table 7

Mean±standard deviation of TAA and EAA for modified Zhang et al.,9 baseline, ground truth (Obs11&2∩Obs21&2), and proposed method in the superficial and deep plexuses of the normal and diabetic subjects.

It can be observed from Table 7 that the proposed method has the closest values of TAA and EAA to the ground truth in each of the superficial and deep plexuses compared to the other methods. In addition, the EAA values in normal images show that the least physiologic regions (false positives) were extracted by the proposed method.

The baseline method has the maximum of TAA and EAA in each normal and diabetic group, which indicates a high false positive rate in this method. And the modified Zhang et al.⁹ method has the least TAA and EAA values due to the high false negative rate in this method.

The proposed method has been implemented in MATLAB 2016a and tested on a Laptop with Intel Core i7-4510U at 2.00 GHz and 8 GB RAM. The average processing time of the algorithm for each image was $\sim 3.31$ and 4.12 s in the superficial and deep plexuses, respectively.

4.1.

Efficiency of the Proposed Method in Automatic Segmentation of Abnormal CNP Regions

In comparison with some state-of-the-art methods, the method proposed here is more flexible and extensible. The proposed method is flexible and extensible because it is possible to extend it based on the new findings in the medical studies of OCTA images by changing the similarity and removing measures in the merging and removing algorithms. The study of CNP regions and their changes in the en-face angiograms of OCTA images is an open medical research field. For example, it is observed in the superficial and deep plexus of diabetic patients that many of the microaneurysms and vessel irregularities are happened in the vicinity of the abnormal CNP regions.^19,47 Therefore, our result can be improved by extending the proposed method based on these new features.

Another advantage of the proposed postprocessing step is its ability to remove some artifacts created by the vesselness filter in the preprocessing step. By using this filter, the speckle structure is enhanced alongside improving the contrast of the input image [Figs. 4(b) and 4(c)].

As it was shown in Figs. 7(a) and 7(b), these artifacts did not affect the final result because by using merging and removing algorithms in the postprocessing step, they are merged with the nonperfusion regions or they are removed.

As it was shown in the evaluation results, the proposed method has more performance compared to the modified Zhang et al.,⁹ and the baseline approach used in the Schottenhamml et al.¹⁰ In the proposed method by Zhang et al.,⁹ the conception of the abnormal CNP regions is almost considered by computing VD map, but finally, abnormal regions are extracted by thresholding this map. In this method, in order to decrease false positive in normal group, the considered threshold value is relatively large that causes some abnormal CNP regions not to be segmented or segmented incompletely. In the baseline approach, after vessel network extraction, all that remains is considered as abnormal CNP regions. This approach not only cannot calculate a trustable amount of the area as a quantification factor of the abnormal CNP regions but also cannot show an accurate map of their position. The reason is that wherever there is the smallest path between CNP regions through vessel network discontinuity, this approach connects the regions and represents them as a large CNP region; this path might occur due to the vessel segmentation error or it was a normal discontinuity in vessel network. Therefore, the efficiency of this approach is highly dependent on the accuracy of the vessel segmentation process while providing an accurate vessel segmentation algorithm in 2-D projected OCTA images is difficult due to the low quality of OCTA images and motion artifacts.

OCTA imaging devices are progressing toward improving signal-to-noise ratio and enhancing image quality but the software that is included in the commercial versions of this device is not properly capable of eliminating noise and managing the distortions caused by the patient movement. Although this motion artifact can be reduced by using the tracking module, imaging with this mode takes longer time than patient tolerance. Therefore, proposing a method that is less sensitive to low-quality images is valuable. The core of the proposed method (the watershed algorithm) is sensitive to noise, but we could significantly reduce this sensitivity by computing foreground and background marker images in the preprocessing step and using merging and removing procedures in the postprocessing step.

Overall, our proposed automatic abnormal CNP segmentation method has the following advantages:

4.2.

Efficiency of the Proposed Method in Automatic Quantification of Abnormal CNP Regions

Different parameters are used for quantifying capillary perfusion alternation in the macular disease that can be categorized into two groups. The first group is vessel density-based parameters used in Refs. 11, 48, and 49. The second group is nonperfusion area-based parameters used in Refs. 9, 10, 47, and 49. A study has been recently published, which shows that the nonperfusion area-based parameters may be more sensitive to microvascular impairment compared to the vessel density-based parameters.⁴⁹

All nonperfusion area-based parameters for clinical studies can be computed using the abnormal CNP map provided by our proposed method. In addition, in comparison with other nonperfusion area-based methods, our proposed method has the following advantages: