3D curve regularization

In this paper, we study the regularization of 3D curves connecting two points. We propose an energy-based formulation which is an extension to 3D of the geodesic active contours introduced in 2D by Caselles et al. in 1997. By minimizing this energy we try to minimize the curve length but keeping the curve close to the original one. The energy depends on a regularization parameter which determines the smoothness of the regularized curve. We show the Euler-Lagrange equation of the proposed energy using the arc-length parameterization of the curve. We interpret the Euler-Lagrange equation in terms of the Frenet–Serret frame and we study some qualitative properties of the energy minima. We apply the steepest-descent method to approximate the local minima of the energy using an evolution equation. We propose a numerical scheme to solve the evolution equation and we present some experiments on real data in the context of aortic centerline regularization.


Introduction
In the context of medical applications, the development of new image modalities and medical treatments is a source of interesting mathematical problems. The current available 3D scans of the human body allow the implantation of stents in the vascular system. A critical issue in stent implantation is the length of the stent required for the treatment. This length is usually measured based on the length of the centerline of the artery between 2 points, but the estimation of this centerline is usually noised which can provide an inaccurate estimation of the length of the stent. Inspired by this practical problem we study in this paper the regularization of 3D curves connecting 2 points. We will use an extension to 3D of the Qualitative Properties of Nonlinear Partial Differential Equations, dedicated to Professor Ildefonso Díaz on the occasion of his 70th birthday.
LetC : [0, 1] → R 3 be a 3D curve satisfying the following boundary condition: for some p 0 , p 1 ∈ R 3 . We consider the minimization problem given by the energy: where g : R 3 → R,C q (q) is the tangent vector to the curveC at the pointC(q) given by the first derivative ofC(q) at q, and . is the L 2 norm. In the context of 3D curve regularization, given an original curveC 0 , we define p 0 =C 0 (0), p 1 =C 0 (1) and g as where d(x,C 0 ) is the Euclidean distance of a 3D point x toC 0 and w ≥ 0 is a regularization parameter. Notice that for all w ≥ 0 the function g(.) is nonnegative and in the case w = 0, the global minimum of the energy E(C) is attained inC ≡C 0 (E(C 0 ) = 0). However, if w > 0, then, a regularization effect is introduced because we penalize high values of the curve length and the minima of energy E(C) are expected to be a regularization ofC 0 . The larger the value of w, the stronger the regularization effect. Equation (2) is an extension to 3D of the geodesic active contour formulation in 2D, introduced by Caselles et al. in [1]. As shown in [1], E(C) corresponds to a new curve length definition obtained by weighting the Euclidean element of length ds = C q (q) dq by g(C(q)). Notice that if we denote by C(s), the curveC(q) reparameterized using the arc-length s, then where |C | is the curve length. In this paper, we will show that if C(s) is the arc-length reparametrization of a local minimum of E(C), then C satisfies the Euler-Lagrange equation with C(0) = p 0 and C(| C |) = p 1 . ∇g represents the gradient of the function g, < ., . > the scalar product of vectors, and C s and C ss represent, respectively, the first and second derivative of the curve C with respect to the arc-length parameter s. We will use this Euler-Lagrange equation to approximate numerically the local minima of E(C) and we will present some experiments on real data in the context of aortic centerline regularization. The rest of the paper is organized as follows: in Sect. 2, we present some related works. In Sect. 3, we study the Euler-Lagrange equation of the energy E(C) and its interpretation in terms of the Frenet-Serret frame. In Sect. 4, we present some qualitative properties of the minima of E(C). In Sect. 5, we present a numerical scheme to approximate local minima of E(C) using the steepest-descent method. In Sect. 6 we present some numerical results in the context of aortic centerline regularization using real data. Finally, in Sect. 7, we draw some conclusions.

Related work
The energy E(C) is an straightforward extension to 3D of the geodesic active contours introduced by Caselles et al. in 2D [1]. The elegant energy formulation of the geodesic active contours have been extensively used in the literature in the last years. However, in the original geodesic active contours, the authors study a completely different problem: the detection of object contours. In particular, the shape of the function g in the energy E(C) is completely different to the one we propose in the Eq. (3). On the other hand, in the geodesic active contours, the curves are closed and the authors do not use the boundary condition (1) to connect two fixed points. An extension of the geodesic active contours using the boundary condition (1) in 2D has been proposed in [2].
In [3] and [4], the authors study the problem of 3D curve smoothing using the evolution equation: where N is the normal and k the curvature. This is an heuristic formulation of the smoothing process which does not take into account the correct formulation of the Euler-Lagrange equation of energy E(C) that we study in this paper.
In [5] and [6], some energies for curvature penalized minimal paths are proposed to regularize 2D curves. In [7], the authors propose a minimal path approach for tubular structures segmentation in 2D images with applications to retinal vessel segmentation. In [8] , the authors introduce a method for the automatic estimation of the aorta segmentation and the centerline estimation.

Euler-Lagrange equation of E(C)
Theorem 1 Let us consider a family of curvesC : Proof : we follow the technique showed in [1] (Appendix B) adapted to the 3D case. By deriving E(C) with respect to u, we obtain: integrating by parts in the second term and taking into account the boundary conditions we obtain: Taking into account that d dq we obtain the theorem statement, that is: As a consequence of this theorem, we obtain that the Euler-Lagrange equation of the energy E(C) is given by the Eq. (5). That is Next, we are going to reformulate this equation in terms of the Frenet-Serret frame which associates to each point of the curve, the orthonormal basis T(s), N(s), and B(s) defined as follows (see, for instance, [9]) : is the unit vector tangent to the curve (notice that as s is the arc-length parameterization of the curve, then C s (s) ≡ 1 and T(s) = C s (s)).
is the normal unit vector. (κ(s) = dT ds (s) is the curvature which, intuitively, measures how far is the curve to be a straight line.
• B(s) = T(s) × N(s) is the binormal unit vector (the cross product of T(s) and N(s)).
In the case that C ss (s) = 0, we have that the Frenet-Serret frame is well-defined, C ss (s) = κ(s)N(s) and then, using that we can express the Euler-Lagrange equation in terms of the Frenet-Serret frame as:

Some qualitative properties of the minima of E(C)
First we show that if the energy E(C) decreases with respect to the energy for the original curve, then the length of the curve also decreases.
Proof This result is a straightforward consequence of the following inequality: Next, we study the asymptotic state of the minima of E(C) when w → ∞.
Proposition 3 LetC(q)= p 0 (1 −q) + p 1 q be the curve corresponding to the straight segment joining p 0 and p 1 In particular, when w → ∞,C w approximates the segment joining p 0 and p 1 .

Approximation of the local minima of E(C) by the steepest-descent method
According to the result of theorem 1 and the steepest-descent method, to connect an initial curve C 0 with a local minimum of the energy E(C) we solve the evolution equation: C t (t,s) = g(C(t,s))C ss (t,s) − (∇g(C(t,s))− < ∇g(C(t,s)), C s (t,s)> C s (t,s)), (10) the asymptotic state of this evolution equation when t → ∞ corresponds to a local minima of E(C). Notice that the vector ∇g(C(s))− < ∇g(C(s)), C s (s) > C s (s) is the projection of the vector ∇g(C(s)) in the plane orthogonal to C s (s). Moreover, if C ss (s) = 0, then C ss (s) = κ(s)N(s) which is also orthogonal to C s (s), therefore the evolution equation makes always move the curve C in the orthogonal plane to the tangent vector C s . In the curve evolution Eq. (10) we use as initial value C(0, s) = C 0 (s) and the boundary condition C(t, 0) = p 0 and C(t, | C(t) |) = p 1 . Therefore we look for a local minima of E(C) close to te original curve C 0 . Notice that we always initialize C(0, s) as the original curve and we do not have to deal with the problem of the curve initialization which appears, for instance, in object segmentation using active contours. We also emphasize that the energy E(C) can have, in general, several local minima and that global minima could be quite different to the one we obtain looking for the local minima closer to the original curve C 0 . The complexity of the structure of the local minima of E(C) strongly depends on the complexity of the geometry of the initial curve C 0 . For instance, if C 0 is an straight segment joining p 0 and p 1 , then, C 0 is the global minima of E(C). However, for a general curve, the situation can be much more complex. Assume, for instance, that w is small, p 0 and p 1 are very close but | C 0 | is much larger than p 1 − p 0 , then we can expect that the global minima is close to the segment joining p 0 and p 1 and that this global minima is different of the local minima closer to C 0 . In any case, the fact that we use always the original curve as initialization of the steepest descent method strongly reduces the possibility to get trapped in spurious local minima.
We use a basic explicit finite difference scheme to approximate numerically the solution of the evolution Eq. (10) with the initial value C(0, s) = C 0 (s) and the boundary condition C(t, 0) = p 0 and C(t, | C(t) |) = p 1 . Numerically, a 3D curve C is given by an ordered collection of 3D points {C i } i=1,..,χ (C) , where χ(C) represents the number of points of the curve. By the boundary condition we always have that C 1 = p 0 and C χ(C) = p 1 . Notice that in the evolution Eq. (10) the curves are parametrized using the arc-length, therefore we assume that which means that we use a curve parameterization with a constant arc-length h > 0. In what follows, we assume, without loss of generality (by scaling the size of C i ) that h = 1. Given a discretization step dt, we denote by C n,k an approximation of C(n · dt, k). We compute C n,k using the following explicit finite difference scheme: for k = 2, .., χ(C) − 1 and n = 0, 1, ... C n,k s and C n,k ss are approximations of C s (n · dt, k) and C ss (n · dt, k) computed as: We use some algorithms proposed in [4] to deal with two technical issues related with this scheme: the first one is that we need to compute g(C n,k ) = d(C n,k , C 0 ) + w and ∇g(C n,k ), which requires the computation of the 3D distance to a curve. The second one is that after each iteration we have to reparametrize the discrete curve to preserve the condition C n,k − C n,k−1 = 1 for all n and k = 2, .., χ(C) − 1. In both cases we use algorithms proposed in [4] to deal with these issues.
We use as stopping criterion of the iterative scheme (12) the condition

Numerical experiments
We will present some experiments in the context of centerline regularization of blood vessels extracted from a 3D Computed tomography (CT) scan. There are different types of algorithms to extract the centerline of blood vessels (see, for instance, [8,10]). These algorithms segment the vessel lumen and the vessel centerline is obtained from this segmentation using, for example, the centers of the maximal spheres included in the vessel lumen. The method that we propose can be used to regularize the centerline extracted by any existing method in the literature. It can be considered a post-processing of the obtained centerline. In our knowledge this approach is new and until now the problem of regularizing the 3D curves obtained by these centerlines had not been studied. Centerline regularization is necessary because normally its estimation may include irregularities caused by poor CT quality, or because the patient's pathologies, such as aneurysms or thrombi, distort the geometry of the vessel lumen and therefore negatively affect the centerline estimation. Correct estimation of the centerline is very important in applications such as stent implantation. Since the stent is a regular and elastic object, the correct estimation of the length of the stent requires an estimation of the centerline through a regular curve.
In our experiments, we use, as the original curve C 0 , the centerline of an aortic lumen computed from a real CT scan using the technique introduced in [8]. The main goal of these experiments is to explore the evolution of the energy, E(C n ), and the length, | C n |, of the curve under the action of the iterative scheme (12), and the influence of w in the results. We denote by C ∞ the final estimate of the iterative scheme (12) at the stopping time t ∞ = n · dt obtained using (13).
In Fig. 1, we compare the curves C 0 and C ∞ for different values of w. As expected, we observe that the larger the value of w, the stronger is the regularization of the curve. For w = 10 and w = 25 we can not appreciate visually in Fig. 1 any difference between the curves C 0 and C ∞ . To visualize such differences we show in Fig. 2 some images corresponding to zoom of the curves around some particular points. We can clearly appreciate in this figure the regularization effect of the minimization of the energy E(C). Fig. 3 We show, for w = 10, 25, 50, 100, the evolution of the energy E(C n ) and the length of the curve C n , solution of the iterative scheme (12), using dual axis charts (on the left the energy and on the right the length).
In the horizontal axis we show the value of t = n · dt In Fig. 3 we show, for different values of w, the evolution of the energy E(C n ) and the length of the curve C n , solution of the iterative scheme (12). We observe a nice convergence behavior of the iterative scheme: the energy and the length of the curve decrease sharply at the beginning, then the energy stabilizes and finally the stopping time t ∞ is reached using criterion (13).
In Table 1, we show, for different values of w, the initial and final values of the energy, the length of the curve (in millimeters) and the stopping time t ∞ of the iterative scheme. We observe a significant reduction of the value of the energy E(C). The reduction of the length of the curve and the stopping time strongly depend on the value of w.
An important practical issue in the applications is the choice of the regularization parameter w. To illustrate the influence of the parameter w in the curve regularization process we have used, in the experiments, high values of w and we have observed that the result obtained is the expected one taking into account the mathematical formulation of the model. In real applications, such as the regularization of blood vessels, the value of w that would be taken is, usually small, below w = 10, since what is required, in general, is to obtain a regular curve but close to the original one. In any case, w is a parameter of the algorithm that will have to be adjusted according to the noise level of the original centerline. As explained above, the noise that the centerline presents can depend on the quality of the CT scan, the pathologies that the patient presents or the particular algorithm used to estimate the centerline.

Conclusion
3D curve regularization is an interesting mathematical problem which appears in a natural way in medical applications. To deal with this problem, we extend to 3D the elegant energy formulation of the geodesic active contours introduced by Caselles et al. in [1]. We show the Euler-Lagrange equation of the energy, some qualitative properties of the minima and we use the steepest-descent method to design an iterative numerical scheme to approximate local minima of the energy based on the associated evolution equation. We present some experiments on real data which show that the proposed iterative scheme works quite well. The energy value decreases sharply across iterations until the convergence of the algorithm.
We also found that, as expected, the regularization parameter w determines the regularity and length of the final curve obtained as the asymptotic state of the iterative scheme. The larger the value of w, the smoother the regularized curve and the shorter its length.

Supplementary information
We provide an ASCII text file with the point coordinates of the 3D curve used in the experiments corresponding to the aortic centerline extracted from a CT scan.
article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.