Attraction Propagation: A User-Friendly Interactive Approach for Polyp Segmentation in Colonoscopy Images

The article raised a user-friendly interactive approach-Attraction Propagation (AP) in segmentation of colorectal polyps. Compared with other interactive approaches, the AP relied on only one foreground seed to get different shapes of polyps, and it can be compatible with pre-processing stage of Computer-Aided Diagnosis (CAD) under the systematically procedure of Optical Colonoscopy (OC). The experimental design was based on challenging distinct datasets that totally includes 1691 OC images, and the results demonstrated that no matter in accuracy or calculating speed, the AP performed better than the state-of-the-art.


Introduction
The colorectal cancer now has been certified as the third most common cancer throughout the world. There will be close to 2 million new cases annually, and as the Top 4 disease mortality, in each year, around 600 thousands of people died from it [1,2]. According to statistics, most cases of colorectal cancer are pathological changed from colorectal polyp, therefore, to identify the types of colorectal polyp at early stage and to target it to solve the problem may reduce the death rate, and even cure patients [3,4]. Colonoscopy as the golden standard method has been widely used in screening and colorectal polyp diagnosis [5]. However, even a well-trained experienced medical staff still needs to spend tremendous time on suffering the optical colonoscopy (OC) images and colorectal polyp analysis. So to combine the Computer-Aided Diagnosis (CAD) system into the screening and colorectal polyp diagnosis could alleviate the burden of medical staffs' routing works.
Generally, a CAD system based on OC image consists of several procedures: image preprocessing, segmentation, feature selection, and classification. Among those procedures, image segmentation is a key procedure in the CAD system. Accurate segmentation of abnormal regions can improve the performance of the CAD system significantly [6]. Because current diagnoses based on endoscope have no calibration and are more operator-dependent than other diagnoses(such as X-ray or MRI), we may encounter a high inter-observer variation Although weighted graph method is valid for various types of image segmentation, it encounters small cut problems during image segmentation, especially when low numbers of foreground and background seeds are marked. The underlying cause of these problems lies in that the weighted graph or diffusion distance [37] is calculated only from adjacent information (such as Laplacian matrix) as solving the optimal problem of graph cuts, and the transition probabilities are difficult to propagate to the region far from the seeds. These existing problems greatly reduce the accuracy and quality of image segmentation. In recent years, the traditional weighted graph method has been extended to some new versions [37,38]. In particular, Couprie et al. [37] suggest a unified framework for GC, RW, SP and Watersheds. By combining the RW with the Maximum Spanning Forest, their algorithm produces an improved segmentation of a higher speed. However, in some cases, this algorithm is also affected by the small cut problem in the same manner as the RW.
In this study, we propose an interactive segmentation algorithm called attraction propagation (AP) for OC image segmentation, which compensates for some of the shortcomings of interactive methods such as AC, LS, GC, and RW. We introduce an attraction scheme based on a shape probability region and reconstruct an image graph for segmentation. In contrast to existing interactive methods, this new implementation scheme AP can solve the small cut problem in an effective manner by boosting the transition probability of pixels in attraction region. When encountering weak boundary problems, AP also produces more accurate segmentation results than many other methods such as AC, LS, RW and PW. Now we summarize the main contributions of this proposed method: • We provide a feasible strategy for the initialization of only one foreground seed. This strategy can be very easily integrated into the current clinical diagnosis and OC image acquisition.
• We introduce a shape probability region to offer an attraction region effectively, and propose a flexible framework for polyp segmentation in OC images.
• We give a fast minimization algorithm for the OC image processing.
• We build an open database of 800 OC images for assessing the performance of segmentation algorithms.
The remainder of this paper is organized as follows. Section 2 describes the proposed AP algorithm. Section 3 presents the experimental results obtained after the OC image segmentation. Section 4 discusses the sensitivity of parameters about our algorithm. In Section 5, we summarize the proposed algorithm.

Methods
In this section, we firstly give a user-friendly initialization of seed placement, and then introduce related works of weighted graph method. Secondly, we propose the attraction propagation algorithm with three forms. Finally, we provide some experimental examples to illustrate the effectiveness of our algorithms.

The initialization of seed placement
To segment OC images through interactive approaches, the primary problem lies in how to handle the initialization of seed placement. For the conventional initialization of seeds, users need to provide a large amount of interactive operations, which are laborious, exhausting and hard to be integrated into the current clinical diagnosis. In this paper, we provide a feasible and convenient strategy for initializing seeds. A schematic view of this strategy is given in

Related work
Given an image I, the segmentation problem can be formulated on an undirected and connected graph G = (V, E) [39], where vertices (nodes) v 2 V and edges e 2 E V × V, and the cardinalities of V and E are |V| and |E|, respectively. An edge that spans two vertices v i and v j is denoted by e ij . A weighted graph assigns a weight to each edge. The weight of an edge, e ij , is denoted by w(e ij ) or w ij , which is nonnegative (i.e., w ij ! 0). The edge set E comprises pairs of pixels, which are neighbors in the image. A Gaussian function of the L1 distance is defined between pixels with different image intensities, for which the weight w ij can be defined as where g i indicates the image intensity at pixel i and β is a free parameter. In weighted graph scenario, image segmentation is considered on an image domain (for the case with two labels). Every unlabeled node is given a label according to different energy function. The labels of nodes decide a segmentation of graph. Among these methods, the RW and PW are two powerful and robust tools [34,37], so we choose them in experiments for discussing and illustrating the effects of graph methods. Fig 3c and 3d show the segmentation results of RW and PW by choosing the seeds shown in Fig 3b. Under new strategy for initialization seed (see the previous subsection), these method may encounter the small cut problem.
To solve the small cut problem, we try to use the shape prior knowledge for enhancing the segmentation results. Actually, there are already several works discussing how to utilize the shape prior knowledge to increase the robustness and accuracy of medical image segmentation. For instance, Aslan et al. [40] use a 3D shape model obtained from a training data set to overcome any in homogeneity in CT images of bones. Chowdhury et al. [41] apply a probabilistic variation of the traditional graph cut algorithm and address the problem of segmentation of cerebral white matter from T1-weighted MRI data. To process cardiac MRI, Grosgeorge et al. [42] propose a segmentation method based on a statistical shape model obtained with a principal component analysis. Inspired by those works, we propose a new segmentation scheme called as AP which uses the linear system to give an optimization solution. By learning the shape priori-knowledge in OC databases, AP offers to determine attraction region and boost transition probabilities for nodes in the region.

The attraction propagation (AP)
In this subsection, we propose the attraction propagation (AP) to segment the input images. The attraction propagation term i ; t c i Þ is introduced in the weighted graph method to control the transition probability. The new segmentation model is as where the x c i is the transition probability of category c and c represents the category. In our problem c = f or c = b, where f and b represent the foreground and the background respectively. A is a set of vertices in the attraction region (AR). We will give a detailed description on how to obtain the attraction region in the next subsection. In our model, a i ; t c i and w ij are three adaptive parameters. The summand in the attraction propagation term can choose several forms according to various problem-specific domains. In this paper, we give three basic forms Pearson distance [43]: Inner product:  L2 norm square: Basing on the AP model (2), we can compute the transition probability x c i , and assign a label to every unlabeled node for getting the final segmentation. By expanding the attraction propagation and the weighted graph terms, we rewrite them in the matrix forms where ⋅ denotes an element-by-element multiplication and τ −1 is an element-by-element reciprocal of a column vector τ. I is a column vector where each element is one.

The attraction region and some details of AP
In this subsection, we first illustrate how to obtain the attraction region A from the shape prior. And then we give the detailed description of formulas (6)- (9). A flow chart of the attraction region is shown in Fig 5a and a detailed procedure is given as follows: Step 1. Select a subset in the OC ground truth image database for generating the attraction region.
Step 2. Calculate the centroids of ground truth images and calibrate these centroids to the centre position.
Step 3. Count the number Fr(m, n) of pixels in the overlapping regions of calibrated images. Normalize the number Fr(m, n) by Where m and n represent the coordinates, Ns is the number of ground truth images. We call SP(m, n) as the shape probability matrix (see Fig 4).
Step 4. Choose a threshold matrix by where Γ is a threshold and 0 Γ 1. We use the non-zero region in the threshold matrix as the shape prior region (see Fig 4). An illustrative example of generating the shape probability matrix is shown in Fig 5b for Ns = 3. (6)- (8). Given one foreground seed, we introduce an attraction region A generated by shifting the shape prior region onto the seed such that the position of seed is same as the centroid of the region (see Fig 4).

Now we illustrate concrete definitions of symbols in formulas
The matrix A is a diagonal |V| × |V| attraction propagation incidence matrix and defined as follows The attraction factor a i is used to adjust the attraction intensity of pixels in the A. In this work, the attraction factor a i = exp(−|g i − g s |), where g s is the image intensity of seed and s is a foreground seed. In formulas (6)-(8), the range base τ in R |V| for controlling the scope of probabilitiesis defined by where T is a matrix with the same size of input image. Fig 6a gives an intuitive example to demonstrate matrix T. It is constructed by shifting the threshold matrix T in formula (11) onto the foreground seed such that the center coordinate of T is same as the position of seed. Vec(Á) operator stacks the entries of a 2-dimensional matrix into a column vector. We hope to boost the probabilities of pixels in the attraction region when computing the probabilities in the foreground, and depress those probabilities in the attraction region when considering the transition probabilities in the background as shown in Fig 6b. In formula (9), the |E| × |V| edge-node incidence matrix K is defined as and the |E| × |E| matrix W is a constitutive diagonal matrix with the square weights of each edge along the diagonal. Furthermore, the Laplacian matrix L of the graph is defined as follows Now we give the detailed description of formulas (6)-(9). Pearson distance: We denote A ¼ 1 2 At À1 I T Á I and L þ A as P, which is partitioned into marked (seed nodes) and unmarked (unseeded nodes) blocks.
Inner product: L2 norm square: We denote L + A as S, which is partitioned into marked (seed nodes) and unmarked (unseeded nodes) blocks.
For eachcategory c, we define the set of labels for the seed points as a function Q(v j ) = c for any v j 2 V M . Furthermore, we define a |V M | × 1 vector (where |Á| denotes cardinality) for each category c at node v j 2 V M as The minimization of (16), (18), (19) with respect to x c U is given by the linear system: Pearson distance: Inner product: L2 norm square: for one label, or for all labels, where the matrix X has K (in our problem, K = 2) columns for each x c U , where M has columns for each m c and T has columns for each t c U Therefore, we solve the probabilities matrix X using K sparse linear systems. The overall segmentation scheme is described in Table 1.
The above attraction scheme AP has many advantages and properties, such as high computational speed and robust segmentation. AP is formulated on a general graph with prioriknowledge shape, which can represent any dimension or topology, and it holds the segmentation accuracy of the abnormal region and reduces the estimates of the background region. An example of the results obtained by the AP is shown in the last three rows of Fig 7. Even when only one seed point is selected in the foreground region, the AP can improve the segmentation accuracy more effectively than other weighted graph methods. More numerical results are shown in the following Section.

Experiments and Results
To the best of our knowledge, there are no literatures offering a large-scale analysis and validation the performance of the automatic and interactive segmentation algorithms for OC image. It is still largely unknown about whether interactive/semi-automatic methods perform better Step1: Compute edge weights according to exp(−β|g i − g j |).
Step2: Compute the Laplacian matrix with formula (15), the attraction propagation incidence matrix with formula (12) and the range base with formula (13).
Step3: Decompose the Re-built Laplacian matrix and compute the transition probability by solving the sparse linear system: Step4: For any node v i , classify it as belonging to segment k (k represents the classify label) if x k U > x k 0 U for all k 0 6 ¼ k.  than automatic methods or vice versa. And then, almost all of the interactive method can realize automation, where the polyp was present in the center of the image (or fixed position in the image). So in this section, we present a comprehensive performance evaluation using some well-known segmentation algorithms (including automatic and interactive segmentation algorithms) for 1691 OC images with abnormal regions. In total, these images are derived from the three databases. Further details about the three different datasets are given below. In this study, the foreground seed is single and the backgroud seeds form a frame as shown in Section 2.1. And for each OC image, the same seed placement is employed by all the interactive segmentation methods. The performance of these algorithms is measured by annotated area covered (AAC) and Dice similarity coefficient (DICE) [44].
where S a is the resulting from algorithm annotation and A a is an image section resulting from manual annotation. O a is the number of common pixels between S a and A a . In our simulation, we split the data into two parts. One is for training the attraction region (AR) and the other is for testing the performance of our algorithm.

Colonoscopy image segmentation on database I
The database I is a comprehensive data set applied in the field of colonoscopy image detection and we use it to test the proposed approaches. This database comprises 358 OC images, each of which contains 500 × 574 pixels. Fifteen related approaches are used for comparison: Active Contours for Colonoscopy (ACC) [20], Chan-Vese-Segmentation (CV) [18], Distance Regularized Level Set Evolution (DRLSE) [28], Entropy Rate Clustering (ERC) [11], Fuzzy C Means (FCM) [7], Mean Shift (MS) [8], Normalized Cuts (NC) [9], Quick Shift (QS) [10], Power watershed (PW) [37], Random Walks (RW) [34], Sector Accumulation-Depth Of Valleys Accumulation (SA-DOVA) [12], Template Matching (TM) [17] and our method. Fig 8 shows an example of the output from each method using the database I. Our method (Fig 8a) identifies the polyp edges better than other methods shown in Fig 8. The underlying reason of this result is that the attraction strategy can take advantage of the structural properties of polyp (such as geometry symmetry and local swell) to improve the precision. The output of ACC (Fig 8b) appears to be too conservative because it has high precision but it fails to identify a larger number of polyp pixels. The output of CV (Fig 8c) does not retain the shape of the original polyp. The output of DRLSE (Fig 8d) is a suitable fit for OC image segmentation but the method is too slow to be used for real-time segmentation. The output of ERC (Fig 8e) is an effective method for OC images. However, a fixed parameter is not available that works well for all OC images. The output of NC (Fig 8h) is not a bad fit in our example, but it consumes a large volume of memory. Specially, the result of PW (Fig 8i) is very well in this example, but PW is unstable and may fail in other cases, which causes a weak performance result as shown in Table 2. SA-DOVA (Fig 8l) obtains a high AAC but this method could also decrease the precision. Weak boundary of polyp is a major factor that leads to the invalidation of SA-DOVA. The other algorithms (FCM, MS, QS and RW) have various deficiencies in this case and we do not discuss them here. Table 2 shows the performance of twelve algorithms on the database I according to two different measures: AAC and DICE, where the bold numbers show the top three results for each measure. In this experiment, we choose fixed parameters for our method (β = 310, Γ = 0.8; β = 340, Γ = 0.6). For algorithm ACC, CV, DRLSE, ERC, FCM, MS, NC, QS and TM, it is difficult to choose a group of fixed optimal parameters for all OC images. So we select optimal parameters for each image in order to improve the performance of these algorithms. In Table 2, DRLSE and PW get a low AAC and DICE. This result may be caused by the low quality of OC images in the database. For sparse foreground seeds, RW encounters small cut problems which cause the low values of its AAC and DICE. Algorithms FCM, MS, SA-DOVA and TM are affected by weak boundary and types of abnormal regions. Algorithms ERC, NC and AP have relatively high AAC and DICE.

Colonoscopy image segmentation on database II
The database II comprises 533 OC images of resolution 288× 384. Fourteen approaches are used for comparison and related results are shown in the Table 3. Compared to database I, most Attraction Propagation for Polyp Segmentation in Colonoscopy Images algorithms perform better on database II for the AAC and DICE since database I contains many low quality images. Most of the algorithms keep a similar rank ordering on performance as Table 2. As shown in Table 3, both AAC and DICE for PEA-AP are in the top three. Specially, the AAC of AP is 47% better than RW, which illustrates the importance of adding the attraction propagation term since the larger AAC value means more retrieval polyp pixels. As shown in Tables 2 and 3, the DICE of general segmentation methods (such as MS, ERC and NC) is close to the SA-DOVA, ACC or TM, which are exclusively used in OC images. And the interactive methods (such as DRLSE, PW and RW) are no better than automatic methods (ERC and NC). So we intend to provide a further measure and analysis in the following Section.

Colonoscopy image segmentation on database III
To promote the development of the colonoscopy image recognition field, we produce a rich set of images for researchers to analyze, i.e., an open database called the Northeast Normal University Colonoscopy (NNUC) dataset, which contains 800 OC images that cover more abnormal region types than databases I and II. The images in NNUC are derived from colonoscopy images of 400 patients, where experts print the screen after identifying potentially cancerous areas during operations. Fig 9 shows the different types of polyp appearances found in the dataset. The size of the OC images varies in different cases due to the diversity of the recording instruments employed. Thus, the areas of the images are cropped in order to reduce the effects of irrelevant information on the image segmentation process. After the cropping process, the images are resized to 408 × 474. Fig 9 shows some of the preprocess results, where the first row comprises raw images, the second row shows the preprocessing results, and the third row shows the polyp masks, which are regarded as the ground truth in the evaluation. Segmentation by retrieval. To evaluate the performance of different segmentation algorithms for OC image in a more comprehensive and systematic manner, segmentation is regarded as the retrieval of polyp or non-polyp pixels. All of the segmentation methods return the designated positions of the polyps in the images as binary masks. Thus, we perform a pixelby-pixel comparison of the polyp masks and the resulting masks to evaluate the segmentation methods. The segmentation methods assign each pixel as a polyp or a non-polyp pixel. Based on these comparisons, we obtain four possible outcomes: true positive (TP), true negative (TN), false positive (FP), or false negative (FN). We denote positive as a polyp pixel and negative as a non-polyp pixel. Thus, rather than using the two standard metrics in the evaluation, we employ six characteristic values to analyze the impact of the parameters on the segmentation results obtained with images from the NNUC database. The six characteristic values can be computed as: Recall ¼ TP=ðTP þ FNÞ; ð31Þ Specificity ¼ TN=ðTN þ FPÞ; ð33Þ Comparison of the results obtained with different segmentation algorithms. In this section, we compare the quality of the segmentation results obtained by using the methods tested in our evaluation. Thirteen related approaches are selected for this evaluation: ACC, DRLSE, ERC, FCM, MS, NC, PW, QS, RW, TM and AP. We apply these methods to all of the 800 images in the database III. Fig 10 shows some results of the output from each method using image from the database III. The optimal segmentation results are presented in Table 4, where the bold numbers show the top three scores for each characteristic. As shown in Table 4, among all the algorithms tested, our AP method achieves the highest Accuracy, Recall, F2measure, and DICE. It has seven characteristics in the top three. Although the Precision with AP is less than that with some other algorithms, the Recall (according to [44], the value of AAC is the same as that of Recall) is 22% better compared with PW. The DICE of AP is 28% better compared with that of RW. This demonstrates that our method can find more polyp pixels with small errors, and has greater flexibility for adapting to variation in the shape of objects.

Comparison of the processing time
We compare the time costs of various algorithms in Fig 11. The experiments are performed by using Matlab R2009b on an Intel Core(TM)2 Duo CPU with 2.93 GHz and 2 GB memory except PW (PW uses the available C++ libraries for operations on Linux-OS). As shown in Fig  11, MS is the fastest method, which required 0.02-0.03s, but its performance is fair. The solution times are slow with ACC and QS, i.e., 4.8-12.9s and 5-7s. The processing time with AP is close to that with RW, it costs 0.79-2.6s to calculate the solution for each image. This indicates that our method can be used for online processing.

Discussion
In the proposed method, there are two parameters, a weighted parameter β and a probability threshold Γ, to control the performance of segmentation algorithms. We discuss the sensitivity on the two parameters for the proposed method, and use the same initialization strategy for processing each OC images. The five characteristic values can be used to compare the robustness of our method. It is because that the Recall is the ratio of correctly detected polyp pixels to all polyp pixels. High values indicate a high number of pixels available for classification of the polyp. The Precision is the number of correctly detected polyp pixels divided by all polyp pixels. A high value indicates a low number of polyp pixels incorrectly regarded as nonpolyp pixels. The accuracy is the fraction of all pixels which were classified correctly. The specificity is the number of correctly detected non-polyp pixels divided by all non-polyp pixels. The F2measure is complementary metric on segmentation performance and it combines Precision and Recall into an only measure. First, we consider the performance of our algorithms based on various weighted parameter β. Fig 12 shows quantitative evaluation of the segmentation results with Γ = 0.7 in 5-fold validation. It can be seen that the Accuracy and the Specificity obtained by PEA-AP and L2-AP are not sensitive as β increases. Although the AAC is 0.70-0.82 and the F2measure is 0.73-0.79, the performance of the proposed scheme is acceptable. Also we note that INN-AP is sensitive as β increases in the Fig 12. This is because that compared to two other algorithms, the INN-AP only considers the shift operation for boosting or inhibiting the probability. Next, we vary the parameters β and Γ at the same time. The experiment results on the Precision, Recall, Accuracy and F2measure are shown in Fig 13. For PEA-AP and L2-AP, the Accuracy is large than 95%. The average AAC and F2measure are roughly from 0.70 to 0.80. This means that these two algorithms are robust for parameter selection.

Conclusion
In this study, we propose a flexible framework for interactive image segmentation in OC image segmentation. Compared with other segmentation methods, our new framework AP provides an efficient, real-time, robust solution. We introduce the shape probability region and the attractive strategy, which allow our AP algorithm to start with only one foreground seed, and automatically detect the contours of polyps. After adjusting the threshold parameter, the AP algorithm can restrain the transition probabilities of pixels in non-polyp regions without any pre-processing of the raw OC images. Moreover, it can effectively improve the accuracy of the segmentation of polyp regions, although they are very noisy with uneven brightness.
In order to evaluate the performance of our segmentation methods, we built an open database of 800 OC images that contained multiple types of abnormal regions. Furthermore, we use seven characteristics to compare the quality of various algorithms, and we demonstrated Attraction Propagation for Polyp Segmentation in Colonoscopy Images the processing efficiency with various database. Our experimental results show that AP algorithm is effective and practical for polyp segmentation in colonoscopy images.