Global Minimum for a Finsler Elastica Minimal Path Approach
Abstract
In this paper, we propose a novel curvature penalized minimal path model via an orientationlifted Finsler metric and the Euler elastica curve. The original minimal path model computes the globally minimal geodesic by solving an Eikonal partial differential equation (PDE). Essentially, this firstorder model is unable to penalize curvature which is related to the path rigidity property in the classical active contour models. To solve this problem, we present an Eikonal PDEbased Finsler elastica minimal path approach to address the curvaturepenalized geodesic energy minimization problem. We were successful at adding the curvature penalization to the classical geodesic energy (Caselles et al, 1997; Cohen and Kimmel, 1997). The basic idea of this work is to interpret the Euler elastica bending energy via a novel Finsler elastica metric that embeds a curvature penalty. This metric is nonRiemannian, anisotropic and asymmetric, and is defined over an orientationlifted space by adding to the image domain the orientation as an extra space dimension. Based on this orientation lifting, the proposed minimal path model can benefit from both the curvature and orientation of the paths. Thanks to the fast marching method, the global minimum of the curvaturepenalized geodesic energy can be computed efficiently.
We introduce two anisotropic image datadriven speed functions that are computed by steerable filters. Based on these orientationdependent speed functions, we can apply the proposed Finsler elastica minimal path model to the applications of closed contour detection, perceptual grouping and tubular structure extraction. Numerical experiments on both synthetic and real images show that these applications of the proposed model indeed obtain promising results.
Keywords:
Minimal Path Eikonal Equation Curvature Penalty Anisotropic Fast Marching Method Image Segmentation Tubular Structure Extraction∎
1 Introduction
Snakes or active contour models have been studied considerably, and used for object segmentation and feature extraction for almost three decades, since the pioneering work of the snakes model proposed by Kass et al (1988). A snake is a parametrized curve (locally) that minimizes the energy:
with appropriate boundary conditions at the endpoints and . and are the first and secondorder derivatives of the curve respectively. The positive constants and relate to the elasticity and rigidity of the curve and, hence, weight its internal force. This approach models contours as curves locally minimizing an objective energy functional that consists of an internal and an external force. The internal force terms depend on the first and secondorder derivatives of the curves (snakes), and, respectively, account for a prior of small length and of low curvature of the contours. The external force is derived from a potential function , which depends on image features such as the gradient magnitude, and it is designed to attract the active contours or snakes to the image features of interest such as object boundaries.
The drawbacks of the classical active contours or snakes model (Kass et al, 1988) are its sensitivity to initialization, the difficulty of handling topological changes, and the difficulty of minimizing the strongly nonconvex path energy as discussed by Cohen and Kimmel (1997). Regarding initialization, the active contours model requires an initial guess that is close to the desired image features, and preferably enclosing them because energy minimization tends to shorten the snakes. The introduction of an expanding balloon force allows the model to be less demanding on the initial curve given inside the objective region (Cohen, 1991). The issue of topology changes leads, on the other hand, to the development of active contour methods, which represent object boundaries as the zero level set of the solution to a PDE (Osher and Sethian, 1988; Caselles et al, 1993; Malladi et al, 1994; Caselles et al, 1997; Yezzi et al, 1997).
The difficulty of minimizing the nonconvex snakes energy (Kass et al, 1988) leads to important practical problems because the curve optimization procedure often becomes stuck at a local minimum of the energy functional, which makes the results sensitive to curve initialization and image noise. This limitation is still the case for the level set approach on geodesic active contours (Malladi et al, 1995; Caselles et al, 1997). To address this local minimum sensitivity issue, Cohen and Kimmel (1997) proposed an Eikonal PDEbased minimal path model to find the global minimum of the geodesic active contours energy (Caselles et al, 1997), in which the penalty associated to the secondorder derivative of the curve was removed from the snakes energy. Thanks to this simplification, a fast, reliable and globally optimal numerical method allows to find the energy minimizing curve with prescribed endpoints, namely, the fast marching method (Sethian, 1999), which is based on the formalism of viscosity solutions to Eikonal PDE. These mathematical and algorithmic guarantees of Cohen and Kimmel’s minimal path model (Cohen and Kimmel, 1997) have important practical consequences, leading to various approaches for image analysis and medical imaging (Peyré et al, 2010; Cohen, 2001).
In the basic formulations of the minimal pathsbased interactive image segmentation models (Appleton and Talbot, 2005; Appia and Yezzi, 2011; Mille et al, 2014), the common proposal is that the object boundaries can be delineated by a set of minimal paths constrained by userprovided points. In (Li and Yezzi, 2007; Benmansour and Cohen, 2011), vessels were extracted under the form of minimal paths over the radiuslifted space. Therefore, each extracted minimal path consists of both the centerline positions and the corresponding thickness values of a vessel.
To reduce the user intervention, Benmansour and Cohen (2009) proposed a growing minimal path model for object segmentation with unknown endpoints. This model can recursively detect keypoints, each of which can be used as a new source point for the fast marching front propagation. Thus this model requires only one userprovided point to start the keypoints detection procedure. Kaul et al (2012) applied the growing minimal path model to detect complicated curves with arbitrary topologies and developed criteria to stop the growth of the minimal paths. Rouchdy and Cohen (2013) proposed a geodesic voting model for vessel tree extraction by a voting score map that is constructed from a set of geodesics with a common source point.
Recently, improvements of the minimal path model have been devoted to extend the isotropic Riemannian metrics to the more general anisotropic Riemannian metrics by taking into account the orientation of the curves (Bougleux et al, 2008; Jbabdi et al, 2008; Benmansour and Cohen, 2011). Such orientation enhancement can solve some shortcuts problems suffered by the isotropic minimal path models (Cohen and Kimmel, 1997; Li and Yezzi, 2007). Kimmel and Sethian (2001) designed an orientationlifted Riemannian metric for the application of path planning, providing an alternative way to take advantage of the orientation information. This isotropic orientationlifted Riemannian metric (Kimmel and Sethian, 2001) was built over a higher dimensional domain by adding an extra orientation space to the image domain.
Bekkers et al (2015) considered a datadriven extension of the subRiemannian metric on SE(2), which shows that the subRiemannian structure outperforms the isotropic Riemannian structures on SE(2) in retinal vessel tracking. The numerical solver used by Bekkers et al (2015) is based on a PDE approach with an upwind discretization and iterative evolution scheme, which requires expensive computation time. To solve this problem, Sanguinetti et al (2015) used the fast marching method (Mirebeau, 2014a) as the Eikonal solver to track the subRiemannian geodesics. The subRiemannian geodesic model (Petitot, 2003) reintroduced curvature penalization to the geodesic energy, similar to the Euler elastica bending energy (Nitzberg and Mumford, 1990) considered in this paper, yet it differs in two ways: firstly, the Euler elastica bending energy involves the squared curvature, a stronger penalization than Petitot’s subRiemannian geodesic energy which is roughly linear in the curvature. Secondly, minimal geodesics for Petitot’s subRiemannian model occasionally feature cusps, which sometimes may not be desirable for the applications of interest. In contrast, Euler elastica curves can avoid such cusps.
Schoenemann et al (2012) proposed a model to address the problems of using curvature regularization for regionbased image segmentation by a graphbased combinational optimization method. This curvature regularization model can find a solution that corresponds to the globally optimal segmentation which is initializationindependent, and has proven to obtain promising segmentation and inpainting results especially for objects with long and thin structures. Ulen et al (2015) proposed a curvature and torsion regularized shortest path model for tubular structure segmentation, where the curvature and the torsion were approximately computed by BSplines. The solution of their energy functional including curvature and torsion penalization terms can be efficiently obtained by using line graphs and Dijkstra’s algorithm (Dijkstra, 1959).
Tai et al (2011) presented an efficient method to solve the minimization problem of the Euler elastica energy, and demonstrated that this fast method can be applied to image denosing, inpainting, and zooming. Other approaches of interest using the curvature penalization include the image segmentation models such as (ElZehiry and Grady, 2010; Schoenemann et al, 2011; Zhu et al, 2013) and the image inpainting model introduced by Shen et al (2003).
1.1 Motivation
In contrast with the classical snakes energy (Kass et al, 1988), Eikonal PDEbased minimal path methods are firstorder models, which do not penalize the secondorder derivative of a curve, i.e., the curvature, and therefore do not enforce the smoothness of the geodesic, leading sometimes to undesired result as shown in Fig. 1, in which we would like to extract a boundary that is as smooth as possible between the two given points indicated by red and green dots. In Fig. 1a we show the edge saliency map. Figs. 1b and 1c are the minimal paths obtained by using the isotropic Riemannian metric (Cohen and Kimmel, 1997) and the anisotropic Riemannian metric (Bougleux et al, 2008), respectively, both of which are unable to find expected smooth boundaries and suffer from the shortcut problem due to the lack of curvature penalization in these metrics. In contrast, the minimal path model presented in this paper reintroduces the curvature, in the form of weighted Euler elastica curves as studied in (Nitzberg and Mumford, 1990; Mumford, 1994). Therefore, the geodesics extracted by the proposed metric can catch the smooth object boundary, as shown in Fig. 1d, with arrows indicating the corresponding tangents at the given positions.
1.2 Contributions
The contribution of this paper is three fold:

Firstly, we propose a novel globally minimized minimal path model, namely, the Finsler elastica minimal path model, with curvature penalization and Finsler metric. We establish the connection between the Euler elastica curves and the minimal paths with respect to a Finsler elastica metric. With an adequate numerical implementation, leveraging orientation lifting, asymmetric Finsler metrics and anisotropic fast marching method, the proposed model still allows to find the globally minimizing curves with prescribed points and tangents.

As a second contribution, we present the mathematical convergence analysis of the proposed Finsler elastica metrics, and of the associated Finsler elastica minimal paths. Furthermore, we discuss numerical options for geodesic distance and minimal paths computation, and settle for an adaptation of the fast marching method proposed by Mirebeau (2014b).

Finally, we provide two types of image datadriven speed functions that are computed by steerable filters. These speed functions are therefore orientation dependent, by which we apply the proposed Finsler elastica minimal path model to interactive closed contour detection, perceptual grouping and tubular structure extraction.
More contributions have been added regarding the original conference paper presented in (Chen et al, 2015), such as the applications of interactive closed contour detection for image segmentation and perceptual grouping, respectively.
1.3 Paper Outline
In Section 2 we briefly introduce the existing minimal path models, the concept of Finsler metric, and algorithms for distance computation and path extraction. The relationship between the Euler elastica bending energy and the Finsler elasica metric is analyzed in Section 3. In Section 4 we introduce two datadriven speed functions which are dependent of orientations. The applications of the Finsler elastica minimal path model are presented in Section 5. Experiments and Conclusion are presented in Sections 6 and 7, respectively.
2 Background on the Minimal Path Models
2.1 CohenKimmel Minimal Path Model
Let () denote the image domain. The classical CohenKimmel minimal path model (Cohen and Kimmel, 1997) was designed to find the global minimum of the geodesic energy which measures the length of Lipschitz paths as follows
(1) 
where denotes the canonical Euclidean norm. is a positive constant and is a potential function which usually depends on the image gradient magnitudes or intensities as suggested by Cohen and Kimmel (1997).
A geodesic , joining the source point to a point , is a path that globally minimizes the length (1):
where is the collection of all Lipschitz paths with and .
The minimal action map associated to a source point , is the global minimum of the length (1) among all paths involved in the collection :
(2) 
The CohenKimmel minimal path model considers an isotropic Riemannian metric, which is defined by
(3) 
where denotes the scalar product over and is a symmetric positive definite tensor field which is proportional to the identity matrix :
(4) 
The minimal action map defined in (2) satisfies the following Eikonal PDE (Cohen and Kimmel, 1997):
(5) 
where is the Euclidean gradient of with respect to the position in the domain.
A geodesic parameterized by its arclength can be obtained through the gradient descent ordinary differential equation (ODE) till reaching the point :
(6) 
where denotes positive collinearity. Then the geodesic is obtained by reversing the path with and .
2.2 Anisotropic Riemannian Metric Extension
Bougleux et al (2008) and Jbabdi et al (2008) extended the isotropic Riemannian metric (3) to the anisotropic Riemannian case invoking a symmetric positive definite tensor field
(7) 
where is the dimension of the domain and the vector fields are related to the image data. are the image datadriven functions associated to .
Based on the tensor field (7), the anisotropic Riemannian metric can be expressed as:
(8) 
The minimal action map associated to the anisotropic Riemannian metric satisfies the Eikonal PDE:
(9) 
where we denote that .
The geodesic is obtained by reversing the geodesic which is the solution to the following ODE
(10) 
Benmansour and Cohen (2011) proposed an anisotropic radiuslifted Riemannian metic, as a special case of the metric , to address the problem of extracting both the centerline positions and the thickness values of a tubular structure simultaneously.
2.3 Isotropic orientationlifted Riemannian Metric Extension
In order to take into account the local orientation in the image, it is possible to include orientation information in the energy minimization. For this purpose, the image domain can be extended to an orientationlifted space by product with an abstract orientation space (Kimmel and Sethian, 2001), i.e., and the problem is to find a geodesic in the new lifted space . Each point in the orientationlifted path is thus a pair composed of a point in the image domain and an orientation , i.e., .
For any point and any vector , where denotes the orientation variation, the isotropic orientationlifted Riemannian metric is expressed by
(11) 
where is an orientationdependent speed function relying on the image data and is a constant.
Following the form of (3), the metric defined in (11) can be written by using a tensor field
(12) 
where is a diagonal matrix. Note that can be extended to a positive scalar function that is dependent of the image data. For simplicity, we set to be a positive constant in this paper, as suggested in (Kimmel and Sethian, 2001; Péchaud et al, 2009).
The curve length of an orientationlifted curve with respect to the isotropic orientationlifted Riemannian metric can be expressed as
(13) 
where , for all .
The idea of orientation lifting was applied to tubular structure extraction by Péchaud et al (2009) with additional radius lifting, typically accounting for the radius of a tubular structure in the processed image. Orientation lifting often improves the results from (Cohen and Kimmel, 1997; Li and Yezzi, 2007), but suffers from the fact that for the curve length of an orientationlifted curve , nothing in (2.3) constrains the path direction to align with the orientation vector associated to , where . In other words, this isotropic orientationlifted Riemannian metric cannot penalize the curvature of the geodesic, a point which is addressed in this paper.
2.4 General Minimal Path Model and Finsler Metric
The general minimal path problem is posed on a bounded domain equipped with a metric depending on positions and orientations . This metric defines at each point an asymmetric norm
(14) 
These norms must be positive whenever , homogeneous, and obey the triangular inequality. However, in general we allow them to be asymmetric: such that
(15) 
One can measure the curve length of a regular curve with respect to the metric :
(16) 
The minimal action map is defined by:
(17) 
which is the unique viscosity solution to an Eikonal PDE (Lions, 1982; Sethian and Vladimirsky, 2003):
(18) 
where is the dual norm of defined for all by
(19) 
Based on the definition of the dual norm in (19), the corresponding optimal direction map is then obtained by
(20) 
Again, the geodesic is obtained by reversing the geodesic with and , where is tracked through the following ODE involving the minimal action map and the optimal direction map
(21) 
Numerically, the ODE expressed in (21) is solved by using Heunâs or RungeKutta’s methods, or more robustly using the numerical method proposed by Mirebeau (2014a).
The metric considered in this paper combines a symmetric part, defined in terms of a symmetric positive definite tensor field , and an asymmetric part involving a vector field :
(22) 
for all and any vector . The fields and should obey the following smallness condition to ensure that the Finsler metric is positive:
(23) 
Equation (22) defines an anisotropic Finsler metric in general. This is an anisotropic Riemannian metric if the vector field is identically zero, and an isotropic metric if in addition the tensor field is proportional to a diagonal matrix like in (12). Therefore, the general Eikonal PDE (18) and the geodesic back tracking ODE (21) reduce respectively to equations (5) and (6) for the isotropic case, and respectively to equations (9) and (10) for the anisotropic case.
2.5 Computation of the Minimal Action Map
In order to estimate the minimal action map , presented in (17) and (18), a discretization grid of the image domain is introduced  or of the extended domain in the case of an orientationlifted metric. For each point , a small mesh of a neighbourhood of with vertices in is constructed. For example, can be the square formed by the four neighbours of in the classical fast marching method (Sethian, 1999) on a regular orthogonal grid used in (Cohen and Kimmel, 1997). In contrast with Sethian’s classical fast marching method which solves the discrete approximation of the Eikonal PDE itself, an approximation of the action map , with the initial source point , is calculated by solving the fixed point system (Tsitsiklis, 1995):
(24) 
where the involved HopfLax update operator is defined by:
(25) 
where denotes the piecewise linear interpolation operator on the mesh , and lies on the boundary of . The equality , replacing in (24) the Eikonal PDE: of (18), is a discretization of Bellman’s optimality principle, which is similar in spirit to the Tsitsiklis approach (Tsitsiklis, 1995). It reflects the fact that the minimal geodesic , from to , has to cross the mesh boundary at least once at some point ; therefore it is the concatenation of a geodesic from to , which length is approximated by piecewise linear interpolation, and a very short geodesic from to , approximated by a segment of geodesic curve length . The dimensional fixed point system (24), with the number of grid points, can be solved by a single pass fast marching method as described in Algorithm 1 in Appendix.
The classical fast marching methods (Sethian, 1999; Tsitsiklis, 1995) using the square formed neighbourhood have difficulty to deal with the computation of geodesic distance maps with respect to anisotropic metrics, especially when the anisotropy gets large. An adaptive construction method of such stencils was introduced in (Mirebeau, 2014a) for anisotropic 3D Riemannian metric, and in (Mirebeau, 2014b) for anisotropic 2D Finsler metric, providing that the stencils or mesh at each point or satisfies some geometric acuteness property depending on the local metric . Such adaptive stencils based fast marching methods lead to breakthrough improvements in terms of computation time and accuracy for strongly anisotropic geodesic metrics. When the above mentioned geometric properties do not hold, the fast marching method is in principle not applicable, and slower multiple pass methods must be used instead, such as the Adaptive Gauss Siedel Iteration (AGSI) of Bornemann and Rasch (2006). The present paper involves a 3D Finsler metric (37), for which we constructed stencils by adapting the 2D Finsler metric construction method proposed by Mirebeau (2014b). Although these stencils lack the geometric acuteness condition, we found that the fast marching method still provided good approximations of the paths, while vastly improving computation performance. In Section 3.3, the complexity will be discussed more.
Note that whenever we mention fast marching method in the next sections, we mean the fast marching method with adaptive stencils proposed by Mirebeau (2014b).
3 Finsler Elastica Minimal Path Model
In this section, we present the core contribution of this paper: the orientationlifted Finsler metric embedding curvature penalty term, defined over the orientationlifted domain , where ) denotes the angle space with periodic boundary conditions and denotes the image domain.
3.1 Geodesic Energy Interpretation of the Euler Elastica Bending Energy via a Finsler Metric
The Euler elastica curves were introduced to the field of computer vision by Nitzberg and Mumford (1990) and Mumford (1994). They minimized the following bending energy:
(26) 
where is a regular curve with nonvanishing velocity, is the curvature of curve , is the Euclidean curve length of , and is the arclength. Parameter is a constant. is an image datadriven speed function.
For the sake of simplicity, we set , yielding the simplified Euler elastica bending energy
(27) 
where the general case will be studied in Section 3.4.
The goal of this section is to cast the Euler elastica bending energy (27) in the form of curve length with respect to a relevant asymmetric Finsler metric. We firstly transform the elastica problem into finding a geodesic in an orientationlifted space. For this purpose, we denote by
(28) 
the unit vector corresponding to .
Let be a regular curve with nonvanishing velocity and be its canonical orientation lifting. For any , is defined as being such that:
(29) 
According to the definition of in (28), one has
(30) 
where denotes the vector that is perpendicular to a vector . It is known that
(31) 
where is the curvature of path . Then we have the following equations:
Thus, using (29) and (31), we have
(32) 
which yields to
(33) 
Using equations (27), (29) and (33), one has
(34) 
where the Euclidean arclength is defined as
By the definition of , for any we have and
(35) 
where we define the Finsler metric on the orientationlifted domain by
(36) 
for any orientationlifted point , any vector in the tangent space, and where denotes positive collinearity. Note that any other lifting of would by construction of (36) have infinite energy, i.e., .
3.2 Penalized Asymmetric Finsler Elastica Metric
The Finsler metric (36) is too singular to compute the global minimum of (27) by directly applying the numerical algorithm such as the fast marching method (Mirebeau, 2014b). Hence we introduce a family of orientationlifted Finsler elastica metrics over the orientationlifted domain , depending on a penalization parameter as follows:
(37) 
for any orientationlifted point and any vector , and where is the unit vector associated to which denotes the position of in the orientation space .
As , the penalized Finsler elastica metric can be expressed as:
(38) 
The term will vanish if vector is positively proportional to . Therefore, one has for any and any
The metric (37) has precisely the required form formulated in (22), with tensor field as:
(39) 
and vector field
(40) 
for any . Note that the definiteness constraint (23) is satisfied:
The anisotropy ratio characterizes the distortion between different orientations induced by a metric. Letting , and , the anisotropy ratio of the Finsler elastica metric (37) can be defined by:
(41) 
where the norm . As an example, for the Finsler elastica metric (37) with and , we can show that in (41) is obtained by choosing and , so that .
Moreover, one can define the physical anisotropy ratio of the Finsler elastica metric by replacing by and the variables and in (41), respectively. In this case, for any , the physical anisotropy ratio is equal to and only depends on .
A crucial object for studying and visualizing the geometry distortion induced by a metric is Tissot’s indicatrix defined as the collection of unit balls in the tangent space. For point and , we define the unit balls for metrics and respectively by
(42) 
and
(43) 
Then any vector can be characterized by
(44) 
where we introduce and as follows:
Using (44), one has
(45) 
Thus is a flat 2D ellipse embedded in the 3D tangent space, and containing the origin on its boundary. Particularly, when , turns to a flat 2D disk of radius as shown in Fig. 2a.
On the other hand, when , a short computation shows that any vector is characterized by a quadratic equation
(46) 
where are all . Hence is an ellipsoid, for instance see Fig. 2b with , almost flat in the direction of due to the large factor , which converges to the flat disk in the Haussdorf distance as .
Tissot’s indicatrix is also the control set in the optimal control interpretation of the Eikonal PDE (18). The Haussdorf convergence of the control sets guarantees that the minimal action map and minimal paths for the metric converge towards those of as . Elements of proof of convergence can be found in Appendix B.
3.3 Numerical Implementations
Numerically, anisotropy is related to the problem stiffness, hence to its difficulty. In Table 1, we show the computation time and the average number of HopfLax updates required for each grid point by the adaptive stencils based fast marching method (Mirebeau, 2014b) for and different values of on a grid. This experiment was performed with a C++ implementation running on a standard 2.7GHz Intel I7 laptop with 16Gb RAM.
We observe on Table 1 a logarithmic dependence of computation time and average number of the HopfLax updates per grid point with respect to anisotropy. These observations agree with the complexity analysis of the fast marching method presented in (Mirebeau, 2014b), yielding the upper bound , depending polylogarithmically on the anisotropy ratio (41), and quasilinearly on the number of discretization points in the orientationlifted domain . In contrast, numerical methods such as (Sethian and Vladimirsky, 2003) displaying a polynomial complexity in the anisotropy ratio would be unworkable. The iterative AGSI method (Bornemann and Rasch, 2006), on the other hand, requires hundreds of evaluations of the HopfLax operator (25) per grid point to converge for large anisotropies, which also leads to prohibitive computation time, thus impractical. For or , the average numbers of the HopfLax updates per grid required by the AGSI method are approximately and , respectively, while the numbers of HopfLax from the fast marching method are only and , respectively, as demonstrated in Table 1.
In Fig. 3a, we show different Finsler elastica minimal paths, computed by the fast marching method (Mirebeau, 2014b), with (see (37)) and different values of . The arrows indicate the initial source and end points tangents. The cyan point denotes the initial position and the blue point indicates the end position. In Figs. 3b and 3c, we show the Finsler elastica minimal paths for different values of , with and respectively. In this experiment, the angle resolution is and the image size is . When , the metric is constant over the domain and degenerates to the isotropic orientationlifted metric (11), since the coefficient in front of the term in (37) vanishes. Hence the minimal geodesics are straight lines, see Fig. 3a, that do not align with the prescribed endpoints tangents. From Fig. 3, one can point out that as and increasing, curvature penalization forces the extracted paths to gradually align with the prescribed endpoints tangents and take the elastica shape.
1  10  20  30  100  200  1000  

time  13.9s  25.3s  27.3s  27.7s  31.7s  33.9s  36.8s 
number  3  5.49  6.06  6.49  7.27  7.82  8.12 
3.4 Image DataDriven Finsler Elastica Metric
We use in Section 3.1 for the sake of simplicity. In the general case, the metric (36) and its approximation (37) should be respectively replaced by and . Furthermore, in order to take into account the orientation information, we use an orientationdependent speed function to replace . In this case, the datadriven Finsler elastica metric can be defined by
(47) 
and minimizing Eq. (26) is approximated for large values of by minimizing
where is the orientationlifted curve of . The datadriven Finsler elastica metric is asymmetric in the sense that for most vectors , one has
(48) 
This asymmetric property can help to build a closed contour passing through a collection of orientationlifted points as discussed in Section 5.1.
The minimal action map associated to the datadriven Finsler elastica metric and an initial source point is the unique viscosity solution to the Eikonal PDE (18) (Lions, 1982; Sethian and Vladimirsky, 2003):
(49) 
where is the dual norm of defined by (19). The geodesic distance value between and depends on both the curvature and the orientation information of the minimal paths. When is sufficiently large, the spatial and angular resolutions are sufficiently small, the fixed point system (24) is properly solved, and the minimal paths are properly extracted by (21).
4 Computation of DataDriven Speed Functions by Steerable Filters
In this section, we introduce two types of anisotropic speed functions for boundary detection and tubular structure extraction, both of which are based on the steerable filters.
4.1 Steerable Edge Detector
Jacob and Unser (2004) proposed a new class of edge detection filters based on the computational framework and the steerable property. Letting be a 2D isotropic Gaussian kernel with variance and , the computational steerable filter with order can be expressed as
(50) 
where and are the orientationdependent coefficients which can be computed in terms of some optimality criteria. Particularly, when , the steerable filter becomes the classical Canny detector (Canny, 1986). For higher order steerable filters, for example, or , the orientationdependent responses of the filters will be more robust to noise. Therefore, we set for the relevant experiments. Regarding the details of the computation of , we refer to (Jacob and Unser, 2004).
A color image is regarded as a vector valued map , for each . A multiorientation response function of a color image can be computed by the steerable filter (50)
(51) 
For a gray level image , we have the simple computation of :
(52) 
The multiorientation response function is symmetric in the sense that for any orientation , one has
4.2 MultiOrientation Optimally Oriented Flux Filter
Optimally oriented flux filter is used to extract the local geometry of the image (Law and Chung, 2008). The oriented flux of an image , of dimension , is defined by the amount of the image gradient projected along the orientation flowing out from a 2D circle at point with radius :
(53) 
where is a Gaussian with variance . is the outward normal of the boundary and is the infinitesimal length on . According to the divergence theory, one has
for some symmetric matrix :
(54) 
where is an indicator function of the circle .
Let be the eigenvalues of symmetric matrix and suppose that the intensities inside the tubular structures are darker than the background regions. If a point is inside the tubular structure, one has and , where is the optimal scale map
(55) 
The scale normalized factor in (55) provides the responses of the optimally oriented flux filter a scale invariant property as discussed in (Law and Chung, 2008).
As shown in (Benmansour and Cohen, 2011), the optimally oriented flux filter is steerable, which means that we can construct a multiorientation response function for any by:
(56) 
where is a unit vector associated to .
In addition, the vesselness map , which indicates the probability of a pixel belonging to a vessel, can be calculated by:
(57) 
The vesselness map will be used to compute the isotropic Riemannian metric in the following section.
4.3 Computation of DataDriven Speed Functions
A requirement on the speed function used by the datadriven Finsler elastica metric defined in (47) is that it should depend on the position and the orientation. Therefore, based on the orientationdependent response function , defined in Section 4.1, the speed function that is used for object boundary detection can be computed by
(58) 
for any and any orientation , where is defined as .
Similarly, we define the speed function for tubular structure extraction by using the function (56)
(59) 
where and are positive constants. In this paper, we use for all the experiments.
5 Closed Contour Detection and Tubular Structure Extraction
We use the following convention in the remaining part of this paper: if is a point in , then we use to denote the orientationlifted point that has the same physical position with but the opposite direction, where .