Active contours driven by difference of Gaussians

In this paper, a novel edge-based active contour method is proposed based on the difference of Gaussians (DoG) to segment intensity inhomogeneous images. DoG is known as a feature enhancement tool, which can enhance the edges of an image. However, in the proposed energy functional it is used as an edge-indicator parameter, which acts like a balloon force during the level-set curve evolution process. In the proposed formulation, the internal energy term penalizes the deviation of the level-set function from a signed distance function and external energy term evolves the contour towards the boundaries of the objects. There are three main advantages of the proposed method. First, image difference computed using the DoG function provides the global structure of an image, which helps to segment the image globally that the traditional edge-based methods are unable to do. Second, it has a low time complexity compared to the state-of-the-art active contours developed in the context of intensity inhomogeneity. Third, it is not sensitive to the initial position of contour. Experimental results using both synthetic and real brain magnetic resonance (MR) images show that the proposed method yields better segmentation results compared to the state-of-the-art.

which global region-based methods are reformulated by replacing the global means with the image local information 24 . These methods can segment intensity inhomogeneous regions, unlike their global counterparts. However, these methods are sensitive to the position of the initial contour. Moreover, they also have high computational cost due to the complicated local information in their formulation.
A local active contour model for segmenting images with intensity inhomogeneity was proposed by Zhang et al. in 22 . Local image information is used to define a local image fitting (LIF) energy functional, which can be interpreted as a constraint on the differences between the fitting image 20,21 and the original image. Furthermore, a new method is used to regularize the level-set function by using Gaussian kernel filtering after each iteration. In addition, re-initialization of the level-set curve is not required. However, this method is sensitive to the initial position of contour and the final contour doesn't strictly follow the object boundary.
Alternatively, a region-based active contour method is formulated in the context of intensity inhomogeneity by utilizing local intensity means 23 . It uses a signed pressure force (SPF) function based on a local fitted image in its energy formulation in order to segment images with intensity inhomogeneity. A Gaussian kernel is used to smooth the level-set function after every step. Therefore, this method does not require re-initialization. This method has high time complexity compared to the traditional global region-based methods 12,13 .
A variational level-set approach for bias correction and segmentation (VLSBCS) for images corrupted with intensity inhomogeneity was proposed by Li et al. in 17,25 . The computed bias field is intrinsically ensured to be smooth by the data term in the variational formulation, without any additional effect to maintain the smoothness of the bias field. A local statistical active contour model (LSACM) for image segmentation in the presence of intensity inhomogeneity was proposed by Zhang et al. in 26,27 . In this work, the inhomogeneous objects are modelled as Gaussian distributions of different means and variances. A statistical energy functional is then defined for each local region, which combines the bias field, the level-set function, and the constant approximating the true signal of the corresponding object. Both of these methods are able to segment and correct intensity inhomogeneous images. However, these are also sensitive to the initial position of contours and have high time complexity.
A hybrid region-based active contours driven by local and global fitted image models (LGFIM) was proposed in the context of intensity inhomogeneity 28 . Local and global intensity information were used to both correct and segment the inhomogeneous regions. In this paper, two SPF functions i.e., local and global were also devised to stabilize the gradient descent solution. This method is not sensitive to the initial position of the contour; therefore, it provides similar segmentation results irrespective to the initial contour position. Moreover, it is also able to properly segment intensity inhomogeneous images. However, it has a high time complexity. Figure 1 shows intensity inhomogeneous image segmentation using active contour method. A traditional global region-based active contour method such as, CV 12 model is unable to segment intensity inhomogeneous object as shown in Fig. 1(a). In turn, Fig. 1(b) shows that the local region-based active contour method is able to properly segment intensity inhomogeneous object. However, local active contour methods have high time complexity because of the intensity mean computations in the local neighbourhood during the contour evolution process. Edge-based active contour methods can segment intensity inhomogeneous objects; however, these methods do not guarantee the level-set curve will stop at the object boundary, if it is blur or when the intensity difference between the background and object boundary is not clear.
In this paper, an edge-based active contour method driven by the difference of Gaussians (DoG) function is proposed in the context of intensity inhomogeneous image segmentation. The Gaussian image difference computed by the DoG function provides edge information of the global structure of the given image. This edge information is used as a balloon force in the proposed energy functional to evolve the level-set curve throughout the image structure. An energy penalizing term from 15 is used to regularize the curve and to maintain the level-set as a signed distance function, which also removes the need for computationally expensive re-initialization of level set. This work has three main contributions. First, it is able to segment the global structure of an image unlike traditional edge-based active contours 15,16 . Second, it is able to properly segment intensity inhomogeneous images with a low time complexity compared to local region-based methods [21][22][23][24] . Third, it is not sensitive to the initial position of the level-set curve.
The selection of the standard deviation parameter of the DoG smoothing kernels is critical, specially when the image is highly affected by noise (the parameter difference should be high). The proposed method is applied to both synthetic and real images to show the segmentation accuracy and robustness of the method.
Proposed method. Active contours are dynamic curves that evolve toward the object boundaries to partition an image into non over-lapping regions. In traditional active contours, the curve C is represented by the zero level-set, such that φ = | = C x y x y ( , ) ( , ) 0 of a level-set function φ (x, y). The evolution of the level-set function φ can be written in the following general form: which is referred to as the level set equation 9 . Function F is called the force function. For image segmentation, function F depends on the image data and the level set function φ.
Let I : Ω → R 2 be a given image, C a curve at which the level-set function φ(x, y) is zero x y x y ( , ) ( , ) 0. The energy functional E is defined as: In traditional edge-based active contour methods 14,29,30 , it is necessary to re-initialize (reshape) the level-set as a signed distance function during the curve evolution to properly follow and capture the object boundaries. Therefore, it is necessary to keep the evolving level-set function as an approximate signed distance function during the evolution, especially in a neighbourhood around the zero level-set. It is well-known that a signed distance function must satisfy the desirable property that φ ∇ = 1. An energy term P(φ) is proposed in 15 as a metric to characterize a function φ to a signed distance function in Ω ∈ R 2 , which helps to penalize the deviation of φ from a signed distance function during its evolution. The internal energy E int (φ) is defined as: where α is the scaling parameter of E int , which penalizes the energy leakage. In (2), E ext is the external energy of a function φ, which is defined as follows: where μ > 0 and v are constants, and the terms L(φ) and A Γ (φ) are defined as: , 1 2 where H ε (φ) is the regularized version of the Heaviside function: Parameter ε controls the smoothness of the Heaviside function. For ε → 0, the Heaviside function is the ideal unit step function. In (6), Γ σ σ , 1 2 is the difference of Gaussian function, which is used to replace a traditional edge indicator function. In (5), δ ε (φ) is the smooth version of the Dirac function, defined as: where parameter ε controls the width of the Dirac function. For ε → 0, the Dirac function is the ideal unit impulse. In this paper, a level-set method based on the difference of Gaussians (DoG) is proposed. A DoG function, which is equivalent to the Mexican Hat function, is a feature enhancement tool that involves the subtraction of a blurred version of an original image from another, less blurred version of the original. As a feature enhancement algorithm, the difference of Gaussian functions can be utilized to increase the visibility of edges and other details present in a digital image. The difference of Gaussians algorithm removes high frequency detail that often includes random noise, thus rendering this approach one of the most suitable for processing images with a high degree of noise. A major drawback of the application of the algorithm is an inherent reduction of the overall image contrast that results. In this work, it is employed as an edge-detector, which works as a balloon force in the external term of the proposed energy formulation during the level-set curve evolution. Let I : Ω → R 2 be an input image. The where σ 1 and σ 2 are the standard deviations of the first and second Gaussian kernels, respectively, where σ 1 < σ 2 . Figure 2 shows a 1D difference of Gaussians (DoG) function. It shows that the DoG function is zero when the slopes of both Gaussian functions intersect with each other. It helps to extract edges even when the images contain intensity inhomogeneous objects.
Finally, from (2), the proposed energy functional, which uses the difference of Gaussians (DoG) function Γ σ σ x ( ) , 1 2 to extract edge information, is defined as: In the above equation, the first term detects the edges using the DoG function. The second term regularizes the region and the third term penalizes the energy leakage. The first variation of the Gateaux derivative 31 of the functional E is denoted by φ ∂ ∂ E , which has the following relationship with the evolution equation: The above equation is the gradient descent flow that minimizes the energy functional E. For a particular functional E(φ) explicitly defined in terms of φ, the Gateaux derivative can be computed and expressed in terms of function φ and its derivatives 31 .
By calculus of variations 31 , the Gateaux derivative (first variation) of the functional E in (10) can be written as: where ∇ is the Laplacian operator. Therefore, the function φ that minimizes this functional satisfies the . The steepest descent process for minimization of the functional E yields the following gradient flow: In this paper, the spatial partial derivatives φ ∂ ∂x and φ ∂ ∂y are approximated by the central difference. The approximation of (13) using a central difference scheme can be written as: , is the approximation of the right hand side in (13) by the above difference scheme. The difference equation in (14) can be expressed as the following iteration: where τ is the time step used in the above numerical implementations. There is a close relation between the time step and the scaling parameter of the energy penalization term i.e., τ × α ≤ 0.2.
In level-set methods, it is essential to initialize the level-set function φ as a signed distance function (SDF) φ 0 . In the proposed formulation, not only is the re-initialization procedure completely eliminated, but the level-set function φ no longer needs to be initialized as an SDF. The initial level-set function φ 0 is defined as: In (16), ρ > 0 is a constant (ρ = 1 in this work) and t = 0 define the initial condition of the level-set function φ 0 = φ(x,t = 0). Ω 0 is the inner region of the initial level set φ 0 , Ω is the image domain and ∂Ω 0 the boundary of level set φ. Figure 3 shows a 1D profile from the middle column of the difference of Gaussian regularized image shown in Fig. 3(d). Figure 3(b) and (c) are Gaussian regularized images, which are subtracted to produce the DoG based edge profile of the given image, as shown in Fig. 3(d). Figure 3(e) shows a middle profile comparison between Fig. 3(b)-(d). It shows that when there is an intersection between the 1D middle profile of Fig. 3(b) and (c), then the edge is detected (a transition from high to low or vice versa), which is shown with the red line in Fig. 3(e).

Results
In this section, segmentation results using both synthetic and real images are discussed. The proposed method is implemented using MATLAB and run on a 3.4 GHz Intel Core-i7 with 16 GB of RAM, testing it on both synthetic images and real brain magnetic resonance (MR) images of 250 × 250 pixels with 256 grey levels (8bpp). The parameters used in all experiments in this section are: μ = 0.001 × 255 2 , v = −40, σ 1 = 1, σ 2 = 2, α = 1.0, ε = 1 and the time step τ = 0.1. Figure 4 shows segmentation results of different state-of-the-art methods using a flower image with different contrast variations until the flower object becomes inhomogeneous with respect to the image background. It shows that the DRLS method is able to properly segment images in the first four rows, but it is unable to segment the image in the last column. In turn, the CV method is able to properly segment the images in the first two rows and fails to properly segment the remaining images. Although the RSF method is able to segment objects in all images, the petals of the flowers in the images of the last two rows are not properly segmented. In turn, the qualitative comparison shows that the proposed method yields the best segmentation result for all the objects. Table 1 shows a CPU time comparison among the evaluated methods tested in Fig. 4. For the images in the first and second rows, the CV method is the fastest among all segmentation algorithms. In turn, the proposed method is the fastest for the images in the last three rows. Figure 5 shows that the DLSR method could segment the first three images properly. The CV method could not properly segment the objects in all of the images. The RSF method could segment the first two images and the last properly. However, in the last image, the final contour has missed some details of the small objects enclosed in the main object. Moreover, some details of the big object in the right corner are also missed. In turn, the proposed method yielded the best segmentation results for all of the images.  Table 2 shows a CPU time comparison among the evaluated methods tested in Fig. 5. It shows that the CV method is the fastest among all segmentation algorithms for the images in the first and the last rows. For the second and fourth images, the proposed method is the fastest. In turn, for the image in the third row, the DRLS method is the fastest. Figure 6 shows that only the proposed method and LGFIM are able to accurately segment all of the intensity inhomogeneous objects. LIF and LSACM are also able to segment all the objects. However, for both methods, the level-set curve around the boundaries of the objects is not quite smooth, which results in information loss. In turn, the VLSBCS method is able to properly segment the images in the first two rows and fails in the last three rows. Table 3 shows a CPU time comparison between the evaluated local active contour methods shown in Fig. 6. It shows for all images the proposed method yields the lowest CPU time; therefore, it is the fastest.
In Fig. 7, two images with and without noise are used to show the segmentation capability of the proposed method in the presence of noise. Images shown in the first and third rows are affected by the Poisson noise. Whereas, image in the second row is affected by the Gaussian noise with mean = 0.01 and variance = 0.2. The propose method is able to properly segment images with Poisson noise by using the default values of σ 1 and σ 2 as discussed earlier in this section. However, in case of Gaussian noise different values of σ 1 and σ 2 are used i.e., In order to segment images with intense noise, bigger values of σ 1 and σ 2 are used. Moreover, regularization parameter μ is also increased (μ = 0.01 × 255 × 255). It concludes that the proposed method is not affected by the presence of noise. It is able to properly segment the object both with and without background noise. Figure 8 shows different images with weak boundaries, salt and pepper and Gaussian noise. First row shows circular object with weak boundaries in both white and black background. It shows that in both cases, the proposed method is able to properly segment objects with weak boundaries. Second row shows image with five different objects with and without salt and pepper noise of d = 0.3. Where d is the noise density. Proposed method is able to properly segment all five objects with weak boundaries, when there is no noise in the image. In second case, when salt and pepper noise of d = 0.3 is applied; the proposed method could only segment four objects out of five. The fifth object was dissipated by the noise. In last row, both salt and pepper and Gaussian noise are applied to an irregular object with weak boundaries. In the first image, a salt and pepper noise of d = 0.2 is applied. In the second image, Gaussian noise of zero mean and variance = 0.3 is applied. The last row shows that the proposed method yields acceptable segmentation results for both cases of noise. Figure 9 shows the impact of the position of the initial contour on the evaluated methods. The segmentation results produced by the proposed method are not affected by the initial position of the contour, unlike the other evaluated methods, which are sensitive to the initial contour. Figure 10 shows the brain MR image segmentation problem using the evaluated methods. RSF and the proposed method are able to segment the detailed anatomical structure of the brain region, whereas the remaining methods fail to do so. There is region overlap in the segmentation result of RSF method shown in the second row, which concludes the proposed method yields the best segmentation.  Table 2. CPU time comparison among the state-of-the-art methods shown in Fig. 5.
which is the region of interest missed by the proposed method during the segmentation process.
From the above subsets, different similarity metrics are computed. In particular, the Jaccard index (JI) 32 , the Dice coefficient (DSC) 33 and the Matthews correlation coefficient (MCC) 34 are frequently used in set comparison that is, to compute the segmentation accuracy when the ground truth of the region of interest is available. In this paper, the three set similarity metrics (JI, DSC and MCC) are computed for the quantitative analysis. They are defined as:    For the maximum segmentation accuracy, the values of JI, DSC and MCC should be close to 1 (ideally 1). The Hausdorff distance (HD) 35 is another similarity metric, which is used to compute the accuracy between two  sets. It provides a symmetric distance measure of the maximal discrepancy between two labelled contours and is defined as: max max min ( , ) , max min ( , ) g G s S s S g G where G and S are the ground truth and computed contours, respectively, and d(g,s) denotes the Euclidean distance. For the maximum segmentation accuracy, the HD value should be close to 0 (ideally 0). Figure 11 shows the segmentation accuracy comparison using the Matthews correlation coefficient (MCC) from        Table 4. Segmentation accuracy analysis comparison using Jaccard index and Dice coefficient similarity metrics. Table 4 shows the segmentation accuracy comparison of the proposed method with the state-of-the-art using the JI and DSC similarity metrics. Both the mean and standard deviation (mean error) of the evaluated metrics are considered in the result compilation. For Fig. 4, both the proposed method and RSF yield similar values for both JI and DSC. In turn, the proposed method yields best segmentation result for the Fig. 5.
This section also shows segmentation results using 2D brain MR images from a public database of 20 brain anatomical models 36,37 . All images have 250 × 250 pixels and 8 bits per pixel. As a practical application, brain MR images are segmented into white matter (WM) and gray matter (GM) regions, which can be helpful to psychologists to pinpoint psychological diseases and to surgeons during brain surgery.
In order to partition a brain MR image into WM and GM regions, the segmentation result is split into two regions based on two phases: φ > 0 and φ < 0. The WM and GM regions represent the brain region, which is the region of interest, while the regions outside the brain (e.g., skull, fat and vacuum) can be taken as irrelevant regions. Therefore, we manually extracted the brain area to segment the WM and GM regions, removing the other irrelevant regions out of second row, the third and fourth images show the ground truths of the WM and GM regions, respectively. Figure 13 shows the accuracy analysis of the region of interest in the brain MR images. A total of 100 2D slices from 20 brain anatomical models 37 were used. Five 2D slices from every patient were considered. The WM and GM regions for all methods were computed as depicted in Fig. 13. The segmentation accuracy corresponding to the WM and GM regions presented in Fig. 14 was obtained using percentage accuracy in terms of Jaccard index from (17). Figure 14 and Table 5 show that the proposed method yields the best segmentation accuracy in most cases for both the WM and GM regions.
Conclusions and future work. In this paper, a novel edge-based active contour method is proposed that uses a difference of Gaussians (DoG) function as an edge-indicator in its formulation. DoG function uses differences of two smooth images to extract edge information in an image, which acts as a balloon force in the energy functional to evolve the level set curve. In the proposed formulation, the internal energy term penalizes the deviation of the level set function from a signed distance function (SDF) and external energy term evolves the contour towards the boundaries of the objects.
The inclusion of DoG function in the energy formulation work as a global edge extractor and is able to segment all the regions in an image. Moreover, it is also able to properly segment images with intensity inhomogeneity. The results show that proposed method yields best segmentation results using both synthetic and real images as compared to the discussed state-of-the-art methods.
The proposed method segments all regions in an image globally; therefore, it cannot be used in the selective segmentation of a particular region in an image. In future, we would like to formulate a method to segment a selective region in an image, which can be a good help to segment particular regions in medical applications. The proposed method can properly segment intensity inhomogeneous images; however, it is unable to correct the bias of the inhomogeneous regions. In future, we would also like to formulate an edge-based active contours which can also correct the inhomogeneity of the regions.  Table 5. Segmentation accuracy of WM and GM regions (in terms of Jaccard index × 100) using VLSBCS, LSACM and proposed methods.