Active Arcs and Contours

—The level set method [31] is a commonly used framework for image segmentation algorithms. For edge detection and segmentation models, the standard level set method provides a ﬂexible curve representation and implementation. However, one drawback has been in the types of curves that can be represented. In the classical level set method, the curve must enclose an open set ( i.e. loops or contours without endpoints). Thus the classical framework is limited to locating edge sets without endpoints. Using the curve representation from [34], [35], we propose a segmentation and edge detection method which can locate arcs ( i.e. curves with free endpoints) as well as standard contours. Within this new framework, we propose a variational segmentation model to detect general edge structures. The proposed energy is composed of two terms, an edge set regularizer and an edge attractor. The proposed energy is related to the Mumford and Shah model [27] for joint segmentation and restoration in terms of an asymptotic limit. Numerical results are given on images with a variety of edge structures.


I. INTRODUCTION
Edge detection and image segmentation remains a largely open problem in image processing. The focus of edge detection is to locate important edges between objects within an image or to locate features which have sharp changes in intensity. One point of difficulty is in the large class of possible edge structures that can occur. Much of the literature defines the edge structure to be composed of curves without endpoints, but in this work, we apply the curve definitions from [35] to extend some classic models to a larger class of edges. In particular, we propose using a level set method for curves that may include free endpoints to locate edges. In this work we will refer to curves that have endpoints as either arcs, curves with free endpoints, or free curves.
The original Active Contours model (Snakes) [16] proposed an energy minimization method for the location of the edge set. Let f be the given (possibly degraded) image defined on a domain Ω ⊂ R 2 and let Γ(s) : [0, 1] → R 2 be a parameterized curve representing the edge set, then the model is (essentially) as follows: where λ 1 , λ 2 > 0 are parameters. The first integral is defined by two terms related to the first and second derivative of the curve and is a regularizer, defining a particular smoothness in the edge structure. The second integral attracts the curve to places where the gradient of the image is large. In this way, the Snakes model is an energy based edge detector -looking for a smooth edge set along high values of |∇f |. A general gradient based edge detector is defined as a nonnegative function of the gradient of the image, g(|∇f |), such that g(x) → 0 as x → ∞. The Geodesic Active Contours (GAC) model defines an energy using the function g as follows [5], [17]: The model is formulated using the level set method, proposed by Osher and Sethian [31]. The level set method provides an implicit representation for curves by defining them as the zero level set of a Lipschitz continuous function φ : Ω → R. Using the level set framework allows the curve to undergo changes in topology and allows for the formation of corners. The level set method restricts the class of possible edge sets to curves made up of segments without endpoints or that terminate at the boundary of the domain. Equation (2) is related to the problem of finding geodesics in a Riemannian space, defined by a metric g. Therefore, minimizers of equation (2) will have minimal lengths with respect to the edge-based metric and reside along the points where g (|∇f |) is small, i.e. where |∇f | is large. In the works [23] and [39], Finsler spaces were used in place of the Riemannian spaces, thus extending GAC to Finsler Active Contours.
Alternatively, Mumford and Shah based methods jointly segment and restore images. The original Mumford and Shah model (MS) [27] finds a piecewise smooth image u and an edge set Γ by minimizing the following functional: The first term is a regularizer on the reconstructed image u, while the second term encourages u to remain close to the initial data f . The third term is a regularizer on the edge set Γ (the Hausdorff measure or length). Minimizing curves Γ of equation (3) are known to be composed of C 1,1 segments whose endpoints either terminate perpendicularly to the boundary of the domain, at a triple junction, or at points that do not connect back to the curve (arcs). For more theoretical results on minimizers, see for example [11], [14], [25], [26].
One approximation to the MS model, based on Γconvergence, is the Ambrosio and Tortorelli model (AT) [2]. To represent the edge set, the model uses a function v ∈ [0, 1], which acts as a smooth indicator function on the edge set. The Ambrosio and Tortorelli functional is as follows: where > 0 and η = o( ). As → 0, equation (4) Γconverges to the weak formulation of equation (3). The third integral in equation (4) is an elliptic approximation to the length term from equation (3). This approximation has the advantage of capturing all types of edges theoretically possible in the MS model. The disadvantage of the elliptic approximation to the edge set is that the width of the edges are a function of the parameter , which makes the edge set thicker than in the level set based methods.
Using the level set method, Chan and Vese proposed a piecewise constant version of the MS model [7], where the reconstructed image is equal to c 1 inside the region enclosed by the curve and c 2 outside the region enclosed by the curve. The Chan and Vese functional (CV) is as follows: where the curve Γ is approximated using one level set function. The model fits a given image to two constants, while also minimizing the size of the boundary between the two regions defined by the constants. This model was also extended to piecewise constant vector valued images in [6], to 4 regions and to piecewise smooth solutions in [38] and [37], and to multilayer piecewise constant images in [8].
The level set method for image segmentation has seen much success in the literature, with a wide range of applications. However, classical level set method is unable to capture edge structures with free endpoints. Using the work of [34] (which is based on [35]) for image segmentation using arcs, we propose a level set based variational edge detector derived from the Mumford and Shah model. The model is able to capture all the edge structures from the MS model, but has the advantage of doing so quickly. For more on arcs and general curve structures see [9], [18]- [21], [24].
The paper is organized as follows. Section II describes the proposed model and its connection to existing models. Theoretical and analytical remarks are provided in section III. Lastly section IV and section V present the numerical method and experimental results, respectively.

A. Level Set Method with Free Endpoints
From [34], [35], we can represent a curve Γ as the intersection of a zero level set of a Lipschitz continuous function φ : Ω → R with the positive part of a second level set ψ : Ω → R. The second level set function acts as a partitioning function on the first. With this definition for Γ, Γ := {x ∈ Ω φ(x) = 0, ψ(x) > 0}, the curve is now able to have loops, segments terminating at the boundary of the domain, and segments with free endpoints.
Applying this definition for the edge set to the Length functional, a term which appears in many segmentation models, yields the following [12]: where H is the Heaviside function, which is 1 for positive arguments and 0 otherwise. Taking a smooth approximation H for the Heaviside function, equation (6) becomes: where δ = H is a smooth approximation of the Dirac delta function. To minimize equation (7), the first variation is embedded in the following dynamic scheme: The partial differential equation for φ defines a mean curvature flow for the zero level set of φ in the region defined by the positive part of ψ, while the differential equation for ψ defines a Hamilton-Jacobi equation for the intersection of the zero level sets of φ and ψ.

B. Description of the Model
In [27], the authors showed that the limit as µ → ∞ (with ν inversely related to µ) of the Mumford and Shah model (in the sense described in Section III) is as follows: where ν ∞ > 0 is a constant. The first term penalizes the length of the curve, while the second term encourages the curve to align with places of high gradient (which corresponds to the edges). Generalizing the functional above, we propose the following energy: where g is a strictly positive function which can be a general feature detector only based on the given image and its derivatives (see [5], [17]) and Γ is a Lipschitz curve. Next, using the definition from Section 2, Γ := {x ∈ Ω φ(x) = 0, ψ(x) > 0}, and assuming the Dirac delta and Heaviside functions are replaced by their smooth versions (dropping the subscripts), we have the following level set version of equation (8): Minimizing the equation above may not result in a wellbehaved curve; thus, to ensure sufficient continuity on the curve, a Lipschitz regularizer is used: where β 1 , β 2 ≥ 0. In practice, β 2 can be close to or (in some cases) equal to zero and 0 < β 1 << ||g|| L ∞ is sufficient. Equation (9) is non-convex, so to find minimizers we find the Euler-Lagrange equations and embed them in a time dependent system (the gradient descent method). The resulting system of equations is: is the renormalized infinity Laplacian (see Appendix for derivation of the non-regularized system). The infinity Laplacian smoothes the level set functions in the normal direction to each of its level sets, keeping the functions Lipschitz. For theoretical results, see [10], [13], [22], and, for a convergent discretization, see [30].

C. Relation to Active Contours
The proposed model is similar in nature to both the Geodesic Active Contours model [5] and the Active Contours (Snakes) model [16]. In the general form, the Snakes energy is composed of two main terms -a curve regularizer (elasticity) and an edge attractor (based on the L 2 norm of the gradient of the original image). In the proposed model, the Lipshitz regularizer acts as the curve regularizer, the L 2 norm of the normal derivative of the original image is an edge attractor, and the geodesic length behaves as a length regularizer and an edge attractor.
The energy proposed in this paper can be seen as an extension to the GAC energy, with the addition of the alignment term. One crucial difference is the GAC energy would not be able to detect edges that contain arcs. In the case of arcs, the second term in our energy prevents the curve from shrinking along the tangents of the edge set, allowing for non-trivial minimizers.

D. Relation to Mumford and Shah Model
Recall the full Mumford and Shah Energy below.
As in the piecewise constant case, the proposed model can be seen as a special case with a particular choice of u. Takeũ to agree with the intial data f outside of a small region of the edge set Γ. Let (r, s) be the curvilinear coordinates defined by the normals and tangents of Γ and takeũ to have the following expression in the region near the edge set. u(r, s) := f (r, s) + sign(r)e − |r| ∂f ∂r (0, s) Choosing µ := 1 and ν = 2 ν ∞ yields the following expression [27]: From here we can see that the E M S,∞ (Γ) energy is a natural limit to the MS model. In a similiar way, the proposed functional E(Γ) is a limit of the MS model with geodesic length.

III. THEORETICAL REMARKS
In this section, several theorems and remarks present wellposedness like conditions of the model.

A. Consistency of the Level Set Formulation
In terms of the level set formulation of the energy (equation ( 9)), the approximation of the length term is consistent, in the limit.
and let φ and ψ be Lipschitz continuous, then For the proof of the theorem above see [34]. The theorem above justifies the differential substitution ds = |∇H (φ)|H (ψ)dx, which is used to transform equation (8) to the first term in equation (9).

B. Limiting Relation to MS Functional
As mentioned in Section II D, the proposed model can be seen as a special case of the MS model. Let ν ∞ = µν 2 ; then the following theorem provides the relationship between the Mumford and Shah functional and equation (8) when g = ν ∞ [27].
Theorem III.2. If f is Lipschitz, Γ is the union of finitely many Lipschitz curves, and neither Γ nor the ∂Ω contain cusps, then as µ → ∞ the following holds: This shows for sufficiently smooth initial conditions, that this edge based model is more closely related to the MS function than other edge detectors. Also, the theorem gives one possible method for solving the full MS energy. First locating the edges using the proposed method, then smoothing the initial image away from the edge set.

C. Well-posedness for fixed endpoints
In the case where the endpoints are fixed, the unregularized model provides a well behaved minimizing curve. We can write the energy in the form of an integral of a Finsler metric along the curve Γ [27]: where dx and dy are the tangent plane coordinates and Since for any constant λ, F (x, y, λdx, λdy) = |λ|F (x, y, dx, dy), it is indeed a Finsler metric.
Minimizing the proposed energy (without the Lipschitz regularization) with fixed endpoints is a well-posed problem. Assuming sufficient regularity on Γ and f , the existence and uniquess of minimizers can be deduced from the Hopf-Rinow theorem [3]. This provides some intuition on the proposed model, in that the ill-posedness comes from the endpoint behavior or the curve continuity. For example, in the case of curves without endpoints, our model attains well behaved minimizers for sufficiently large c > 0, where g ≥ c.

D. Behavior of Minimizing Curves
When the endpoints are not fixed, equation (8) is not generally well-posed. In this unregularized model, two issues may arise. First, if ||∇f || 2 ∞ ≤ min x g then the minimizer is the empty set. However, if there exists a pointx such that |∇f (x)| 2 > g(x), then the minimizing curve will be composed of an infinite number of curves in a small neighborhood ofx with inf E(Γ) = −∞. Thus, minimizing the original energy for a general curve Γ is ill-posed. In fact, in [27] the authors conjectured that the problem is actually well-posed in the sense of currents and not curves.
With the proposed regularization, intuitively both cases can be handled. In the case when there exists a point such that the integrand is negative, the minimizing curve is made up of infinitely many nearly parallel curves. However, the Lipschitz regularization prevents the creation of extraneous curves since, as the number of components of the minimizing curve grows, the Lipschitz norm grows also. In the case when the global minimizer is the empty set (which is also true of the GAC model), there is both practical and theoretical support for the existence of an approriate edge set.
The non-convexity of the energy allows for local minimizers and thus non-trivial solutions. This is true for curves without endpoints, but in the case of arcs the local minimizer can shrink tangentially (to the edge) until it becomes the empty set. However, in practice, if the second level set function ψ has only partially converged, the method yields approriate edge sets. This is essentially equivalent to detecting an approximate location of the endpoints and fixing them. This tranforms the problem into locating the endpoints with a few iterations and then fixing the endpoints and minimizing a well-posed problem (see Section III.C).

IV. NUMERICAL METHOD
The proposed model is non-convex, and the resulting system of partial differential equations is non-linear. To ensure a smooth flow of the system and better posed equations, we implement a Sobolev gradient descent method. This is done by preconditioning the Euler-Lagrange equations by the smoothing operator (I − ∆) −1 , and then embedding the new equation in a dynamic scheme. In particular, for this problem we have the following preconditioned equations: where the function G = [G 1 , G 2 ] is defined as with Neumann boundary conditions for both φ and ψ. In practice, with the addition of the Sobolev gradient, the regularizer on ψ is not necessary, so one may take β 2 = 0. For more on Sobolev gradients, see [4], [15], [28], [32], [33], [36]. The Euler-Lagrange equations (defined as the function G) are discretized using forward differences for the gradient and backwards differences for the divergence. Also, the discretization is explicit in time. The terms that appear in the denominator are replaced with a regularized version to avoid dividing by zero.
To iterate forward in time, we introduce the auxiliary variable v = [v 1 , v 2 ] and solve v j − ∆v j = G j at every time step using the Gauss-Seidel sweeping method (which converges within a small number of iterations). Then using v, the level sets are marched forward in time by a forward Euler step: φ n+1 = φ n + ∆t v 1 and ψ n+1 = ψ n + ∆t v 2 .
Although the Euler-Lagrange equations are completely explicit, the Sobolev gradient inversion is solved using a semiimplict method, which adds stability. The table below outlines the complete algorithm.
The method can be sensitive to the intialization, thus by applying a simple edge detector one can locate a region of interest. To measure convergence, the difference in the energy between iterations or the difference in the curves' location can be used, both yielding similar results. In practice, this alternating minimization scheme needs only a few iterations to converge.

Algorithm
Initialize φ 0 , ψ 0 while Not Converged do Step 1: Compute G 1 (φ n , ψ n ), Solve for v n 1 Step 2: Iterate Forward to φ n+1 Step 3: Compute G 2 (φ n+1 , ψ n ), Solve for v n 2 Step 4: Iterate Forward to ψ n+1 end while V. EXPERIMENTAL RESULTS In this section, we present several experimental results on both synthetic and real images. For the smooth versions of the Heaviside and Dirac delta functions, we use H (x) = 1 2 + 1 π arctan x and δ (x) = 1 π 2 +x 2 , with = 1. In Sections V A-C, g (f ) is constant (i.e. g = C). The time step ∆t is taken to be between 0.001 and 0.05. The edge set can be composed of combinations of curves with and without endpoints. To test these cases, two noisy synthetic images are used. Figure 1 includes two arc edge sets and is initialized by using one arc. Within two iterations the curve locates the edge structures. The last few iterations refine the curve and remove non-edge components.

A. Examples
In Figure 2, the edge set is composed of an arc edge and an edge without endpoints. The curve is initialized locally to the edge structure and converges within a few iterations. Both structures are captured using the same level set functions.

B. Edges without Endpoints
In the case where the edge set is a priori known to not include arc edges, the energy (equation (9)) reduces to the following equation which is equivalent to equation (9), with ψ being a fixed positive function. In this case, the proposed model gives different results than existing models. In Figure 3, using multiple curves (same level set function) as the initial condition, the Shapes image is segmented. The flow is qualitatively different than those of CV [7], GAC [5], or Snakes [16] due to both the non-Riemannian metric and the Lipschitz regularizer. These allow the curve to form kinks and corners. Notice that in the intermediary iterations, the curve grows along the edges (see Figure 3 (c)). In Figure 4, the method is applied to an illusory contour which does not exist in the sense of sharp intensity changes. The curve is initialized around the contour and shrinks until capturing the shape.

C. Edge Sets with Arcs
Edge sets with arcs appear frequently in astronomical and physical imaging. We apply our method to locating arc edges in two real images. In Figure 5, the front of a comet is located using an arc as the initial curve. Within a couple of iterations, the front tip is located, while several more are needed to capture the entire comet front.
In Figure 6, the method is applied to a plasma shock front. The edge is of low constrast and has a long tail of decay; however, our method is able to accurately capture the full edge set.
In Figure 7, the results in Figure 5 are compared to solutions from existing methods. The CV and GAC method both locate a loop around the interior of the comet and miss both the comet tip and the low constrast comet sides. On the other hand, the AT [2] edge function and Canny edge detector locate the tip correctly and approximate the edge set by an arc. However, both of these methods locate the high constrast interior region of the comet and thus over-segment the image. In Figure 8, the solution from Figure 6 is compared to the same existing methods. The CV method locates the lighter region in the top right corner of the images and it does not locate the plasma shock. The GAC solution locates the region of high constrast along one of the sides of the shock; however, it does not locate the other side of the shock (the top part) and undersegments the front. Once again, the AT edge function and Canny edge detector both over-segment the solution and do not locate the low constrast regions.
The comparisons show the improvements and accuracy gained by our proposed model.

D. Geodesic Length
In this section, we take g (f ) = C(f +δ), where δ is a small postive constant and C is a positive constant. The parameter δ is needed to keep g strictly positive and thereby acts as a small length regularizer.
In Figure 9, the method is applied to region detection (in particular linear-structure detection). The result (Figure 9 (d)) is qualitatively different from that of the GAC solution ( Figure  9(e)). The GAC is smooth along the entire curve, missing some of the local geometry that our solution captures. Figure 10 displays our method applied to linear-structure detection with multiple intializations. In both cases the final curve locates the bridge structure correctly. Notice that as long as the curve is initialized locally to the desired structure, the different initializations converge to similar final edge sets (see Figure 10 (d) and (h)).

VI. CONCLUSION
This work presents a new image segmentation and edge detection technique using a free curve representation of the edge set. Within this free curve framework, we propose a

ACKNOWLEDGMENT
This research was made possible with support by the Department of Defense (DoD) through the National Defense Science and Engineering Graduate Fellowship (NDSEG) Program and by NSF DMS-0714945. The author would also like to thank Professor Luminita Vese for her helpful discussions.

APPENDIX A DERIVATION OF EULER-LAGRANGE EQUATIONS
Let φ, ψ ∈ W 1,∞ and assume that the absolute value, Dirac delta, and Heaviside functions are replaced by their differentiable approximations and recall the unregularized energy below.