Automatic quantification of choroidal neovascularization lesion area on OCT angiography based on density cell-like P systems with active membranes

: Detecting and quantifying the size of choroidal neovascularization (CNV) is important for the diagnosis and assessment of neovascular age-related macular degeneration. Depth-resolved imaging of the retinal and choroidal vasculature by optical coherence tomography angiography (OCTA) has enabled the visualization of CNV. However, due to the prevalence of artifacts, it is difficult to segment and quantify the CNV lesion area automatically. We have previously described a saliency algorithm for CNV detection that could identify a CNV lesion area with 83% accuracy. However, this method works under the assumption that the CNV region is the most salient area for visual attention in the whole image and consequently, errors occur when this requirement is not met (e.g. when the lesion occupies a large portion of the image). Moreover, saliency image processing methods cannot extract the edges of the salient object very accurately. In this paper, we propose a novel and automatic CNV segmentation method based on an unsupervised and parallel machine learning technique named density cell-like P systems (DEC P systems). DEC P systems integrate the idea of a modified clustering algorithm into cell-like P systems. This method improved the accuracy of detection to 87.2% on 22 subjects and obtained clear boundaries of the CNV lesions.


Introduction
Choroidal neovascularization (CNV) is a manifestation of neovascular age-related macular degeneration (AMD) and a leading cause of blindness for the elderly population [1]. It consists of pathological choroidal vessels extending through Bruch's membrane into the outer retina. Optical coherence tomography angiography (OCTA) [2][3][4] is superior to dye based angiography for imaging CNV because it is noninvasive and it allows three dimensional assessment of CNV [5][6][7][8]. However, superficial retinal flow casts shadowgraphic artifact flow onto the outer retinal layers, preventing clear visualization of CNV flow. Many algorithms designed to remove these projections from outer retinal OCTA images have been reported recently [9][10][11][12]. Although these methods effectively remove most of the projection artifacts, residual noise in the en face angiograms of the outer retinal slab persists and additional image processing is needed to distinguish noise from pixels with true CNV flow. Previously we developed an algorithm inspired by notions of context-aware salience detection of objects [13] to recognize the areas formed by true vessels [14]. Although CNV is generally more "salient" compared to residual projection artifact noise, occasional misclassification has been observed. Aiming at developing a more generalizable method that detects CNV lesion area with higher accuracy and specificity, we developed a new method that identifies CNV by density-based extraction of clusters.
Cluster analysis is an unsupervised machine learning technique that recognizes objects related to each other without any prior knowledge of data sets. It has been widely used in many fields, such as image processing, pattern recognition, and bioinformatics [15,16]. Clustering is the task of finding the best partition to separate objects within the same group that share maximum similarity from objects in different groups that have minimum similarity. Cluster analysis is frequently used to segment pathologies in biomedical images [17]. A very effective modality known as density-based clustering is based on the idea that objects forming a dense region should belong to the same cluster. Density-Based Spatial Clustering of Applications with Noise (DBSCAN), OPTICS, DENCLUE and CURD are the algorithms most commonly used in this modality [18]. Among them, DBSCAN is the most popular method owing to its ability to discover arbitrary shapes without prior knowledge of the number of clusters. These characteristics make DBSCAN a suitable clustering algorithm for CNV detection.
Typically, clustering algorithms are implemented in the classical computational paradigm. However, membrane computation [19,20] is a novel paradigm that has gained popularity in the past few years for its promising features, such as the encapsulation of data, the simple representation of information, and especially its parallelism. Clustering based on P systems (a type of membrane computing) has shown excellent performance in convergence, robustness and parallelism [21][22][23][24][25][26] and has been successfully applied on the segmentation of digital images [27,28].
Inspired by the potential of membrane computation, we propose in this paper a new method based on a density cell-like P system (DEC P system) to accurately extract the CNV lesion area on outer retinal angiograms. We treat the detection of CNV as a clustering problem and the DBSCAN algorithm is modified to suit the membrane computation paradigm rather than classical computation. DEC P systems integrate the idea of a modified adaptive DBSCAN into original cell-like P systems with cell creation [29]. The robustness of the proposed method will be tested by detecting CNV lesion areas with distinct sizes, shapes and positions. Being a parallel computation model, the DEC P system can process all cases simultaneously, which reduces time consumption compared to the traditional saliency algorithm.

Data acquisition
OCTA data was acquired from patients with neovascular AMD, recruited at the Casey Eye Institute of Oregon Health & Science University. Participants were enrolled after informed consent in accordance with an Institutional Review Board/Ethics Committee approved protocol in compliance with the Declaration of Helsinki.
A single eye per patient was imaged. If one of the eyes had dry AMD, the eye with wet AMD was selected to be imaged. If both eyes had wet AMD, the one with sub-foveal CNV was selected. Imaging was performed by a commercial spectral domain OCT system (RTVue-XR; Optovue, Fremont, CA) with center wavelength 840 nm, full-width half maximum bandwidth of 45 nm, and an axial scan rate of 70 kHz [5]. Two 3 mm × 3 mm scans (one in the horizontal priority scanning direction and one in the vertical priority scanning direction) with 2 mm depth were sequentially obtained in one eye of each participant. In the fast transverse scanning direction, 304 axial scans were sampled to obtain a single B-scan and two repeated B-scans were acquired at every location. A total of 304 locations were sampled in the slow transverse direction, completing 608 B-scans within 2.9 seconds. Flow signal was detected by using the split-spectrum amplitude-decorrelation angiography (SSADA) algorithm [2,3]. Motion artifacts were removed by registration and merging of the horizontal priority and vertical priority scan volumes by motion correction technology (MCT TM ) [30]. Available scans with signal strength index (SSI) lower than 55 were considered of inadequate quality and excluded from analysis.

Algorithm overview
A framework of the automated CNV segmentation based on DEC P systems is shown in Fig.  1. Projection artifacts were removed by the reflectance-based projection-resolved OCTA algorithm [12] and retinal layer segmentation was performed by an algorithm based on directional graph search [31,32]. En face outer retinal angiograms were generated by maximum projection of the flow signal between the posterior boundary of the outer plexiform layer and the Bruch's membrane [31]. A pre-processing step consisting of a 5 × 5-pixel Gaussian filter was first applied to smooth angiograms ( Fig. 1, step 2). The distance of every non-zero pixel to all other non-zero pixels was recorded and the sorted distribution of distances was fit by a Gaussian model to extract the parameters ε and MinPts ( Fig. 1, step 3) used by the membrane computing algorithm described in section 2.3. A DEC P system was used to recognize the CNV vascular pattern (Fig. 1, step 4). By this method, several clusters were generated and the one with the maximum number of points was recognized as the cluster containin the pixels in the CNV vessels. Finally, morphological operations were applied to generate the CNV lesion area (Fig. 1, step 5). The software was implemented using Matlab 2017a (Mathworks, Natick, MA). Framework of CNV identification based on a DEC P system. A projection-resolved en face outer retinal angiogram was obtained in step 1 [31]. A Gaussian filter was used to smooth the CNV structure and the position of non-zero pixels are extracted in step 2.
Step 3 sorts the distances of non-zero pixels to every other non-zero pixel and fits them into a Gaussian function to obtain parameters ε and MinPts, as explained in section 2.3.
Step 4 describes the execution of DEC P system for CNV detection, explained in detail in section 2.3.
Step 5 shows the clustering result of the DEC P system method and morphological operations were used to generate the final CNV lesion area.

CNV segmentation based on DEC P systems
We propose DEC P systems as a parallel implementation of DBSCAN to discriminate the areas occupied by CNV from surrounding noise. A DEC P system is a type of cell-like P systems, a computational paradigm that relies on a nested cell structure in which cells contain multisets of objects that evolve according to given rules in a non-deterministically and optimally parallel manner [ Fig. 1, step 4]. The cells forming the DEC P system are enclosed by active membranes for the purpose of clustering optimization. In living cells, membranes are active, meaning that they allow cell generation and dissolution through mitosis (division), autopoiesis (creation) and cell lysis (dissolution). These properties allows the generation of new workspaces and the dynamic removal of useless workspaces. Since DBSCAN does not have a priori information about the number of clusters, its implementation by P systems necessarily implies using active membranes. Here, what we mean by membranes is an abstraction encasing the cells or fundamental units of the computational P-system, which assists the task of clustering CNV pixels by enforcing cell creation, dissolution and communication rules. They are not to be confused with the CNV lesion itself, which has been called a "CNV membrane" in the literature.
Initially, the DEC P system receives a map of positions of CNV-candidate pixels Fig. 1, step 2] of en face outer retinal angiograms after removing projection artifacts and application of a Gaussian filter. The density associated with every candidate pixel is obtained by counting the number of other candidates within a region of specified radius ε . Pixels with an associated density above a certain threshold (MinPts) are promoted to the status of core objects of a cluster. From all clusters found in this first step, new clusters will be created according to evolution and communication rules defined further below. Those rules will be applied depending on four pixel-to-pixel relationships: (1) whether the pixel under examination is a core object and whether its neighboring pixels are (2) directly density-reachable, (3) density-reachable or (4) density-connected to it. Pixel p X is directly density-reachable from q X if p X is within the ε -neighborhood of q X , and q X is a core object. p X is density-reachable from t X with respect to ε and MinPts if there is a chain of core objects  Fig. 2. Additionally, a cluster C is a non-empty subset of X following the criteria: (a) , X p q X ∀ : if q X C ∈ and p X is density-reachable from q X with respect to ε and MinPts, MinPts. C is initialized by core objects, each of which form a cluster by itself in the beginning. Fig. 2. p X and q X are directly density-reachable to each other; p X and o X are densityconnected to each other by q X ; p X is density-reachable from t X through where, the initial cell structure μ is the nested structure  values. We found empirically that due to the characteristics of proximity between CNV pixels, analysis of neighborhoods of radius (5) k ε > pixels did not provide improved performance. For this reason we disregarded values corresponding to bins 6 to 200 and only tested 5 pairs of ( ε , MinPts) values, which alleviates the difficulty of pre-setting parameters in the DEC P systems algorithm. The first five bins contains all the CNV pixels and some of the noise. Through comparing the value of f of the five bins and the distribution of CNV pixels in each case, we can obtain an approximate range of f . For each cell, k was selected according to function f shown in Table.1. The algorithm goes from cell 6 [] to 2 [] with ε increasing from (1) k to (5) k . Cell subjected by evolution rule

Value of ε and MinPts
Value range of f

Results and discussion
Scans of a single eye of 22 patients with neovascular AMD were processed serially by saliency method and parallelly by DEC P systems. The total time consumed in processing all 22 scans by DEC P systems was 37.9 s, which is about 23 times faster than the saliency implementation. CNV lesions with different characteristics in shape, size and position could be detected accurately. When applied on the old database used to compute an 83% accuracy of our previous saliency method [14], DEC P-systems exhibited an accuracy of 87%. In the new database used for the present work, the similarity with manual delineation was much higher for the DEC P systems compared to the saliency method (p<0.001, based on the paired Wilcoxon rank-sum test) [Fig. 3]. Positive and negative error rates were reduced by DEC P systems ( Table 2).
The saliency value at each position is calculated by a proportion of distance, orientation and brightness between non-zero neighboring pixels. This method does not expect the object of interest to occupy a significantly large section of the image, leading to a considerable inaccuracy in Fig. 3(case 3) where only the area within the central vessels of the CNV architecture is detected as visually "salient". Another source of error in the saliency method is associated to the the vascular pixels with lower flow signal than noise pixels [Fig. 3, case 1, case 2], in which case vascular areas could become less "salient". Although the saliency method is capable of preserving relatively low signal vascular pixels if they are proximal to each other, it might fail when that area of the vasculature architecture appears scattered. Thus, it usually cuts the full length of vessels at the boundaries of the CNV lesion [ Fig. 4], resulting in an underestimation of the lesion size that is calculated after morphological dilation. Alternatively, the DEC P systems can better segment the CNV vascular patterns [ Fig. 4] because it operates by clustering the neovascular pixels based on their spatial characteristics, as long as their flow signal is above a pre-defined threshold. Fig. 3. Performance of the DEC P systems and the saliency method compared to manual delineation for three examples. Each row case 1, case 2 and case 3 represent an outer retinal angiogram with CNV. Green and blue arrows point at positions where DEC P systems and saliency differ from manual segmentation respectively. The regions of disagreement are noted in the manual segmentation images with dashed ellipses. Note that the ellipses representing disagreement with the DEC P systems method (green color) are smaller in size than the ones representing disagreement with the saliency method (blue color). Jaccard similarity was larger for the DEC P systems segmentation. Fig. 4. Performance evaluation of the DEC P systems method and the saliency method for extracting CNV vessels in Fig. 3 (case 1). Blue arrows point at positions that saliency method cuts the full length of vessels at boundaries, whereas, DEC P system recognize the boundries well.
A critical step in the DEC P systems algorithm is the appropriate selection of the upper boundary of ε and MinPts values. As every clustering method, DEC P systems need to recognize the characteristics of the unknown object without supervision, and relies on the parameters predefined by the operator to bound the process of cluster population for successful implementation. This value of ε sets the largest of all possible vicinity radii over which density-connected points can be incorporated into clusters. It prevents noise pixels from entering the cells [] p that collect neovascular pixels, given that noise pixels are generally more distant from vascular pixels than adjacent vascular pixels from each other. We have empirically chosen a range of ε and MinPts from  Just as neovascular pixels typically exhibit a range of flow signal values greater than the flow signal of background, it is also possible that noise pixels in the vicinity of neovessels might occasionally exhibit similar decorrelation values to that of vascular pixels. This is owing to the linear dependency of bulk motion flow signal on the local OCT reflectance of tissue [33] and to the typical speckle characteristcs of interferometric imaging methods. The DEC P system method can't accurately discriminate vessel information from noise if noise pixels are bright and found within the vicinity of vessels, such as in the upper-left corner of the CNV area in [Fig. 6]. Since these types of pixels can't be removed by neither clustering based on pixel intensity nor by clustering based on pixel proximity, enhanced removal of projection artifacts is necessary to solve such inaccuracies. On the other hand, instead of flow signal of noise pixels being too high, flow signal in neovascular pixels could be too low (under threshold), also leading to inaccuracies owing to discontinuities. Since there is a thresholding step at the beginning of this method [ Fig. 1], neovascular pixels cannot be incorporated into the cluster even if they are proximal to other neovascular pixels just because they have been filtered out in a previous step. In order to minimize the occurrence of this error, a reflectance-adjusted flow signal thresholding method developed recently for OCTA should be considered in the pre-processing step [34]. Fig. 6. Performance of the DEC P systems method and the saliency method compared to manual delineation in the scan that noise pixels are bright and found within the vicinity of vessels. Green arrows point at positions where DEC P systems differ from manual segmentation. The regions of disagreement are noted in the manual segmentation images with dashed ellipses.

Conclusion
We developed a method based on a DEC P system to recognize CNV on en face outer retinal angiograms automatically. The method improved CNV area accuracy compared to a previous CNV saliency method and CNV lesion areas were segmented more precisely. The algorithm only needs one feature (position) to detect the group of CNV flow pixels, which together with its parallelism could reduce the processing time of a pool of scans. These results suggest a potential for improvement of conventional segmentation algorithms by using membrane computation models. For future works, it is of interest to explore P systems in improving biomedical image processing algorithms in terms of speed, storage, and convergence.