Using the Split Bregman Algorithm to Solve the Self-repelling Snakes Model

Preserving contour topology during image segmentation is useful in many practical scenarios. By keeping the contours isomorphic, it is possible to prevent over-segmentation and under-segmentation, as well as to adhere to given topologies. The Self-repelling Snakes model (SR) is a variational model that preserves contour topology by combining a non-local repulsion term with the geodesic active contour model. The SR is traditionally solved using the additive operator splitting (AOS) scheme. In our paper, we propose an alternative solution to the SR using the Split Bregman method. Our algorithm breaks the problem down into simpler sub-problems to use lower-order evolution equations and a simple projection scheme rather than re-initialization. The sub-problems can be solved via fast Fourier transform or an approximate soft thresholding formula which maintains stability, shortening the convergence time, and reduces the memory requirement. The Split Bregman and AOS algorithms are compared theoretically and experimentally.


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 real-world brain cortical structure [1].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 [2], 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 [3] and non-local shape descriptors [4,5], multi-label image segmentation with preserved topology [6], and min-cut/max-flow segmentation using topology priors [7].
Since the problem of topology preservation can be intuitively linked to the process of contour evolution, many active contour models [8,9,10] 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. [11] 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) [10].In the follow-up work [17], the short time existence/uniqueness and Lipschitz regularity property of the SR model were studied.Later, [5] successfully extended the SR model to non-local topology-preserving segmentation-guided registration.Attempts have also been made [2] 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) [18] 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 [19] 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) [20] and the Augmented Lagrangian Method (ALM) [21].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.

The Original Self-Repelling Snake Model
The original SR model as proposed in [16] is an edge-based segmentation model based on the GAC [10].It adopts the variational level set formulation [22], where the segmentation contour is implicitly represented as the zero level line of a signed distance function [23].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 f (x) : Ω → R be a scalar value image, x ∈ Ω, and Ω is the domain of the image.The standard edge detect function g(x) ∈ [ 0, 1] is given by where s = 1 or 2, ρ is a scaling parameter, and G σ denotes a Gaussian convolution of the image with a standard deviation of σ.The object boundary C 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 d(x, C) is the Euclidean distance between point x and contour C. 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 H(φ) and Dirac functions δ(φ).Since the original Heaviside function is discontinuous and therefore indifferentiable, we adopt the regularization schemes below [22], 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 E(φ) of the SR model can be written as where γ, α, β are penalty parameters that balance three terms.
E g (φ) 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).
E a (φ) 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 [9] .
measures the nearness of the two points x and y, e.g. the further away the points the smaller the repulsion.In (10), h ε (φ(x)) and h ε (φ(y)) denote the narrow bands around the points x and y, where, When the points x and y are further than distance l from the contour, This signifies that the points outside the narrow bands are largely unaffected by repulsion.For −∇φ(x) • ∇φ(y), if the outwards unit normal vectors to the level lines passing through x and y have opposite directions, i.e., the contours passing x and y are merging or splitting, then the functional approaches the maximum value.Thus, the minimization of E r (φ) 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 φ(x) derived from E g (φ) and E a (φ) is where (15) is the geodesic curvature that shifts the contour towards the edges detected by g(x).In the image areas with near-uniform intensity, ∇g(x) → 0, → 0 in those areas, the geodesic curvature term has little effect and the balloon force αg(x) dominates.
Lastly, the evolution equation derived from the repulsion term is To summarize, by applying variational methods to the three energy terms and substituting δ ε (φ(x)) with |∇φ(x)|, the following evolution equations can be derived With regards to the constraint |∇φ| = 1, 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 [23].To solve (17), the original solution adopts the AOS strategy [18].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 semiimplicit schemes are constructed by concatenating the rows and columns of the image respectively [16], where A x1 , A x2 are the two concatenation matrices, v and w are intermediate variables, and T 2 , T 3 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 A l (l ∈ (x 1 , x 2 )), where i, j are two points in the image, N l (i) is the set of nearest neighbors of , and A l is a diagonally dominant tridiagonal matrix.Finally, φ k+1 can be calculated as Since i and j span the entire image, if Ω ∈ R m×n , then A l ∈ R (m×n)×(m×n) Consequently, the variable A 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 O(N ).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.

The Split-Bregman Algorithm for the Self-repelling Snake Model
The Split Bregman method is a fast alternating directional method often used in solving L 1 -regularized constrained optimization problems [19].To design the Split Bregman algorithm for (7), we first introduce a splitting variable w = ∇φ and the Bregman iterator b.We can re-formulate the energy minimization problem as where b 0 = 0, w 0 = 0, and µ is a penalty parameter.The original problem can then be solved as two sub-problems in alternating order for loops k = 1, 2, ..., K.
The sub-problems are, To solve the sub-problem (24), we can derive the following evolution equation of φ via standard variational methods [24], The initial condition and boundary condition are as below, where, With the Heaviside function originally adopted in [16], the newly introduced component δ ε (φ) in the Split Bregman algorithm may be excessively smoothed.
However, since the second term in (30) contains the integral of w(y), it is difficult to construct the iterative scheme for w k .An approximation formula with projection is designed in the next section to address this issue.

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 (i, j) be x i,j = (x 1i , x 2j ) where i = 0, 1, 2, ..., M + 1, j = 0, 1, 2, ..., N + 1 , we get φ i,j = φ(x 1i , x 2j ).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 φ i,j can be approximated as ∂φi,j ∂t . Therefore, from (26), a semi-implicit iterative scheme can be de-signed for φ k+1,s+1 i,j where s = 0, 1, 2, ..., S, such that, until y denotes a point taken from a small window of size 2d × 2d around point x.
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.
In practice, a single iteration is often enough for computing (35).Alternatively, we can directly use the soft thresholding formula to derive w k+1 .For abbreviation, let The formula for w k+1 i,j is The same projection scheme as (36) is used afterwards.After φ k+1 i,j , w k+1 i,j have been obtained, we can derive b k+1 i,j 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 m×n, the parameter A in the AOS solution is size 2×(m×n)×(m×n).However, in the Split Bregman algorithm, the sizes of both w and b are 2 × (m × n) 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 w 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.

Experimental Results
The experiments below demonstrate that the Split Bregman solution of the SR model successfully prevents contour splitting (which causes over-segmentation)   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 5 × 5 or 7 × 7 as according to [16].The time step τ is chosen according to the convergence condition τ ≤ 1 4µ based on the Courant-Friedrichs-Lewy condition [25].
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.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 Gaussian [27].Since the topology is maintained, the number of cells remains the same.In   The algorithm can also be extended to 3D, as seen in Figure 5.

Conclusions
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.

Furthermore, as
the SR is an edge-based model and the repelling force is local, smoothing H(φ) 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 2ε surrounding the contour which in practice stabilizes contour evolution.For the sub-problem (25), if | w(x)| = 0, we obtain the corresponding Euler-Lagrange equation of w(x) as,

Algorithm 1 :
The Split Bregman algorithm for the Self-repelling Snake Model (1) Initialize Calculate g(x) using (1) Initialize φ 0 (x) as signed distance function and set w 0 = For s=0,1,2,...,S Calculate φ k+1,s+1 from (33) End for s when (24) converges Calculate w k+1 from (38) Calculate b k+1 from (23) End for k when (13) converges 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.

Figure 3
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.

Figure 4 ,
Figure 4, the repulsive force prevents cell contours from merging and separates the adhesive cells.