1 Introduction
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 realworld 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 oversegmentation, i.e., the forming of ”holes” due to image complexity 19 , or undersegmentation, 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 nonlocal shape descriptors 2 ; 3 , multilabel image segmentation with preserved topology 4 , and mincut/maxflow 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 selfintersecting, i.e., merging or splitting. For example, Han et al. 10 proposed a simplepoints 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 nonlocal regularization term, which was applied in the tracking of long thin objects in remote sensing images. Building on the previous ideas, the selfrepelling 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 nonlocal repulsion term to the classic geodesic active contour model (GAC)8 . In the followup 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 nonlocal topologypreserving segmentationguided registration. Attempts have also been made 19 to combine the SR with the regionbased ChanVese model, though a direct combination proved less successful than the original SR.
The SR model has intuitive and straightforward geometric interpretations, but its’ nonlocal 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 Osher
32 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 subproblems 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 subproblems, which turns a secondorder optimization problem into two firstorder ones. Solving the new subproblems no longer requires taking complex differentials of the geodesic curvature term. We also replaced the reinitialization 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 EulerLagrange equations or gradient descent equations for the subproblems. In section 4, the discretization schemes for the subproblems 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 SelfRepelling Snake Model
The original SR model as proposed in 16 is an edgebased 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 selfrepulsion 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
(1) 
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 ,(2) 
The level set function is defined as a signed distance function, such that,
(3) 
where is the Euclidean distance between point and contour . As a signed distance function, satisfies the constraint condition below, i.e. the Eikonal equation,
(4) 
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 ,
(5) 
(6) 
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 edgebased 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
(7) 
where , , are penalty parameters that balance three terms.
(8) 
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).
(9) 
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 .
(10) 
describes the selfrepulsion 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,
(11) 
(12) 
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 through
and 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 selfintersection of the contour.Given the energy functional (7) and the constraint (4) , the variational formulation for SR is
(13) 
and the evolution equation of derived from and is
(14) 
where
(15) 
(15) is the geodesic curvature that shifts the contour towards the edges detected by . In the image areas with nearuniform 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
(16) 
To summarize, by applying variational methods to the three energy terms and substituting with , the following evolution equations can be derived
(17) 
With regards to the constraint , the dynamic reinitialization scheme below is adopted in 16 ,
(18) 
The above is a typical HamiltonJacobi equation that can be discretized and solved through an upwind 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 halfpoint difference scheme and the harmonic averaging approximation. The next two terms adopt the upwind scheme. Two semiimplicit schemes are constructed by concatenating the rows and columns of the image respectively 16 ,
(19) 
where are the two concatenation matrices, and are intermediate variables, and are the upwind 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 ,
(20) 
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
(21) 
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 reinitialization step, more memory efficient by using compact intermediate variables, and more concise by bypassing the complex discretization schemes.
3 The SplitBregman Algorithm for the Selfrepelling 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 reformulate the energy minimization problem as
(22) 
(23) 
where , , and is a penalty parameter. The original problem can then be solved as two subproblems in alternating order for loops . The subproblems are,
(24) 
(25) 
To solve the subproblem (24), we can derive the following evolution equation of via standard variational methods 47 ,
(26) 
The initial condition and boundary condition are as below,
(27) 
where,
(28) 
(29) 
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 edgebased 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 subproblem (25), if , we obtain the corresponding EulerLagrange equation of as,
(30) 
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,
(31) 
(32) 
The first order time derivative of can be approximated as . Therefore, from (26), a semiimplicit iterative scheme can be designed for where , such that,
(33) 
until .
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
(34) 
and . For , since , the iterative formula for from (25) can be written as,
(35) 
(36) 
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
(37) 
The formula for is
(38) 
The same projection scheme as (36) is used afterwards. After have been obtained, we can derive directly from (23).
In summary, the SplitBregman 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 narrowbands 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 oversegmentation) and contour merging (which causes undersegmentation). 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) i77700 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 CourantFriedrichsLewy 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 kmeans clustering or detector filters such as the circle Hough Transform or the Laplacian of Gaussian
45 . 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.
6 Conclusions
By introducing an intermediate variable and the Bregman iterative parameter, the SelfRepelling Snake model can be solved through the SplitBregman method. This new solution bypasses the frequent reinitialization 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.
Acknowledgements
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.
References
 (1) D. MacDonald, N. Kabani, D. Avis, A. C. Evans, Automated 3D extraction of inner and outer surfaces of cerebral cortex from MRI, NeuroImage 12 (3) (2000) 340–356.
 (2) J. A. Geiping, Comparison of topologypreserving segmentation methods and application to mitotic cell tracking, B.S. thesis, Dept. Math. Comput. Sci., Westfälische WilhelmsUniversität at Münster, Münster, Germany (2014).
 (3) H.L. Chan, S. Yan, L.M. Lui, X.C. Tai, Topologypreserving 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 nonlocal topologypreserving segmentationguided 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, Topologypreserving multilabel 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 mincut/maxflow 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.

(12)
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 topologypreserving 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, Selfrepelling snakes for topologypreserving 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 topologypreserving 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 l1regularized 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 curvaturedependent speed: algorithms based on HamiltonJacobi 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 PDEbased 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.
Comments
There are no comments yet.