Topology preservation in image segmentation is an external constraint to discourage changes in the topology of the segmentation contour. It is usually applied in problems where the object topology is known a priori. One example is in medical image analysis where the segmentation of the brain cortical surface must produce results consistent with the real-world brain cortical structure 41 . Another example is the segmentation of objects with complicated interiors, noises, or occlusions, where a topological constraint can be used to prevent over-segmentation, i.e., the forming of ”holes” due to image complexity 19 , or under-segmentation, i.e., when the contours of separate objects merge. Much active research undergoes in the area, such as image segmentation and registration using the Beltrami representation of shapes 1 and non-local shape descriptors 2 ; 3 , multi-label image segmentation with preserved topology 4 , and min-cut/max-flow segmentation using topology priors 5 .
Since the problem of topology preservation can be intuitively linked to the process of contour evolution, many active contour models 6 ; 7 ; 8 have been proposed for it. In these models, the contour is affected by various forces until it converges to the final segmentation result. To preserve topology during the contour evolution process, a constraint term is usually added to the variational formulation which prevents the contour from self-intersecting, i.e., merging or splitting. For example, Han et al. 10 proposed a simple-points detection scheme in an implicit level set framework in 2003. Meanwhile, Cecil et al. 12 monitored the changes in the Jacobian of the level set. In 2005, Alexandrov et al. 13 recast the topology preservation problem to a shape optimization problem of the level sets, where narrow bands around the segmentation contours are discouraged from overlapping. Sundaramoorthi and Yezzi 14 proposed an approach based on knot energy minimization to realize the same effect. Rochery et al. 15 used a similar idea while introducing a non-local regularization term, which was applied in the tracking of long thin objects in remote sensing images. Building on the previous ideas, the self-repelling snake model (SR) was proposed by Le Guyader et al. in 2008 16 . The SR uses an implicit level set representation for the curve and adds a non-local repulsion term to the classic geodesic active contour model (GAC)8 . In the follow-up work 17 , the short time existence/uniqueness and Lipschitz regularity property of the SR model were studied. Later, 3 successfully extended the SR model to non-local topology-preserving segmentation-guided registration. Attempts have also been made 19 to combine the SR with the region-based Chan-Vese model, though a direct combination proved less successful than the original SR.
The SR model has intuitive and straightforward geometric interpretations, but its’ non-local term leads to complications in the numerical implementation. To the best of our knowledge, the SR model has always been solved through the additive operator splitting (AOS) 48 strategy. On the one hand, the derivation of gradient descent equations is complicated and requires the upwind difference discretisation scheme. On the other hand, though the AOS is stable, the memory requirement grows quadratically with the size of the image. In this paper, we propose an alternative solution using the Split Bregman method that aims towards a more concise algorithm and less memory usage.
The Split Bregman algorithm was first proposed in computer vision by Goldstein and Osher32 for the total variation model (TV) for image restoration. By introducing splitting variables and iterative parameters, it transforms the original constrained minimization problems into simpler sub-problems that can be solved alternatively. The Split Bregman algorithm is shown to be equivalent to the Alternating Direction Method of Multipliers (ADMM) 33 and the Augmented Lagrangian Method (ALM) 34 . In this paper, we introduce an intermediate variable to split the original problem into two sub-problems, which turns a second-order optimization problem into two first-order ones. Solving the new sub-problems no longer requires taking complex differentials of the geodesic curvature term. We also replaced the re-initialization of the signed distance function with a simple projection scheme. As a result, the optimization of the level set function is simplified. In addition, to address some problems arising from the Split Bregman solution, we replaced the Heaviside representation of the level set in 16 with one that performed better in our algorithm.
The paper is organized as follows. In section 2, we review and provide some intuition to the original SR model. In section 3, we design the Split Bregman algorithm for the SR model and derive the Euler-Lagrange equations or gradient descent equations for the sub-problems. In section 4, the discretization schemes for the sub-problems are presented for the alternating iterative optimization. In section 5, we provide some numerical examples and comparisons of results. Finally, we draw conclusions in section 6.
2 The Original Self-Repelling Snake Model
The original SR model as proposed in 16 is an edge-based segmentation model based on the GAC 8 . It adopts the variational level set formulation 42 , where the segmentation contour is implicitly represented as the zero level line of a signed distance function 11 . An energy functional is minimized until convergence is reached and the segmentation contour is obtained. The energy functional comprises three terms, two of which are taken from the GAC model and contribute to edge detection and the balloon force respectively, while the last one accounts for the self-repulsion of contour as it approaches itself.
The definition of the SR model is as follows. Let be a scalar value image, , and is the domain of the image. The standard edge detect function is given by
where or 2, is a scaling parameter, and
denotes a Gaussian convolution of the image with a standard deviation of. The object boundary is represented by the zero set of a level set function ,
The level set function is defined as a signed distance function, such that,
where is the Euclidean distance between point and contour . As a signed distance function, satisfies the constraint condition below, i.e. the Eikonal equation,
To represent the image area and contour, we introduce the Heaviside function and Dirac functions . Since the original Heaviside function is discontinuous and therefore indifferentiable, we adopt the regularization schemes below 42 ,
These schemes are different from the ones chosen for the original model in 16 . In particular, does not regularize the entire image domain, which improves stability of edge-based models. The effect is more apparent in our Split Bregman algorithm, as we will discuss further in section 3.
Given the above, the energy functional of the SR model can be written as
where , , are penalty parameters that balance three terms.
is the geodesic length of the contour. The total variation of the Heaviside function, or the total length of the contour, is weighted by the edge detector in (1).
is the closed area of the contour also weighted by the edge detector. It contributes to a balloon force that pushes the segmentation contour over weak edges 7 .
describes the self-repulsion of the contour 16 . measures the nearness of the two points and , e.g. the further away the points the smaller the repulsion. In (10), and denote the narrow bands around the points and , where,
When the points and are further than distance from the contour, . This signifies that the points outside the narrow bands are largely unaffected by repulsion. For
, if the outwards unit normal vectors to the level lines passing throughand have opposite directions, i.e., the contours passing and are merging or splitting, then the functional approaches the maximum value. Thus, the minimization of prevents the self-intersection of the contour.
Given the energy functional (7) and the constraint (4) , the variational formulation for SR is
and the evolution equation of derived from and is
(15) is the geodesic curvature that shifts the contour towards the edges detected by . In the image areas with near-uniform intensity, , . Since in those areas, the geodesic curvature term has little effect and the balloon force dominates.
Lastly, the evolution equation derived from the repulsion term is
To summarize, by applying variational methods to the three energy terms and substituting with , the following evolution equations can be derived
With regards to the constraint , the dynamic re-initialization scheme below is adopted in 16 ,
The above is a typical Hamilton-Jacobi equation that can be discretized and solved through an up-wind difference scheme 11 . To solve (17), the original solution adopts the AOS strategy 48 . The first term on the r.h.s. of (17) is discretized with the half-point difference scheme and the harmonic averaging approximation. The next two terms adopt the up-wind scheme. Two semi-implicit schemes are constructed by concatenating the rows and columns of the image respectively 16 ,
where are the two concatenation matrices, and are intermediate variables, and are the up-wind discretizations of the second and third term of the r.h.s. of (17), the formulations of which are omitted here for brevity. For each ,
where are two points in the image, is the set of nearest neighbors of in the matrix , , and is a diagonally dominant tridiagonal matrix. Finally, can be calculated as
Since and span the entire image, if , then Consequently, the variable greatly increases the memory requirement for the AOS solution. In the last step, (19) is solved via the Thomas algorithm which involves LR decomposition, forward substitution, and backward substitution, with the convergence rate of . In the following section, we will propose another solution to the SR with the Split Bregman method that aims to be faster by replacing the re-initialization step, more memory efficient by using compact intermediate variables, and more concise by bypassing the complex discretization schemes.
3 The Split-Bregman Algorithm for the Self-repelling Snake Model
The Split Bregman method is a fast alternating directional method often used in solving -regularized constrained optimization problems 32 . To design the Split Bregman algorithm for (7), we first introduce a splitting variable and the Bregman iterator . We can re-formulate the energy minimization problem as
where , , and is a penalty parameter. The original problem can then be solved as two sub-problems in alternating order for loops . The sub-problems are,
To solve the sub-problem (24), we can derive the following evolution equation of via standard variational methods 47 ,
The initial condition and boundary condition are as below,
With the Heaviside function originally adopted in 16 , the newly introduced component in the Split Bregman algorithm may be excessively smoothed. Furthermore, as the SR is an edge-based model and the repelling force is local, smoothing over the entire image causes the repelling force to propagate across the image, resulting in unnecessary instability. With the new choice of Heaviside function, the smoothing effect is restricted only to a narrow band of width surrounding the contour which in practice stabilizes contour evolution.
For the sub-problem (25), if , we obtain the corresponding Euler-Lagrange equation of as,
However, since the second term in (30) contains the integral of , it is difficult to construct the iterative scheme for . An approximation formula with projection is designed in the next section to address this issue.
4 Discretization and Iterative Scheme
For the next step in solving (26) and (30), we devise the discretization of the continuous derivatives. Let the spatial step be 1 and time step be , and the discrete coordinates for the pixel be where , , we get . Let the other variables take similar forms. With the first order finite difference approximation, we can obtain the discrete gradient, Laplacian, and divergences respectively as,
The first order time derivative of can be approximated as . Therefore, from (26), a semi-implicit iterative scheme can be designed for where , such that,
which is the discrete approximation of . denotes a point taken from a small window of size around point . The repulsion from points further away is negligible, therefore we only check within a small window. Note that the initial and boundary conditions in (27) still hold.
Next, we will solve (30) iteratively. By temporarily fixing , we can design a concise approximate generalized soft thresholding formula. For abbreviation, let
and . For , since , the iterative formula for from (25) can be written as,
In practice, a single iteration is often enough for computing (35). Alternatively, we can directly use the soft thresholding formula to derive . For abbreviation, let
The formula for is
The same projection scheme as (36) is used afterwards. After have been obtained, we can derive directly from (23).
In summary, the Split-Bregman algorithm proposed in this section has four main advantages. First, the memory requirement is reduced. For an image of size , the parameter in the AOS solution is size . However, in the Split Bregman algorithm, the sizes of both and are only. As the image size increases, the memory usage in the original algorithm increases quadratically while the one in the new algorithm increases linearly. This is an important point when dealing with big images. Second, the numerical solution can be simplified. In (17), the convolution term containing is hyperbolic, which requires the upwind difference scheme. By substituting with the auxiliary variable we can remove the need for complex discretization schemes. Third, the use of a simple projection scheme in place of the initialization step improves algorithm efficiency. Finally, contour evolution is stabilized by confining the smoothing of the Heaviside function to the narrow-bands around the contours.
The proposed algorithm is summarized in Algorithm 1.
5 Numerical Experiments
5.1 Experimental Results
The experiments below demonstrate that the Split Bregman solution of the SR model successfully prevents contour splitting (which causes over-segmentation) and contour merging (which causes under-segmentation). The qualitative performance is comparable to the original solution of SR. Two piratical applications are showcased as well as the adaptation to 3D. All experiments are performed on the PC (Intel(R) Core (TM) i7-7700 CPU @ 3.60GHz 3.60 GHz; 16.0 GB memory). The segmentation program is written in Matlab and runs in Matlab environment R2018b.
In Figure 1, contour splitting is successfully prevented and the topology is preserved. The parameter controls the outwards or inwards driving force, dictates the geodesic length, weighs the repelling force, and weighs the constraint. A large causes the contour to become unstable, as the repulsive force is a highly local term. However, increasing and decreasing the window size narrows the gap between the contours. Typically, the window size is or as according to 16 . The time step is chosen according to the convergence condition based on the Courant-Friedrichs-Lewy condition 44 . Increasing improves the smoothness of the contour but lowers the effectiveness of topology preservation, as it smooths out the repulsive force.
In Figure 2, contour merging is prevented as the fingers of the hand remain separate. In the basic GAC model, the proximity of the contours would cause them to merge despite there being a detected edge, because it reduces the total geodesic length.
Figure 3 compares the effect of two different Heaviside functions. The advantage of the new Heaviside formulation lies in the stabilization of the repelling term, which makes the algorithm more robust. In the experiment in Figure 3, the amount of repulsion was set to very high through . Yet, it still did not disturb the smoothness of the contour nor cause the lose of necessary details. With the same set of parameters and the original Heaviside function, the repulsive force was dissipated and stability could not maintained.
One example of a practical application of the algorithm is adhesive cell segmentation. The centers of cells can be detected via k-means clustering or detector filters such as the circle Hough Transform or the Laplacian of Gaussian45 . Since the topology is maintained, the number of cells remains the same. In Figure 4, the repulsive force prevents cell contours from merging and separates the adhesive cells.
The algorithm can also be extended to 3D, as seen in Figure 5.
By introducing an intermediate variable and the Bregman iterative parameter, the Self-Repelling Snake model can be solved through the Split-Bregman method. This new solution bypasses the frequent re-initialization of the signed distance function and simplifies computation, as well as reduces the memeory requirement. The performance of the Split Bregman solution is similar to that of the original solution in 16 , e.g. both merging and splitting of the segmentation contour can be prevented and the topology is preserved.
Special thanks to Dr. C. Le Guyader (Université de Rouen, France) for sharing her source code for the AOS solution of SR. Many thanks to the reviewers for their valuable comments and suggestions.
- (1) D. MacDonald, N. Kabani, D. Avis, A. C. Evans, Automated 3-D extraction of inner and outer surfaces of cerebral cortex from MRI, NeuroImage 12 (3) (2000) 340–356.
- (2) J. A. Geiping, Comparison of topology-preserving segmentation methods and application to mitotic cell tracking, B.S. thesis, Dept. Math. Comput. Sci., Westfälische Wilhelms-Universität at Münster, Münster, Germany (2014).
- (3) H.-L. Chan, S. Yan, L.-M. Lui, X.-C. Tai, Topology-preserving image segmentation by Beltrami representation of shapes, Journal of Mathematical Imaging and Vision 60 (3) (2018) 401–421.
- (4) N. Debroux, C. Le Guyader, A joint segmentation/registration model based on a nonlocal characterization of weighted total variation and nonlocal shape descriptors, SIAM Journal on Imaging Sciences 11 (2) (2018) 957–990.
- (5) N. Debroux, S. Ozeré, C. Le Guyader, A non-local topology-preserving segmentation-guided registration model, Journal of Mathematical Imaging and Vision 59 (3) (2017) 432–455.
- (6) J. Waggoner, Y. Zhou, J. Simmons, M. De Graef, S. Wang, Topology-preserving multi-label image segmentation, in: 2015 IEEE Winter Conference on Applications of Computer Vision, IEEE, 2015, pp. 1084–1091.
- (7) Y. Zeng, D. Samaras, W. Chen, Q. Peng, Topology cuts: A novel min-cut/max-flow algorithm for topology preserving segmentation in N–D images, Computer vision and image understanding 112 (1) (2008) 81–90.
- (8) M. Kass, A. Witkin, D. Terzopoulos, Snakes: Active contour models, International journal of computer vision 1 (4) (1988) 321–331.
- (9) L. D. Cohen, On active contour models and balloons, CVGIP: Image understanding 53 (2) (1991) 211–218.
- (10) V. Caselles, R. Kimmel, G. Sapiro, Geodesic active contours, International journal of computer vision 22 (1) (1997) 61–79.
- (11) X. Han, C. Xu, J. L. Prince, A topology preserving level set method for geometric deformable models, IEEE Transactions on Pattern Analysis and Machine Intelligence 25 (6) (2003) 755–768.
T. C. Cecil, Numerical methods for partial differential equations involving discontinuities, Ph.D. thesis, University of California, Los Angeles (2003).
- (13) O. Alexandrov, F. Santosa, A topology-preserving level set method for shape optimization, Journal of Computational Physics 204 (1) (2005) 121–130.
- (14) G. Sundaramoorthi, A. Yezzi, Global regularizing flows with topology preservation for active contours and polygons, IEEE Transactions on Image Processing 16 (3) (2007) 803–812.
- (15) M. Rochery, I. H. Jermyn, J. Zerubia, Higher order active contours, International Journal of Computer Vision 69 (1) (2006) 27–42.
- (16) C. Le Guyader, L. A. Vese, Self-repelling snakes for topology-preserving segmentation models, IEEE Transactions on Image Processing 17 (5) (2008) 767–779.
- (17) N. Forcadel, C. Le Guyader, A short time existence/uniqueness result for a nonlocal topology-preserving segmentation model, Journal of Differential Equations 253 (3) (2012) 977–995.
- (18) J. Weickert, B. T. H. Romeny, M. A. Viergever, Efficient and reliable schemes for nonlinear diffusion filtering, IEEE transactions on image processing 7 (3) (1998) 398–410.
- (19) T. Goldstein, S. Osher, The split bregman method for l1-regularized problems, SIAM journal on imaging sciences 2 (2) (2009) 323–343.
- (20) T. Goldstein, B. O’Donoghue, S. Setzer, R. Baraniuk, Fast alternating direction optimization methods, SIAM Journal on Imaging Sciences 7 (3) (2014) 1588–1623.
- (21) C. Wu, X.-C. Tai, Augmented Lagrangian method, dual methods, and split Bregman iteration for ROF, vectorial TV, and high order models, SIAM Journal on Imaging Sciences 3 (3) (2010) 300–339.
- (22) H.-K. Zhao, T. Chan, B. Merriman, S. Osher, A variational level set approach to multiphase motion, Journal of computational physics 127 (1) (1996) 179–195.
- (23) S. Osher, J. A. Sethian, Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations, Journal of computational physics 79 (1) (1988) 12–49.
- (24) T. F. Chan, J. J. Shen, Image processing and analysis: variational, PDE, wavelet, and stochastic methods, Vol. 94, Siam, 2005.
- (25) S. Osher, R. Fedkiw, Level set methods and dynamic implicit surfaces, Vol. 153, Springer Science & Business Media, 2006.
- (26) J. Duan, Variational and PDE-based methods for image processing, Ph.D. thesis, University of Nottingham (2018).
- (27) H. Kong, H. C. Akakin, S. E. Sarma, A generalized Laplacian of Gaussian filter for blob detection and its applications, IEEE transactions on cybernetics 43 (6) (2013) 1719–1733.