Abstract-The geodesic active contour model (GAC) is a commonly used segmentation model for medical image segmentation. The level set method (LSM) is the most popular approach for solving the model, via implicitly representing the contour by a level set function. However, the LSM suffers from high computation burden and numerical instability, requiring additional regularization terms or re-initialization techniques. In this paper, we use characteristic functions to implicitly approximate the contours, propose a new representation to the GAC and derive an efficient algorithm termed as the iterative convolutionthresholding method (ICTM). Compared to the LSM, the ICTM is simpler and much more efficient and stable. In addition, the ICTM enjoys most desired features (e.g., topological changes) of the level set-based methods. Extensive experiments, on 2D synthetic, 2D ultrasound, 3D CT, and 3D MR images for nodule, organ and lesion segmentation, demonstrate that the ICTM not only obtains comparable or even better segmentation results (compared to the LSM) but also achieves dozens or hundreds of times acceleration.
Index Terms-Geodesic active contours, segmentation, convolution, thresholding
A CTIVE contours have been widely used in various segmentation tasks [1] and image modalities [2] , such as organs [3] in magnetic resonance (MR) scans, tumors in computed tomography (CT) scans ( [4] , [5] ), and ultrasound images ( [6] , [7] ). Basically, there are mainly two types of active contour models: edge-based active contours (e.g., [8] , [9] , [10] ) and region-based active contours (e.g., [11] , [12] , [13] , [14] ). Edge-based active contours are driven by edge indicator functions that are commonly defined by image gradients. The contour evolution is expected to stop on boundaries with high gradient magnitude. Region-based active contours are driven by regional information that can be defined by intensity statistical information inside and outside of the contour. The contour is expected to evolve to a position where the regional information inside and outside of the contour reaches a balance.
This paper focuses on edge-based active contours that were firstly proposed by Kass et al. [15] in 1988 (also termed as the Manuscript received July 1, 2020. Jun Ma is with Department of Mathematics, Nanjing University of Science and Technology, Nanjing, 210094, P. R. China (junma@njust.edu.cn).
Dong Wang (corresponding author) is with Department of Mathematics, University of Utah, Salt Lake City, UT, 84112 USA (dwang@math.utah.edu).
Xiao-Ping Wang is with Department of Mathematics, Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong (mawang@ust.hk).
Xiaoping Yang (corresponding author) is with Department of Mathematics, Nanjing University, Nanjing, 210093, P. R. China (xpyang@nju.edu.cn). snake model). The contour (explicitly represented by a parametric curve) evolves by the image gradient to the boundary of the desired object. Because one only tracks the explicit curve, the snake model is very efficient and requires low memory, allowing a fast evolution of an accurate boundary. However, the snake model suffers from numerical instabilities, and is difficult to automatically handle topological changes of the curve during the evolution (typically, it works only for a single closed curve) ( [16] , [9] ).
To improve the evolution of the contour of the snake model, Caselles et al. proposed geodesic active contours [8] in 1997. The contour is implicitly represented by a level set function that automatically allows topological changes including splitting and merging and simultaneous segmentation of single or multiple objects. The key idea of the level set method (LSM), introduced by Osher and Sethian [17] , is to represent a curve as the zero level set of a graph function defined in a higher dimensional space. Nowadays, the LSM has been widely used in many applications including computer vision, computational geometry, fluid dynamics, material science, and so on (see [18] and references therein for more details). In particular, using the LSM to implicitly represent the contour and approximately solve the active contour models becomes the most popular choice [19] .
Even the geodesic active contour model (GAC) allows to change the topology of the curve during evolution, it still suffers from numerical instability. In fact, even if the analytical model could generate a correct interface for all times, it might happen, for either the analytical or the numerical reason, that the gradient of the level set function would become "too small" or "too large" on the interface (i.e., on the zero-level set). The gradient of a level set function being too small will result in the location of the zero level set (the interface) being sensitive to perturbation. If the gradient is too large, one loses accuracy in the interface representation. To avoid this problem, the level set function is periodically reinitialized as a distance function from the interface, allowing to keep the norm of the gradient close to the unity and avoiding ill-conditioning. However, the reinitialization procedure usually involves many tricks, for example, it is hard to decide when it should be applied. Li et al. [9] proposed a penalty term to keep the regularity of the level set function during evolution. The core idea is to use the intrinsic property of the signed distance function: the magnitude of the gradient of the signed distance function equals one. A penalty term is then introduced to penalize this constraint. In addition, they proposed to keep the magnitude of the gradient of the level set function to be 1 in a neighborhood of the zero level set and the value of the level set function to arXiv:2007.00525v1 [cs.CV] 1 Jul 2020 be a constant at locations far away from the zero level set.
Compared to the classical snake model, GAC has achieved significant improvements by using the level set function to represent the contour. However, it is still inefficient because the level set is evolved by a time-dependent partial differential equation that has large computational burden to obtain a solution. Moreover, the solutions are usually obtained by explicit schemes that have stability issues regarding the choice of the time step.
Recently, Wang et al. ([20] , [21] ) proposed to use a characteristic function to represent the contour in region-based active contour models, where the perimeter of the contour is approximated by a heat kernel convolution with a characteristic function. Furthermore, they derived an iterative convolutionthresholding method (ICTM) to minimize a general energy functional with general fidelity terms, which enjoys the unconditionally energy-decaying property. Numerical experiments in ( [20] , [21] ) have shown that the ICTM is simple, efficient, and applicable to a wide range of region-based segmentation models. In [22] , Wang proposed a new ICTM-type method for reconstructing surfaces from point clouds, which can be interpreted as an anisotropic evolution of curves.
An alternative way to represent contours is to employ characteristic functions of regions enclosed by the contours (beyond level sets), but how to use a characteristic function to represent the contour in the geodesic active contour model is still an unsolved problem because it is difficult to represent the dynamic boundary curve evolution by a non-differentiable binary function. This paper aims to answer the following two questions:
• how to use a characteristic function to represent the geodesic active contour model? • how to develop a more efficient and stable minimization algorithm beyond the LSM? Motivated by ( [20] , [21] , [22] ), in this paper, we propose to use a characteristic function to implicitly represent the evolving contour, approximate the geodesic active contour energy functional and derive an efficient algorithm to minimize the energy. The main contributions are summarized as follows:
1. We propose to use characteristic functions rather than the widely used level set functions to represent the contours in the GAC energy functional. 2. We develop an iterative convolution-thresholding method to minimize the energy functional of GAC. Our algorithm is more efficient and stable than the classical gradient descent methods in the level set-based evolution, which is the most popular algorithm in this area. 3. We theoretically prove that our algorithm enjoys the energy-decaying property. 4. Experiments on synthetic, ultrasound, CT and MR images for the tumor and organ segmentation imply that our method is superior than traditional LSM in terms of efficiency (running time), which achieves dozens or hundreds of times acceleration. The paper is organized as follows. In Section III we introduce the characteristic function-based representation of the contours, give an approximation of the energy functional and derive an efficient algorithm for minimizing the energy. In Section IV we show the performance and verify the efficiency of our method by extensive experiments. In Section V we discuss the intuitive understanding on the ICTM and other potential applications. We draw some conclusions in Section VI.
For an image I on a domain Ω, there are usually three methods to represent the evolution of object boundary curves in it: 1) parametric curves, 2) level set functions, and 3) characteristic functions. Table I presents the main features of above three representation methods. In the following, we first review the parametric curve-and level set-based approaches.
1) Parametric curve representation: Let C(q) : [0, 1] → R 2 be a parametric curve. To find the object boundary, the classical snakes model [15] defines the following energy functional associated the curve C
where α, β, and λ are real positive constants. The first two terms belong to the internal energy that controls the smoothness of the contours, while the third term is the external energy that drives the contours towards the boundary of the object. The snakes model is a pioneer that formulates the image segmentation problem as an energy functional minimization problem. Using a parametrized planar curve to represent an object contour allows a fast evolution with a spline function method [16] . However, it suffers from the fixed topological property. For example, if there are more than one objects in a given image and the initial segmentation contour surrounds the objects, the snake model can not detect all objects. In other words, the classical snakes model can not directly deal with topological changes.
2) Level set representation: To address the drawback of the snakes model, Caselles et al. [8] proposed the wellknown GAC that is a geodesic computational problem in a Riemannian space whose metric is defined by the image information.
where g : [0, +∞) → R + is an edge indicator function. In general, g is defined by
where G σ is a Gaussian kernel with a standard deviation σ that is used to smooth the image. It is easy to see that g takes smaller values on the object boundary where the gradient magnitude is larger. Due to |C (q)|dq = ds, we write 
where |C| is the length of C and s is arc-length parameter that is intrinsic and related to the geometry of the curve. The energy functional can then be represented by embedding the dynamic contour C(s, t) as the zero level set of a time dependent level set function φ :
where δ is the Dirac delta function on the set φ = 0.
One of the most significant advantages of level set-based contour representation is that it can handle topological changes (e.g., merging and splitting) in a natural way, which is not allowed in parametric-based contour representation.
In practice, an area term is usually introduced to speed up the motion of the zero level set during the evolution, which is important when the initial contour is far away from the desired object boundaries. The energy functional is then defined by
where H is the Heaviside function. The corresponding L 2 gradient flow with respect to φ can be derived as
where δ is a smooth approximation of δ, which is defined by
and > 0 is a hyper-parameter that controls the band width of the non-zero region. Although the level set methods have the desired property on handling topology changes, their applications suffer from issues on numerical instability. To be specific, a level set function is usually defined as a signed distance function and typically develops irregularities due to numerical errors during evolution, which could destroy the stability of the level set evolution.
3) Distance regularized level set: To address the numerical instability problem, reinitialization was periodically used in earlier level set methods [11] to force the level set to be a signed distance function during the contour evolution. To avoid the above problem, Li et al. proposed the following distance regularized level set evolution (DRLSE) method
where µ > 0 is a weight hyper parameter and p(x) is a potential function that is used to keep the signed distance regularity of the level set function. Typically, the potential function p(x) can be defined as a single-well potential
or a double-well potential
In general, double-well potential is the default setting because it is more robust than single-well potential. With the smoothed Dirac delta function δ and Heaviside function H , one can derive the gradient flow of the energy functional (9) as
x .
Although the LSM has many advantages than the classical parametric curve-based methods and the distance regularization approach makes the LSM more stable, they still have a main drawback: low computational efficiency. This is because solving the level set-based GAC needs to update the level set function according to a partial differential equation (12) , which has large computational burden.
The characteristic function is defined by
where Γ is the object boundary and Ω Γ denotes the region inside Γ. It provides an alternative way to implicitly represent a curve, which not only owns the adaptively topological change property of the LSM but also is more computationally efficient and stable. In this paper, we propose to formulate the energy functional of GAC beyond existing level set methods. In particular, a characteristic function is introduced to approximate the energy (6), which allows us to design a more efficient algorithm for GAC. As shown in [23] , a general boundary integral can be approximated using the characteristic functions u by:
where τ is a free parameter, * represents convolution between two functions, and
).
Here n is the dimension of the Euclidean space and Γ could be an interface when n = 2 or a surface when n = 3. Esedoglu and Otto [24] established a novel framework on modelling and simulating the multiphase mean curvature flow with arbitrary surface tensions based on this approximation.
To keep the symmetry of the formulation with respect to u and 1 − u, combining with the area term, we propose to approximate the geodesic active contour energy functional (6) by
Thus, we arrive at the following energy functional minimization problem u * = arg min
where
and BV (Ω, R) denotes the space of functions with bounded variation.
B. Algorithm for problem (17) It is easy to see that the feasible set B of the energy functional minimization problem (17) is non-convex. To address this problem, we relax B to its convex hull K := {u ∈ BV (Ω, R)|u ∈ [0, 1]} and derive the following relaxed minimization problem
Furthermore, we prove the equivalence between (17) and (18) in the following lemma. Lemma 3.1: The original problem (17) is equivalent to the relaxed problem (18) , in the sense that if u * is a solution of (17), so is (18) . Vice versa.
Proof: On the one hand, letû = arg min
Then, it is obvious that
because B K.
On the other hand, letũ = arg min u∈K E τ (u), we can use reduction to absurdity to provẽ
Assume it is not true, then there exists a set A ⊆ Ω with nonzero measure and a > 0 such that the minimizer u * satisfies u * (x) ∈ (a, 1 − a), ∀x ∈ A.
Let u t = u * + tχ A where χ A is the characteristic function of A, we have u t ∈ K for any |t| < a. Directly computing the first and the second derivatives of
with respect to t, we have
Due to √ g ≥ 0 and √ g = 0 only on a set with zero
< 0 especially at t = 0 (i.e., u * ). However, this result contradicts with the assumption that u * is a minimizer ( d 2 E τ (u t ) dt 2 ≥ 0 at u * ). Thus, we haveũ ∈ B, and then
Finally, based on (19)- (21), we obtain
Next, we derive an iterative method to solve the relaxed problem (18) based on the sequential linear programming. At the k-th iteration u k , we compute the linearization (the first order Taylor expansion) of E τ (u) at u k as
where
Then, we can obtain the k + 1-th iteration u k+1 by solving the following linearized problem:
Following, we carry out the linearized minimization problem (25) in a pointwise manner. That is,
at each point x ∈ Ω. Due to the fact that the minimizer of a linear functional over a convex set must be reached at the boundary, we obtain the following solution
We summarize the proposed method in Algorithm 1 and term it as the iterative convolution-thresholding method (ICTM) because the main operators are only simple convolution and thresholding.
The iterative convolution-thresholding method (ICTM) for geodesic active contours Input: Image edge indicator function, τ > 0 and initialization u 0 ∈ B. Output: Segmentation results u * ∈ B;
while not converged do (1) Convolution. Fix u k , compute
For the convergence criteria, we stop the iteration if
where tol is a given error tolerance (1 × 10 −5 in this paper). Compared to the LSM, the proposed ICTM enjoys following advantages:
1. No requirement for the addition regularization term; The LSM in (9) needs an auxiliary regularization term to maintain the numerical stability, while our algorithm is intrinsic stable during iterations as demonstrated in the following section. 2. Fewer hyper-parameters; Specifically, the LSM has three model hyper-parameters 1 : α, λ, and µ (can be reduced to two by some normalizations) and two algorithm hyperparameters: time step and band width in δ and H , while our method reduces the number of hyperparameters to one model parameter λ and one joint model-algorithm hyper-parameter τ . 3. Less computational burden; The LSM needs to solve a partial differential equation for the evolution of the contour as shown in the gradient flow (12) , while our algorithm only alternates simple convolution and thresholding operations. 1 The standard deviation σ in the Gaussian kernel of edge indicator function (Eq. (2)) is excluded because this hyper-parameter is not always necessary. For example, if an input image is clean, we do not need to use Gaussian filter to smooth the image.
In this section, we prove that the Algorithm 1 is unconditionally stable for any τ > 0, which means the total energy E τ (u) is decreasing during the iteration. Thus, the proposed method can always converge to a stationary segmentation.
Theorem 3.1: Let u k (k = 0, 1, 2, · · · ) be the k-th iteration derived in Algorithm 1. We have
for any τ > 0.
Proof: As for E τ (u) defined in (16) , the linearization of E τ (u) at u k is defined by:
Direct calculation yields that
Because u k+1 is the solution from the sequential linear programming, we have L τ (u k+1 , u k ) ≤ L τ (u k , u k ). Then, we compute
Based on the semi-group property of the heat kernel convolution, i.e.,
we have
Therefore, we are led that 
In this section, we show the performance of Algorithm 1 on four different types of images, including synthetic images (2D), ultrasound (2D), CT (3D) and MR (3D) images, and compare with the well-known level set-based GAC [9] (DRLSE) to show the efficiency of the proposed method. The software code of LSM is obtained from Chunming Li's homepage 2 . For fair comparison, we also implement our algorithm with pure MATLAB code. All the experiments are ran on a Windows 10 laptop with Intel(R) Core(TM) i7-6700HQ CPU @ 2.60GHz, 24 GM RAM. Our code will be publicly available upon acceptance. In all experiment, if not specifically pointed out, we set σ in (3) to be 15 and in (8) to be 1.5.
We first apply the proposed method on two synthetic images to show that using characteristic functions can automatically handle the topological changes, such as merging and splitting, during evolution, which is the most important feature in level set-based methods. The hyper-parameters of the LSM follow the default setting in Chunming Li's code 2 . As for Algorithm 1, we set τ = 2 and choose the parameter λ as -0.3 and 0.3 for merging and splitting experiments, respectively. Figure 1 displays the segmentation results of the proposed ICTM and the LSM on two synthetic images. The first two rows show the merging process while the last two rows show the splitting process, using selected snapshots during the iteration. The corresponding number of iterations and running time are listed at the top of each image. Based on these segmentation results, we observe that 1. Both the ICTM and the LSM can adaptively handle the topological changes. Specifically, the contours can adaptively merge and split during iterations. 2. Both the ICTM and the LSM can achieve same (or similar) segmentation results on the two synthetic images because they are just different approximations to the same model (i.e., GAC). 3. Compared to the LSM, the ICTM achieves the same segmentation results with many fewer iterations and much less running time. In particular, ICTM is approximate 58 times and 27 times faster than the LSM on the contour merging and splitting results, respectively.
To validate the performance and efficiency of our method on real images, we apply it to breast nodule segmentation in ultrasound images. In this experiment, we use the BUSI dataset [25] and select 100 random benign breast ultrasound images to compare the efficiency between the LSM and the ICTM. The same rectangle initialization is generated based on ground truth for images. The width and height of each rectangle is half of the major axis and minor axis of the ground truth's bounding box, respectively. For a fair comparison, we apply grid search to each method to find the best set of hyperparameters that can achieve the best average Dice similarity coefficient (DSC) on this dataset. The main contribution in this work is a faster algorithm. Thus we try our best to achieve the best performance for both LSM and ICTM, and then compare their running time. Specifically, we search α, λ and the time step for the LSM in {1, 2, · · · , 10}, {−1, −2, · · · , −5}, and {1, 2, · · · , 5}, respectively. For ICTM, we search τ and λ in {1, 2, · · · , 4} and {−0.1, −0.2, · · · , −0.4}, respectively. Table II shows the quantitative results of breast ultrasound image segmentation. The proposed ICTM method is slightly better than the LSM with 1.35% improvements in average DSC. For the average running time, the LSM takes about 74 seconds per image while the proposed ICTM only needs 0.66 seconds per image, achieving more than 100 times acceleration. Figure 2 displays 10 random selected images with segmentation results by the LSM and our ICTM, with corresponding DSC (yellow) and running time (white) printed on. The two methods start with the same initialization, and we observe that • The ICTM achieves similar or even better segmentation results than the LSM, indicating the feasibility of applying the ICTM into real image segmentation. • The ICTM requires much less running time than the LSM, implying the high efficiency.
To validate the effectiveness of the ICTM on 3D organ segmentation tasks, we apply it to lung segmentation in COVID-19 CT scans. We use the public COVID-19-CT-Seg dataset ( [26] , [27] ) and select 10 random lung CT scans of earlier COVID-19 patients for segmentation experiments. The number of image slices ranges from 40 to 400. We set the same initialization for both the ICTM and the LSM, which is generated by eroding the ground truth with a sphere structure element (radius=10).
For fair comparison, we apply the similar grid search (as shown in Section IV-B) to tune the hyper-parameters. Quantitative and qualitative segmentation results are displayed in Table II and Figure 3 , respectively. We observe that both two methods achieve comparable average DSC without significant differences, which implies that the ICTM can be an alternative choice (beyond the LSM) for the GAC. As for the efficiency, 
To validate the effectiveness of our ICTM on 3D lesion segmentation tasks, we apply it to liver lesion segmentation in liver MR scans. We randomly collect 20 liver MR scans from a local hospital. Three experienced radiologists manually annotate them, and majority vote is used to generate final labels. The image sizes range from 256 × 256 × 105 to 400×400×120. The initialization of each MR scan is a cuboid inside the tumor that is generated from the ground truth. Figure  4 (a) and (d) show some initialization results (red rectangle) in 2D slices. It should be noted that not all tumor slices have initializations such as the first and third images in Figure 4 (a) because the segmentation method is applied in a 3D manner. In fact, only half of the tumor slices have initializations.
For fair comparison, we also apply the same grid search (as shown in Section IV-B) to tune the hyper-parameters. Quantitative and qualitative segmentation results are displayed in Table II and Figure 4 , respectively. We observe that, again, both two methods achieve similar average DSC without significant differences, indicating that the ICTM can obtain similar results to the LSM for the GAC. However, the ICTM requires much less running time compared to the LSM, which is about 68.8 times acceleration, implying the high efficiency of the ICTM.
In Section IV, we have applied the proposed ICTM on synthetic, ultrasound, CT and MR images to show its effectiveness on nodule, organ and lesion segmentation. Compared to the LSM, the ICTM obtains similar or even better results but achieves dozens or hundreds of times acceleration. Furthermore, we discuss the intuitive understanding on the advantages of the ICTM (e.g., efficiency) and many potential applications especially in the modern deep learning era.
This is contributed to the simple and inherent features (Table I) of using characteristic functions to implicitly represent a contour. Specifically, there are several main reasons:
1. Each iteration in the ICTM (Algorithm 1) is much simpler than each iteration in the LSM as shown in (12). 2. The ICTM directly minimizes the geodesic active contour energy functional, while the LSM usually needs to minimize additional energy term to stabilize the iteration. 3. At each iteration, our ICTM can find the optimal minimizer of the linearized functional. This is because the optimal minimizer of a linear functional over a convex set can be reached at the boundary. Moreover, the minimizer can give a smaller value in E τ because the graph of the functional E τ (concave) is always below its linear approximation. This accelerates the convergence of the ICTM. In the LSM, one needs to solve the level setbased partial differential equation (PDE) with a relatively small time step. This step more or less restricts the decay of the energy (at least not optimal). What's worse, the reinitialization step after it (or adding penalty terms in the level set equation) usually increases the energy, which decreases the value of the energy minimized at each iteration (increases the value of E τ ). This makes the LSM converge slower. In the ICTM, thanks to the concavity of E τ , the minimizer at each iteration automatically gives a new partition (i.e., the minimizer automatically remains at characteristic functions). No reinitialization and related regularization techniques are needed in the ICTM. We plot the energy curves of both methods during iterations for the synthetic image segmentation (see Figure 1 ), and observe that ICTM converges faster than the LSM. In particular, we observe more oscillations in the energy curve for the LSM.
The ICTM is just a method for approximately solving an image segmentation model (i.e., minimizing an objective functional). On one hand, from extensive numerical experiments in this paper, we claim that, our our ICTM is much more efficient than the LSM for various segmentation tasks. On the other hand, the proposed ICTM may not be outside of the application scope of the LSM. This is because both two methods aim to approximately solve the same GAC (6), which mainly determines the accuracy of the segmentation. For better normalization, we apply max-min normalization to the energy values. We zoom in the energy curve of the last 50 iterations to better show the difference between two methods. It can be found that ICTM converges faster than the LSM. In particular, more oscillations occur in the energy curve for the LSM.
C. What is the role of the traditional model-based geodesic active contours in modern deep learning era?
Although deep learning-based segmentation methods have been increasingly popular and dominating current segmentation tasks, these methods generally require much annotated training data that is difficult to obtain in medical images. Deep learning techniques are still open for many mathematical explanations and theories, which may cause incomprehensible segmentation results.
Model-based GAC has interpretable nature and still plays important roles in following three circumstances:
1. Assisting radiologists to annotate medical images (e.g.; tumor segmentation [28] ). 2. Serving as a post-processing method to refine the segmentation results that are generated by deep learning-based approaches (e.g.; breast tumor [29] , liver [30] , and dental root [31] segmentation). 3. Explicitly embedding shape information (e.g.; a left ventricle shape model that is learned by auto-encoder network can be embedded into the GAC [32] ). 4. Reformulating the active contour model as a loss function to guide CNNs to learning richer features, such as Mumford-Shah loss [33] , level set loss [34] , and active contour loss [35] . The proposed ICTM is expected to be applied to some of these situations directly, obtaining a dramatic acceleration. We leave these applications as our future work and will be reported elsewhere.
In this paper, we have proposed an efficient iterative convolution-thresholding method (ICTM) to solve the wildly used geodesic active contours (GAC). The method mainly relies on a characteristic function-based representation for the contour and an integral approximation of the energy functional. A relaxation and linearization approach is used to derive the ICTM method. Extensive numerical experiments on four different types of images are presented to show the performance of the proposed method, indicating a dramatic improvement in the efficiency (compared to the level set-based approaches).