Deriving external forces via convolutional neural networks for biomedical image segmentation

: Active contours, or snakes, are widely applied on biomedical image segmentation. They are curves deﬁned within an image domain that can move to object boundaries under the inﬂuence of internal forces and external forces, in which the internal forces are generally computed from curves themselves and external forces from image data. Designing external forces properly is a key point with active contour algorithms since the external forces play a leading role in the evolution of active contours. One of most popular external forces for active contour models is gradient vector ﬂow (GVF). However, GVF is sensitive to noise and false edges, which limits its application area. To handle this problem, in this paper, we propose using GVF as reference to train a convolutional neural network to derive an external force. The derived external force is then integrated into the active contour models for curve evolution. Three clinical applications, segmentation of optic disk in fundus images, ﬂuid in retinal optical coherence tomography images and fetal head in ultrasound images, are employed to evaluate the proposed method. The results show that the proposed method is very promising since it achieves competitive performance for diﬀerent tasks compared to the state-of-the-art algorithms.


Introduction
Biomedical image segmentation plays an important role in developing computer aided systems for analysis of the images automatically.Numerous biomedical image segmentation methods [1][2][3] have been proposed in last few decades.According to the number of segmented objects in the images, these segmentation methods can be grouped into single object segmentation methods and multi objects segmentation methods.For single object segmentation, parametric active contour models are one of most widely used approaches, whose first version was proposed by Kass et al. in [4].There are several advantages with parametric active contour models, for example, they can offer smooth and closed contours as segmentation results and obtain high accuracy of object boundaries [5].However, for the classical active contour models [4], the active contours can not progress into boundary concavities.Besides, the initial curves are required to be close to object boundaries or else they will converge to wrong places.Many methods have been proposed to handle these limitations, in which gradient vector flow (GVF) [6] is one of most effective among them.Even though GVF is insensitive to initial curves and can drive the curves into boundary concavities, it is sensitive to noise and false edges.Fig. 1 gives an example, in which Fig. 1(a) is the GVF obtained by the method introduced in [6].It is observed that the GVF is disordered, which can not drive the initial curve into object boundary, as shown in Fig. 1(b), where the red curve is initial curve and the blue is final result.Abundant improved versions for GVF have been reported in the literature.For instance, Tang et al. [7] extended the classical GVF snake into multi direction GVF snake so that it can trace the boundary of skin cancer.Zhu et al. [8] used a minimal surface and component-normalized method to improve the GVF.These improved versions, generally, were careful designed based on some domain knowledge, e.g. the intensity distribution of the region of interest.Different from these domain knowledge based methods, in this paper, we propose using GVF as reference to train a convolutional neural network (CNN) [9] to derive an external force for the parametric active contour models.CNNs belong to end-to-end trainable methods.It allows us to avoid making hypotheses during designing the external forces, as a result, the dependence on the prior knowledge is abated, while the quality of the external forces is improved compared with those obtained by the classical methods [6,10] (Quantitative comparison is given in Section 4).Fig. 1(c) shows an example.It is observed that the derived external force is ordered and points to object boundary, which can ensure the active contour models to exclude the interference from noise and drive the initial curve into the right place, as shown in Fig. 1(d).

Related works
As the revival of CNNs in recent years, some researchers have tried to employ CNNs to facilitate the application of active contour models in biomedical image segmentation domain.For instance, Srivastava et al. [11] employed a deep neural network to perform a preliminary optic disk segmentation in fundus images.The preliminary segmentation was then refined further using active shape models.Avendi et al. [12] applied CNNs to automatically detect left ventricle (LV) chamber in MRI dataset.Stacked autoencoders were then used to infer LV shape.The inferred shape was finally incorporated into deformable models to improve the accuracy of segmentation.Li et al. [13] used CNNs to segment myocardium of LV first.Active contour models were then applied to segment endocardium based on the initial segmentation.
The above methods [11][12][13] employed CNNs to highlight the objects directly or indirectly so that active contour models can capture the objects easier.Besides, several different methods are also presented to facilitate the application of active contour models based on CNNs.Such as, Hoogi et al. [14] employed CNNs for adaptive estimation of active contour parameters.Le et al. [15] boosted the classic active contour methods to a new level of learnable deep learning approaches by reformulating active contour as deep recurrent neural networks.Different from these works, in this paper, we employ CNNs to derive an external force for the parametric active contour models.There are two advantages with the proposed method: 1) CNNs belong to end-to-end trainable methods, which can abate the dependence on the prior knowledge during designing the external forces.2) It can ensure the initial curves close to object boundaries based on the finding that the magnitude of the external force derived by CNNs around object boundaries is larger than other places.
To demonstrate the effectiveness of the proposed method, we apply the proposed method on different clinical applications, which are segmentation of optic disk in fundus images, fluid in retinal optical coherence tomography (OCT) images and fetal head in ultrasound images.These three applications play an important role in developing computer aided systems in their respective domains.Such as, optic disk segmentation is an important part in the computer aided systems for diagnosis of glaucoma.Fetal head segmentation is a key step to estimate the circumference, which is an important measurement to estimate the gestational age and monitor the growth of fetus.Many different methods have been designed for these three tasks.Due to the main motivation of this paper is to explore whether the parametric active contour models with the external force derived by CNNs is effective for single object segmentation in biomedical images, thus we do not offer a detailed introduction for these methods.For a comprehensive reading about these methods, we refer readers to papers [16][17][18][19][20][21][22] for optic disk segmentation, papers [23][24][25][26][27][28] for fluid segmentation and papers [29][30][31][32] for fetal head segmentation.
The rest of the paper is organized as follows: the details of the proposed algorithm are described in Section 3. Section 4 presents experimental results, including the comparison with ad hoc algorithms.Discussion and conclusion are made in Section 5.

Methods
We begin by first reviewing the parametric active contour models, then formulating an active contour model with the external force derived by CNNs.Finally, we show how to initial the active contours based on the external force magnitude.

Active contour models: a quick review
Active contours, or snakes [4] are curves defined within an image domain that can move to object boundaries under the influence of internal forces and external forces, in which internal forces are used to hold the curves together and to keep them from bending too much.External forces are used to drive the curves toward target boundaries [6].
The process of using active contour models to move a curve x(s) = [x(s), y(s)], s ∈ [0, 1] into an object boundary is formulated as a problem of energy functional minimization where the weighting parameter α controls the tension of x(s) and β controls the rigidity.x (s) is the first derivatives of x(s) with respect to s and x (s) is the second derivatives.E ext is the external energy term.
To minimize E, a curve x(s) must satisfy the following Euler equation Equation ( 2) can be regarded as a force balance equation where internal force term F int = αx (s) − βx (s), external force term It is noted that internal forces are computed from the curve x(s) itself.External forces are general derived from image data.One of popular external forces is GVF, which is computed as a diffusion of the gradient vectors of a gray-level or binary edge map derived from an image [6].Mathematically, GVF is a vector v(x, y) = [u(x, y), v(x, y)] that minimizes the following energy functional where f is an edge map derived from an image.The parameter µ is a regularization parameter controlling the tradeoff between the first term and the second term in equation ( 4).∇ is the gradient operator.u x , v x and u y , v y are partial derivatives of u and v to x and y respectively.GVF is achieved by solving the following Euler equations where ∇ 2 is the Laplacian operator.For the details of numerical implementations of equations ( 2) and ( 5), we refer readers to papers [4] and [6] respectively.

Deriving external forces via CNNs
As introduced above, GVF is sensitive to noise and false edges in an image.To handle this problem, in this work, we propose using GVF as reference to train a CNN to derive an external force that points to the target boundaries, so that the active contour models can exclude the interference from false edges and noise with the derived external force.Mathematically, an active contour model with the external force obtained from CNNs is formulated as where φ(I, w) is a network, termed GVF-Net for future reference, which is trained by using GVF as reference to derive an external force from a raw image I.
Given a training set of images and their referenced GVFs {I i , v i } N i=1 , the weights w of GVF-Net φ(I, w) are learned to minimize the mean squared error [9] It is noteworthy to point out that the referenced GVFs are computed from the boundary ground truth with equation ( 5), which can guarantee the referenced GVFs to point to the target boundaries.In this work, we employ u-net structure [33] for GVF-Net, which can ensure the output without reducing size after a series of operations.Note that GVF is a matrix with size of m × n × 2, where m × n is the size of an image.As shown in Fig. 3, the left part of GVF-Net is termed contracting path, which processes the input images through a series of convolutional layers.The size of each kernel in each convolutional layer is 3 × 3. Every two convolutional layers is followed by a 2 × 2 max pooling operation with stride 2 for downsampling.The number of feature map channels are doubled after each pooling step.The right part in Fig. 3 is termed expansive path, which upsamples the feature maps by a series of deconvolutional layers.The size of each kernel in each deconvolutional layer is 2 × 2. Each deconvolutional layer has a concatenation with the correspondingly cropped feature maps from contracting path and is followed by two convolutional layers.The final layer is a fully connected layer which used to map each component feature channel as the external force.In total, there are fourteen convolutional layers, three pooling layers, three deconvolutional layers and one fully connected layer in GVF-Net.

Initialization with external force magnitude
Curve initialization is one of key concerns with active contour algorithms.In general, an initial curve should be close to the truth boundaries or else it will likely converge to a wrong place [6].In this paper, the initial curves are obtained based on the external force derived by GVF-Net.We find that the magnitude of the derived external force around object boundaries is larger than other places.Fig. 4 gives an example, in which Fig. 4(a) and Fig. 4(b) are the external force and corresponding magnitude respectively.The dark regions in the magnitude image indicate that the force is very weak, while the gray regions indicate that the force is strong.If an initial curve is located in the strong force regions and surrounds the objects as much as possible, it is very likely to converge to the target boundaries.
Based on these findings, we propose an automatic initialization approach for fluid segmentation and another for optic disk segmentation respectively.As shown in Fig. 4, the first step of the automatic approach for fluid segmentation is to transfer the external force magnitude into a binary image with thresholding method [34].The maximum connected region in the binary image is extracted, as shown in Fig. 4(c).We then calculate the centroid of the maximum connected region.In addition, the major axis length and minor axis length of the ellipse that has the same normalized second central moments as the maximum connected region [35] are also computed.An original ellipse is obtained based on the centroid, major axis length and minor axis length.The original ellipse is then rotated with the interval of 10 degree around the centroid.The one containing the largest number of pixels belonging to the maximum connected region is selected as the initial curve, as shown in Fig. 4(d).For optic disk segmentation, we search a circle as the initial curve with circle hough transform [36] based on the binary images of the external force magnitude.
It is noteworthy to point out that the above automatic initialization approaches are not effective for fetal head segmentation.Hence, an interactive initialization method is proposed.As shown in Fig. 5(a), we draw a polygon as the initial curve based on the external force magnitude.The rules required to follow to draw a polygon are to ensure that each edge of polygon is located in the strong force regions and the polygon surrounds the object.We believe that the task of drawing such a polygon is simple since the strong force regions are obvious in the magnitude images.

Cascaded version
As introduced above, the external force is derived by CNNs directly.We believe that some existing schemes in machine learning community might have the capacity to improve the external force further.In this work, a cascaded scheme reported in [37] is employed to improve the external force.Concretely, the image-GVF pairs are first applied to train the GVF-net.The output of the trained GVF-net paired with the referenced GVF are then used to train another GVF-net.In this scheme, the second time training can be regarded as a fine tuning process, which is used to fine tune the external force obtained by the first GVF-Net.The effectiveness of the cascaded scheme will be verified by the experiments in Section 5.

For optic disk segmentation in fundus images
Three public databases for optic disk segmentation in fundus images are employed to evaluate the performance of the proposed algorithm, which are MESSIDOR database [38], DRIONS database [39] and ONHSD database [40].There are 1200 images in MESSIDOR database, 110 images in DRIONS database and 99 images in ONHSD database.Note that there are 9 images without discriminable optic nerve head in ONHSD database, which will not be employed to evaluate the algorithm.The annotation of optic disk in each image has been marked by experts.Since the resolutions of fundus images in different databases are different, for convenience, the method for optic disk localization in [41] is first used to locate the disk in a fundus image.A 301 × 301 sub image for optic disk segmentation is then determined based on the optic disk location.To train and test the GVF-Net, 60% of the images in MESSIDOR database are randomly selected as training set, 20% as validation set and 20% as test set.60%-20%-20% is a common used proportion to split a small dataset into training set, validation set and test set to train and test a CNN [42].All the images in DRIONS and ONHSD databases are used as test set.

For fetal head segmentation in ultrasound images
The data for fetal head segmentation in ultrasound images is from a grand challenge 2018 [43].The purpose of this challenge is to measure circumference of fetal head.The database contains a total of 1334 two-dimensional (2D) ultrasound images, which is divided into a training set of 999 images and a test set of 335 images.The size of each 2D ultrasound image is 800 × 540 with a pixel size ranging from 0.052 to 0.326 millimeters (mm).In this challenge, the annotations are available for training set only.According to the requirement of sponsors, the boundary of a fetal head should be fitted by an ellipse, which is described by five parameters, namely the center (x,y), semi major axis length, semi minor axis length and angle of rotation.The detected parameters by algorithms on test set are required to upload to sponsors to gain final results.In this work, we employ the direct least squares fitting of ellipses method [44] to fit the fetal head boundaries obtained by the active contour models.An example to show the result of fitting is given in Fig.

For fluid segmentation in OCT images
There are 481 OCT B-scan images for fluid segmentation.These images are collected from different OCT volumes obtained by ZEISS scanner in a local database.The resolution of each B-scan image is 512 × 1024.The fluid regions in these images have been depicted by an expert.To train and test the proposed method, these images are divided into training set, validation set and test set with the ratio of 60%-20%-20%.

Metrics
Different metrics are employed to evaluate the performance of the proposed method according to different tasks.

For optic disk segmentation and fluid segmentation
The performance of the proposed method for these two tasks is evaluated according to area overlap (AOL), Dice similarity coefficient (DSC), accuracy (Ac), true positive fraction (TPF) and false positive fraction (FPF) [17], which are defined as where TP is short for true positive, TN true negative, FP false positive and FN false negative.

For fetal head segmentation
The metrics used to evaluate the methods for fetal head segmentation are difference (DF), absolute difference (ADF), Hausdorff distance (HD) and Dice similarity coefficient (DSC) [29].DF, ADF and HD are defined as where HC R is the head circumference measured by the expert and HC S is the head circumference determined by algorithms.h(S, R) = max s ∈S max r ∈R ||s −r ||, in which R = {r 1 , r 2 , ..., r q } are the pixels in fetal head boundary from the expert, S = {s 1 , s 2 , ..., s p } are the pixels from algorithms.

Training details and parameter setting
A PC, equipped with an Intel (R) Core (TM) i7-4790 M CPU at 3.60 GHZ and 8 GB of RAM capacity, is employed to perform the experiments with MATLAB.We use the MatConNet [45] to train the CNNs.All the parameters of GVF-Net are initialized with Xavier initialization and trained for 66 epochs with the mini-batch size of 3 image instances.Training convergence can be observed within 60 epochs.For other hyperparameters, learning rate is set to 0.02, momentum 0.9 and weight decay 0.0005.Once training done, α and β are parameters which may influence the performance of the proposed method.However, we find that the performance of the proposed method is insensitive to these two parameters and we set α = 0.7 and β = 0 for the following experiments.

Performance analysis
Fig. 6 shows some examples for optic disk segmentation, fetal head segmentation and fluid segmentation respectively.It is observed that the proposed method works well for most of cases in the given examples, which demonstrates that the ability of generalization of the proposed method is high.However, the proposed method might fail for some cases, e.g.Fig. 6(c).Fig. 7 gives an example to explain why the proposed method might fail or succeed, where we can observe that for the case that the proposed method fails to detect the boundary, the external force derived by the GVF-Net dose not point to the target boundary.Conversely, if the derived external force points to the target boundary, it can guarantee the initial curve to be converged into the target boundary.The average values and standard deviations of different performance metrics for optic disk segmentation and fluid segmentation achieved on different databases are summarized in Table 1.One might wonder that why the performance of the method on MESSIDOR database is better than on DRIONS and ONHSD databases.The reason is that GVF-Net for optic disk segmentation is trained by the images in MESSIDOR database only.In the discussion section, we will demonstrate that the performance of the method on DRIONS and ONHSD databases can  be improved by using a few of images in these two databases to fine tune the trained GVF-Net.Table 2 lists the results of quantification for fetal head segmentation in ultrasound images.
A correlation analysis through a scatterplot is performed to assess the performance of the  proposed method in depth.Concretely, the correlation between the object area (the total number of object pixels) determined by the algorithm automatically and ground truth is given, as shown in Fig. 8.Note that the scatterplot for fetal head segmentation is not given here since the ground truth for test set is not available.The different styles of points in Fig. 8 mean the method performs on different types of images, which categorized according to different AOL intervals.It is observed that the points are very close to identity line for those images that the objects are segmented correctly (AOL>0.8).The percentage of images with AOL>0.8 in each database is summarized in Table 3, which illustrates that the proposed method is effective for most of images.

Comparison with hand designed external forces
To analyse the effectiveness of the proposed method, we compare the external force obtained by the GVF-Net with the hand designed external forces, including GVF [6].Table 4 summarizes the results.It is observed that the results obtained by the proposed method are improved significantly compared with those obtained by the active contour models with GVF or VFC as external force, which illustrates that employing CNNs to derive an external force for active contour models is a very effective way.Besides, we also employ two neural networks for comparison.One is the original u-net introduced in [33] and we use the Dice coefficient as the loss function for it.The other is the network introduced in Fig. 3 with the same loss function as GVF-Net.Note that these two neural networks are trained with image-region pairs.The results are summarized in Table 4.It is observed that the proposed method outperforms the u-net whether in terms of AOL or DSC, which verifies the effectiveness of the proposed method further.

Comparison with ad hoc algorithms
In this section, the ad hoc methods, namely designed for a specific task mentioned in this work, are selected for comparison.

For optic disk segmentation
Since most of the existing methods for optic disk segmentation are mainly evaluated according to AOL, thus AOL is selected as the metric for comparison in this section.If a method in the relevant literature has been not evaluated on an analysed database, the result will not be included in the table.
The existing methods for optic disk segmentation can be divided into unsupervised based methods and supervised based methods.Table 5 summarizes the comparison with several unsupervised methods.Even though the results obtained by the proposed method are superior to those obtained by the listed unsupervised methods, it might be unfair for these unsupervised methods since the test set used for evaluating the proposed method and the unsupervised methods is different.However, the images in DRIONS and ONHSD databases are not employed to train the GVF-Net, compared with the mathematical morphology based method [17], our method performs better on these two databases, which indicates that our method is very competitive.
Several supervised based methods are also selected for comparison, including a structured learning based method [22] and a deep learning based method [21].Even though the training set and test set used to evaluate these supervised methods are different, we only employ the images in MESSIDOR database to train the GVF-Net while obtain very competitive results in DRIONS and ONHSD databases compared with the structured learning based method [22].It is noteworthy to point out that K-fold cross validation is employed to construct training set and test set in [22], in which K=12 for MESSIDOR database, K=10 for DRIONS database and K=9 for ONHSD database.Since the code of the existing methods for fluid segmentation is not open-sourced.Thus we re-implement the graph cut based method [25] and the random forest based method [28] and apply them on the same database with the proposed method for comparison.Besides, we also compare with the enface fundus driven method [23] and blob detection based method [27] via citing the best results recorded in the literature directly.Note that DSC is a common used metric to evaluate the performance of these methods, thus we compare with these methods based on DSC metric.Table 6 summarizes the results, where we can observe that the result obtained by the proposed method is very competitive with those obtained by the listed methods, which illustrates that the proposed method is very promising.

For fetal head segmentation
We can make direct comparison among different methods for fetal head segmentation since these methods have been evaluated on the same test set with the same metrics.The results obtained by different methods can be found from the home page of the challenge [43], in which some of them obtained by the state-of-the-art methods are listed in Table 7.It is observed that the proposed method is very competitive with the methods for this task also, which verifies the effectiveness of the proposed method further.

Discussion & conclusion
In this paper, we propose using GVF as reference to train a CNN to derive an external force for the parametric active contour models.The derived external force is also used for initializing the active contours.As a result, the quality of external force influences the performance of the proposed method directly.We set an experiment to verify this conclusion.Concretely, 20% of the images in DRIONS database are selected randomly as training set and 20% as validation set to fine tune the GVF-Net trained on MESSIDOR database.The remaining 60% is used as test set.The same setup is also for ONHSD database.Table 8 summarizes the results.It is observed that the results obtained by GVF-Net with fine tuning are better than those without fine tuning, To improve the external force, we employ a cascaded strategy developed in machine learning community, which has been introduced in Section 3.4.The effectiveness of the cascaded strategy can be verified by the results summarized in Table 7.We believe that other existing schemes in machine learning community might also have the capacity to improve the external force further, e.g. by improving the structure of the network.Besides, integrating the shape prior knowledge about the objects into the learning process [46] would be also another potential way to improve the external force.
In summary, in this paper, we propose a novel algorithm for biomedical image segmentation by deriving an external force via CNNs for the active contour models.Abundant experiments show that the proposed algorithm is very promising since it obtains very competitive results for different tasks.

Fig. 1 .
Fig. 1.An example to demonstrate different external forces obtained by different ways.(a) The GVF derived from the gradient of original image.(b) The detected result (blue curve) with (a).(c) The external force derived by the trained CNN.(d) The detected result (blue curve) with (c).

Fig. 2 .
Fig. 2.An example to demonstrate the process of obtaining the referenced GVF.(a) Raw image.(b) Boundary ground truth of (a).(c) Referenced GVF.

Fig. 4 .
Fig. 4.An example to demonstrate the process of automatic initialization of curve.(a) The external force.(b) The magnitude of (a).(c) Maximum connected region of (b).(d) The initial curve (red) and final result (blue).

Fig. 5 .
Fig. 5.An example to show the interactive initialization.(a) The external force magnitude and initial curve.(b) The fetal heal boundary obtained by active contour model (red curve) and the result of ellipse fitting (blue curve).

Fig. 6 .Fig. 7 .
Fig. 6.Examples for optic disk segmentation(first row), fetal head segmentation (second row) and fluid segmentation (third row), where green curves are ground truth and blue curves are curves detected by the proposed algorithm.

Fig. 8 .
Fig. 8. Scatterplot of region area between manual segmentation and automatic segmentation.

Table 1 .
The average values and standard deviations of different performance metrics achieved on different databases.

Table 2 .
The performance of the proposed method for fetal head segmentation.

Table 3 .
The percentage of images with AOL > 0.8 in different databases.

Table 4 .
Average AOL and DSC obtained by different methods.

Table 5 .
Comparison with different methods for optic disk segmentation with AOL metric.

Table 6 .
Comparison with different methods for fluid segmentation with DSC metric.

Table 7 .
Comparison with different methods for fetal head circumference measurement.

Table 8 .
Comparison with and without fine tuning on DRIONS and ONHSD databases.the quality of the external force plays a leading role in the performance of the proposed method.