Multipass Active Contours for an Adaptive Contour Map

Isocontour mapping is efficient for extracting meaningful information from a biomedical image in a topographic analysis. Isocontour extraction from real world medical images is difficult due to noise and other factors. As such, adaptive selection of contour generation parameters is needed. This paper proposes an algorithm for generating an adaptive contour map that is spatially adjusted. It is based on the modified active contour model, which imposes successive spatial constraints on the image domain. The adaptability of the proposed algorithm is governed by the energy term of the model. This work focuses on mammograms and the analysis of their intensity. Our algorithm employs the Mumford-Shah energy functional, which considers an image's intensity distribution. In mammograms, the brighter regions generally contain significant information. Our approach exploits this characteristic to address the initialization and local optimum problems of the active contour model. Our algorithm starts from the darkest region; therefore, local optima encountered during the evolution of contours are populated in less important regions, and the important brighter regions are reserved for later stages. For an unrestricted initial contour, our algorithm adopts an existing technique without re-initialization. To assess its effectiveness and robustness, the proposed algorithm was tested on a set of mammograms.


Introduction
The extraction of meaningful information from images by means of digital image processing techniques is an important task in many application domains. To identify significant information in an image, one can exploit distinctive features of objects (e.g., shape, location, margin) and uncover significant regions or patterns by analyzing the topological and geometrical properties of the image. In the biomedical sector, medical imaging techniques, such as X-rays, magnetic resonance imaging (MRI), and tomography are used to visualize internal structures of the body. Medical images, such as mammograms, have inherently complex and variable features with blurred object boundaries, which make the use of explicit features of objects in image analysis difficult. Hence, image analysis methods based on isocontour mapping are better suited to complex medical images as mentioned below.
Meaningful information can be extracted efficiently from a digital image on an isocontour map. In biomedical image processing, isocontour mapping is used extensively to perform the topographic analysis of medical images. An isocontour map consists of a set of curves of equal value (e.g., height or intensity). Analyses based on isocontour maps can provide the association between image inclusions. They detect a region of interest (ROI) by analyzing the correlation (or the enclosure relationship) between objects. ROI analysis that highlights suspicious regions in medical images is an essential step in computer-aided diagnosis systems. Thus, isocontour maps can provide a robust topographic representation of medical images for ROI analysis [1].
To our current knowledge there is no deterministic way to find parameters like the number of quantization levels (contour interval and the difference in elevation between successive contours) that yield the best results for isocontour generation. However, the determination of contour generation parameters is an open question and adjustable. Hong and Sohn [1] proposed a multiscale approach for ROI segmentation, which extracts isocontours at multiple scales and analyzes mammographic features in a hierarchical manner from a coarse scale to a fine one. This multiscale approach was necessary because the information provided by isocontour maps with fixed parameters was sometimes either too excessive or scarce due to varying image conditions. This paper aims to produce an adaptive contour map that provides "not too much and not too little" information by adapting active contours spatially during the curve evolution.
A number of active contour models have been developed. Kass et al. [2] proposed a successful method based on variational and partial differential equations (PDE), the well known active contour/snake model, to extract interesting objects in an image. Various active contour models and enhanced versions are employed in various image processing applications, as well as medical images. The active contours are represented as parameterized curves in a Lagrangian framework [2] and the implicit curves are given in an Eulerian framework [3][4][5][6].
Geodesic active contour (GAC) models in [3,4] are geometrically intrinsic and embed the level set function [7], which involves the representation of the implicit curve. The curve evolution with the level set function naturally splits and merges the contours during the evolution, and therefore automatically handles topological changes. The curves evolve based on the minimization of the energy functional from the image, the curve, and the level set function. Energy functionals are used in the energy of edge-based model [2][3][4][5] and the region-based model [6].
The classical active contour model [2][3][4][5][6][7][8], which detects objects in an image, starts with a given initial contour and performs the curve evolution to find the optimal contour. In the algorithm for adaptive contour mapping proposed in this paper, the initial contour divides the image domain into sub-regions in which a new optimal contour is found. In subsequent iterations, the contour would have a different spatial domain from that of the previous contour. This domain segmentation (or the curve evolution) is repeated until the stopping criterion is met, thereby creating an adaptive contour map of the image. The adaptability of the proposed algorithm is governed by the energy term of the active contour model, so it is important to employ one that is both effective and reliable.
The proposed algorithm for adaptive contour mapping is based on two previous active contour models: active contours without edges (ACWE) and level set evolution without re-initialization (LSEWR). From the ACWE concept, we used force image terms to get regional information, whereas in LSEWR we selected penalizing terms to eliminate re-initialization. The ACWE model proposed by Chan and Vese [6] considers the intensity distribution of an image to establish an optimality criterion for segmenting the image into sub-regions; therefore, it is suitable for use in analyzing mammographic intensities. The ACWE model finds an optimal partition from the energy of a region in an image that has a weak edge and heavy noise. The ACWE method converges relatively faster than edge based active contours [2][3][4][5] because the merging of similar regions occurs broadly while contours move narrowly. The level set function partitions a region into an inside and outside of a zero level curve. An extension of the ACWE model is proposed for a multi-phase segmentation. Vese and Chan [8] proposed the multi-phase segmentation model with n level set functions. This method always presents 2 n regions from the combination of each phase by level set functions. The curve evolution with the level set function requires costly re-initialization because the level set function deviates from a signed distance function (SDF) in each evolution. Li et al. [5] proposed the LSEWR model, which consists of an internal energy term that penalizes the deviation of the level set function from an SDF, and thus eliminates re-initialization.
Our algorithm is designed with a similar manner of isocontour mapping to detect an arbitrary number of contours for spatial adaptive isocontour mapping. The existing multi-phase method, which detects contours at multiple level sets, always produces 2 n regions. This indicates that many insignificant features might be included in the contour map, thereby influencing the image analysis results. Our approach divides a region into two sub-regions using the base contour. It then divides one of the segmented sub-regions into two sub-regions in successive iterations. The proposed algorithm detects sub-regions by minimizing the new energy model, restricting it to the characteristic function of a base sub-region. The iterative segmentation process automatically terminates when the stopping criterion is met. Note that only one of the two sub-regions is further segmented in successive iterations. This is associated with the characteristics of the mammographic image in addition to problems in initialization and local optimum of the active contour model.
In mammograms, bright regions contain information that is more significant (e.g., candidate masses). Our algorithm takes advantage of this mammographic characteristic to address the problems in initialization and local optimum of the active contour model. That is, the proposed algorithm starts with the initial contour found in the darkest (low intensity) region so that the local optima encountered during the contour evolution is placed in less important low intensity regions. In terms of initialization, this enables the our approach to start its operation with an initial contour in the darkest region of the

Active Contours
Active contours are a core component of computer vision and medical imaging, and they can be divided into three areas: boundary driven, region driven, and hybrid contours that combine the boundary and region driven areas.
Malladi et al. [4] presented an approach to shape modeling in which a speed term synthesized from the image is used to stop the contour near object boundaries. They used geometric flows for boundary extraction in which distance transforms were used as embedding functions. Caselles et al. [3] proposed a geodesic approach for object segmentation, which allows connecting classical "snakes" based on energy minimization and geometric active contours based on the theory of curve evolution. Paragios et al. [9] proposed an edge-driven bidirectional geometric flow for boundary extraction. They combined the geodesic active contour flow and the gradient vector flow external force for snakes. The resulting motion equation is considered within a level set formulation that can deal with topological changes and important shape deformations. Weickert et al. [10] presented fast algorithms based on the semi-implicit additive operator splitting (AOS) scheme for both geometric and geodesic active contour models.
In 2001, Chan and Vese [6] proposed a new model for active contours to detect objects in a given image, based on the techniques of curve evolution, Mumford-Shah functional for segmentation, and level sets. Their model can detect objects whose boundaries are not necessarily defined by gradients. We minimize an energy that represents a particular case of the minimal partition problem. In the level set formulation, the minimal partition problem becomes evolving the active contour, which will stop at the desired boundary. However, unlike the classical active contour models, the stopping term does not depend on the gradient of an image, but instead is related to its homogeneity. Paragios et al. [11] introduced a frame partition paradigm within the level set space, which can account for boundary and global region-driven information. In 2005, Li et al. [5] presented a new variational formulation for geometric active contours, which forces the level set function to be close to a signed distance function, thus eliminating the need for the costly re-initialization procedure.
Vese and Chan [8] proposed a new multi-phase level set framework for image segmentation using the Mumford and Shah model for application in piecewise constant and piecewise smooth optimal approximations. The proposed method is also a generalization of an active contour model without edge-based two-phase segmentation. They introduced classifications according to a combination of all level sets at a given pixel. Cremers et al. [12] presented a novel variational approach for segmenting the image plane into a set of regions of piecewise constant motion based on only two consecutive frames from an image sequence. They proposed the implementation of this functional using a multi-phase level set framework. Minimizing the functional with respect to its dynamic variables results in an evolution equation for a vector-valued level set function and an eigenvalue problem for motion vectors.
Our algorithm is performed in such a way that it partitions an image into two sub-regions. One of the sub-regions is then iteratively partitioned into two more sub-regions. The sub-regions are determined by the minimization of a new energy model restricted to a characteristic function of a sub-region, and no re-initialization is needed. Our method segments the image into any number of regions, and the process automatically terminates at the stationary solution. In this paper, the proposed approach is The energy of sub-region k w only affects E by the region restriction function k M , which is calculated from k ϕ . The weighted arc length functional L g with univariate Dirac Delta function H δ ′ = in the level set method [7] is: where the edge indicator function g is a positive and strict decreasing function. For example: The zero level curve C is driven into a smooth curve from a complicated curve to minimize the functional L g . By the level set method, the weighted area functional A g is expressed as: The small energy of A g accelerates the curve evolution. By definition, an SDF satisfies a desirable property | | 1 I ∇ = . Li et al. [5] proposed the functional: as a metric. The functional P characterizes how close a function φ is to an SDF. The energy functional ( ) g A φ in Equation (10) is introduced to speed up the curve evolution. Note that when the function g is constant 1, the energy functional in Equation (10)  E c c φ in Equation (4) can be written as the gradient flow: where Δ is the Laplacian operator. Therefore, the function φ , which minimizes this functional, satisfies the Euler Lagrange equation / 0 This gradient flow is the evolution equation of the level set function in the proposed method. The first and second terms on the right hand side of Equation (12) correspond to the energy functional F 1 and F 2 , which partitions a region inside and outside of the zero curve C based on the energy functional values. The third and the fourth terms correspond to the gradient flows of the energy functionals ( ) g L μ φ and ( ) g vA φ , respectively, and are responsible for driving the zero level curve towards the object boundaries. The fifth term, which is associated with the penalizing energy ( ) P α φ , represents the gradient flow: which has the factor ( ) A classic iterative process to minimize the functional E is the following gradient flow with an artificial time t: If one of the sub-regions 1 k w + or 1 k w + is empty, then the formulation degenerates, and therefore the algorithm automatically terminates. Finally, the principal steps of the algorithm are: (7) and initialize φ by 0 φ • Solve the PDE for φ with Equation (14) to obtain 1 k ϕ + • Check whether the solution is stationary. If not, 1 k k = + and repeat from step 2 In this paper, we note that the stationary problem obtained directly from the minimization problem could also be solved numerically using a similar finite difference scheme.
where d ε is a small constant value to prevent dividing-by-zero and , x y D D are difference operators. (12) is approximated by: The Laplacian operator has been implemented in a similar manner.

Selection of Time Steps and Other Constants
In our experiments, we choose the following parameters: τ = where τ is a time step in the numerical implementation of Equation (14). We know that the time step τ and the coefficient α must satisfy 1/ 4 τα < to maintain stable level set evolution [5]. Using a larger time step can speed up the curve evolution, but may cause errors in the boundary location if the time step chosen is too large. There is a tradeoff between choosing a larger time step and accuracy in boundary location. In our case, we use 0.2 τα = .

Initialization of Level Set Function
In outmoded level set methods, it is essential to initialize the level set function φ as an SDF 0 φ . If the initial level set function is expressively different from a signed distance function, then the re-initialization schemes are not able to re-initialize the function to a signed distance function. In our formulation, not only is the re-initialization procedure eliminated, but the level set function φ also no longer requires initialization as a signed distance function by the penalizing energy in [5]. We propose the following functions as the initial level set function 0 φ , where the denoised image ( ) x is the calculated value of mask for every calculated 1 k ϕ + value, that is, for every evolved level set function. Performing from the lowest intensity region, the initial level set function 0 φ is defined as: where 0 ρ > is a constant and 0 κ > is a small constant. In this section, we use 0 φ with 4 ρ ε = and 1 κ = . By definition, this initial level set function 0 φ takes only two values: 4ε − and 4ε .

Results and Discussion
The proposed multi-phase segmentation algorithm has been applied to synthetic and real mammographic images from the mini-MIAS database [16]; the range of intensity in all images is represented from 0 to 255 and the images are 1,024 × 1,024 pixels in size. We show the approximated circle enclosing the abnormality from the database, and the two-phase segmentations in each pass k of Equation (3) with a number of iterations of Equation (14).
In Figures 4 and 5, we illustrate the result of real mammographic images and the approximated region of abnormality from the database. As represented, the proposed algorithm partitions the images, including weak and blurred edges. The recursive segmentation on the higher intensity region finely segments the region. In Figure 5, we consider a more difficult case, which is an abnormality region in a high intensity region. Our algorithm reduces the number of contours in the map from an average of 206 contours to 11 contours.  Figure 6 shows the segmentation results of each pass on a synthetic image with 5% uniform noise. The red regions represent the inside partition, and the other regions indicate the outside partition. Each two-phase segmentation pass is converged after a number of iterations of Equation (14). The proposed algorithm works well on sharp-edged objects.

Conclusions
This paper presents an adaptive contour map that captures topographic information in mammograms characterized by blurred object boundaries. The proposed multipass active contour algorithm for adaptive contour mapping is based on the two-phase piecewise constant segmentation model (ACWE) proposed by Chan and Vese [6] and the variational level set formulation of curve evolution without re-initialization (LSEWR) proposed by Li et al [5]. Our algorithm performs spatial segmentation iteratively to find the optimal contours. It starts with the initial contour selected in the darkest region of a mammographic image so that the brighter regions of the image can be segmented and analyzed in a more refined way at later stages. The proposed algorithm provides an optimized topographic representation of mammograms that can increase the computational efficiency and accuracy of the analysis.
Our algorithm produces an arbitrary number of regions, and it automatically terminates when its stopping condition is met. The proposed algorithm was tested using synthetic and real mammographic images that include masses varying in size and subtlety. The experimental results showed that our approach yields an accurate contour map of both distinctive and subtle masses in mammograms. It also successfully produced an adaptive contour map in synthetic images that have relatively clear edges.
The experimental results show sensitive segmentation on the important region as well as intuitive segmentation structure.