# Iterative Methods for Obtaining Energy-Minimizing Parametric Snakes with Applications to Medical Imaging

Authors: Alexandru Ioan Mitrea; Radu Badea; Delia Mitrea; Sergiu Nedevschi; Paulina Mitrea; Dumitru Mircea Ivan; Octavian Mircia Gurzău
Journal: Computational and Mathematical Methods in Medicine
 (2012)

---

## Abstract

After a brief survey on the parametric deformable models, we develop
an iterative method based on the finite difference schemes in order to obtain energy-minimizing snakes. We estimate the approximation error, the residue, and the truncature error related to the corresponding algorithm,
then we discuss its convergence, consistency, and stability. Some aspects
regarding the prosthetic sugical methods that implement the above numerical methods are also pointed out.


---

## Body

## 1. Introduction

The deformable models represent a powerful researched model-based approach to computer-assisted medical image analysis, their applications in this framework including image segmentation, shape representation and motion tracking. The theory of deformable models is an interdisciplinary scientific domain, which has appeared and developed in the last two decades, in strong connection with practical problems of medicine, image processing, and physics. This theory joins methods, results, and techniques of various mathematical fields, physics and mechanics. The mathematical foundation of this theory represents the confluence of Functional Analysis, Approximation Theory, Differential Equations, Differential Geometry, Calculus of Variations, Numerical Analysis, Linear Algebra, and Probability Theory. The ancestors of the deformable models, in classical sense, are considered Fischler and Elschlager, with their spring-loaded templates, [1], together with Widrow [2], with its rubber mask technique [2, 3].The theory of deformable models, in its modern form, originates from the general theory of continuous multidimensional deformable models in a Lagrangian dynamic of Terzopoulos (1987) [4]. In fact, the deformable curves (2D models) and the deformable surfaces (3D models) gained popularity after their use in computer vision by Kaas et al. [5] and in computer graphics, by Terzopoulos and Fleischer [6] in the mid-1980s. Since then, the deformable models, known also as active contour models or snakes, have been extensively used for many applications both in 2D and 3D.Two general types of deformable models have been developed: firstly, the parametric or variational models, which originate from the papers of Kaas et al. [5] and are based on the minimization of the energy-functional associated to the model, and secondly, the geometric models, which were introduced independently by Caselles et al. [7] and Malladi et al. [8], and are based on the front propagation theory [9].A good survey on deformable models and their applications can be found in [10, 11]. Recent contributions on parametric deformable models have appeared in the papers [12, 13]. On the general topic of numerical methods applied in medical imaging, the recent papers [14, 15] must be mentioned.In this paper, we deal with the deformable parametric models. The basic goal of the theory of parametric deformable models is to determine the energy-minimizing 2D or 3D models, namely, the curves or surfaces which minimize the corresponding energy functional. Two approaches will point out in order to obtain the optimal model. The first approach is based on the Euler-Lagrange-Poisson (ELP) and Euler-Gauss-Ostrogradski (EGO) equations of Calculus of Variations in order to minimize the energy-functional. The second one (the classical approach) consists of using reconstruction methods, such as the interpolation of the sparse data extracted from the image, in order to obtain a representation of the original data. In what follows we develop methods and techniques related to the first approach. Generally, the energy-functional is not convex, so it may have many local minimum. On the other hand, the analytic solution of (ELP) equation has a complicated form or it is inaccessible explicitly. Therefore, a practical and strong approach for finding local minimum of the energy functional is to construct a dynamic system that is governed by the energy functional and allow the system to evolve to the equilibrium state. Dynamic models are valuable for medical image analysis, because most anatomical structures are deformable and continually undergo nonrigid motion “in vivo.” In fact, the user is interested to find a good 2D or 3D contour in a given area. Consequently, a rough prior estimation of the 2D or 3D model is provided, then this initial model undergoes a deformation until reaching a local minimum of the energy functional. This deformation process can be achieved in one of the following ways:(1)
in a Hamiltonian-type approach, by performing a strictly decreasing energy path, for example, via dynamic programming methods [16, 17];(2)
in a Lagrange-type approach, by applying the mechanical principles of Lagrange [3, 18];(3)
by using a friction force, in order to constrain the displacement of the snake [5];(4)
by using the (ELP) evolution equation, associated to the initial (ELP) equation [19].In this paper, we shall adopt the method of the evolution equation. So, a prior estimation of the deformable surface is provided, then it is refined step by step, based on the (EGO) equation and using discretization methods.The paper outline is as follows. The next section is devoted to present 2D and 3D energy-minimizing models, both in their static and dynamic forms. The method for reducing the 3D problem to a 2D modeling is also pointed out, in order to minimize the computational costs of the numerical methods. The third section contains the main theoretical result of the paper. Based on finite difference schemes of explicit type, we derive an (ELP) algorithm for obtaining an energy-minimizing snake in its approximated form, then we estimate its approximation error and we discuss its consistency, convergence, and stability. The last section deals with the behaviour of prosthetic surgical methods and prosthetic medical materials, based on Software tools, which implement the iterative methods developed in the previous sections.


## 2. Energy-Minimizing Models

### 2.1. Energy-Minimizing Snakes (2D Models)

From mathematical point of view, a 2D parametric deformable model (usually known assnake) is provided by a family 𝒜 of parametrized smooth curves satisfying given boundary conditions and an associated energy-functional. More exactly, denote by C2([0,1],ℝ2) the space of all vectorial functions v=(x,y)T so that the scalar functions x=x(s) and y=y(s), 0≤s≤1 are continuous together with their derivatives up to the second-order on the standard interval [0,1], that is, x,y∈C2[0,1]; obviously, we can consider an arbitrary compact interval [a,b] of the real axis instead of [0,1]. The family 𝒜 of admissible deformations consists of all parametrized curves (snakes):(1)(γ):v(s)=(x(s),y(s))T,0≤s≤1,v∈C2([0,1],R2),
such that the values v(0),v(1),v′(0), and v′(1) are given; we adopt the notation |v|2=|x|2+|y|2.In order to find the optimal position of the snake, it is necessary to characterize its state, by means of anenergy-functional,that is associated to the class 𝒜. Let us consider the following data:(i)
theweight-functionsw1(s) and w2(s), which control the elasticity and the rigidity of (γ), respectively; generally, these are nonnegative scalar functions of class C2[0,1],(ii)
theimage intensity function I=I(x,y), which is a real function of class C2(ℝ2),(iii)
thepotential associated to the external forces, represented by a real function P(v)=P(x,y), of class C2(ℝ2). The simplest useful choice for the potential is P(v)=w3I(v), where w3 is a weight-scalar. The most used choice is P(v)=-λ|∇I(v)|, where λ>0 is a given scalar and ∇=(∂/∂x,∂/∂y)T is the Hamilton (nabla) operator; this choice will be used in this paper, too. Note that P can be defined also by P=Gσ0*I, that is, the Gaussian (variance σ0) filtered image of the input image I [10], and(iv)
the vectorial functionk(v)=(k1(v),k2(v))T of class C1(ℝ2,ℝ2) which control the local dilatation or local contraction of (γ) along its normal; usually, we take k(v)=cv, with c∈ℝ.The shape of the snake(γ) subject to the image I(v) is dictated by the energy functional:(2)E(v)=Eint⁡(v)+Eext(v)+Ebal(v),
where the terms of the right hand of (2) are defined as follows.Theinternal energy(3)Eint⁡(v)=Eels(v)+Erig(v)
is obtained by adding the elastic energy(4)Eels(v)=∫01α(s)|v′(s)|2ds,
and the rigid (bending) energy(5)Erig(v)=∫01β(s)|v′′(s)|2ds.The internal energy characterizes the deformation of a stretchy, flexible snake (contour). The values ofw1(s) and w2(s) show the extent to which the snake can stretch or bend at an arbitrary point (x(s),y(s)) of the snake.Theexternal energy,derived from the image, is given by(6)Eext(v)=∫01P(v(s))ds=-λ∫01|∇I(x,y)|2ds,
and it allows to find the edges in an image so that the snake is attracted to contour with large image gradients.Theballoon energyis an energy of constrained-type, defined as(7)Ebal(v)=-∫01det⁡(k(v),v′)ds=-∫01(k1y′-k2x′)ds.
This energy can be added, optionally, by users, in order to expand (or contract) the snake.Denoting by(8)F(s,v,v′,v′′)=w1(s)|v′2(s)|+w2(s)|v′′2(s)|+P(v(s))-det⁡(k(v),v′),P(v)=-λ|∇I(v)|2,
the following expression of energy-functional E(v) is obtained from (2)–(8):(9)E(v)=∫01F(s,v,v′,v′′)ds.By definition, the triple(𝒜,I,E) is said to be a deformable 2D model (snake).The basic goal of a deformable parametric model is to minimize its energy-functionalE(v), which leads to the energy-minimizing snake. The minimization of the snake energy gives rise to the following vectorial Euler-Lagrange-Poisson (ELP) Equation of Calculus of Variations:(10)∂F∂v-dds(∂F∂v′)+d2ds2(∂F∂v′′)=0.
Now, taking into account the relations (8) and (10), we obtain the vectorial (ELP) equation:(11)2w2(s)viv(s)+4w2′(s)v′′′(s)+2[w2′′(s)-w1(s)]v′′(s)-[2w1′(s)I2+Tr⁡(∇k)J2]v′(s)+∇P(v(s))=0,
where I2=(1001), J2=(01-10), and Tr(A) is the trace of a square matrix A.The scalar (ELP) equations, derived from (11), have the form:(12)2w2xiv+4w2′x′′′+2(w2′′-w1)x′′-2w1′x′-(∂k1∂x+∂k2∂y)y′+∂P∂x=0,2w2yiv+4w2′y′′′+2(w2′′-w1)y′′-2w1′y′+(∂k1∂x+∂k2∂y)x′+∂P∂y=0.
On the other hand, we infer from(8)(13)F(s,v,v′,v′′)=w1(x′2+y′2)+w2(x′′2+y′′2)+P(x,y)-k1y′+k2x′,
which leads to(14)∂2F∂(v′′)2=(2w1,2w2)T.
According to the Legendre conditions and the hypothesis w2>0, the relation (14) proves that any solution of the (ELP) equations (11) or (12) provides a minimum for the energy-functional E(v), namely, an energy-minimizing snake.Example 1.
If we choose in (12) w1=1, w2=0.05, and the boundary conditions v(0)=v(1)=(0,5)T, v′(0)=v′(1)=(0.5,0.5)T we obtain the general solution of the (ELP) equation:
(15)x(s)=C1e3.8042s+C2e-3.8042s+C3e2.3511s+C4e-2.3511s,y(s)=C5e3.8042s+C6e-3.8042s+C7e2.3511s+C8e-2.3511s.
Using boundary conditions, we obtain the graph of the curve (γ) in Figure 1.Figure 1Example 2.
If we choose in (12) w1=1, w2=1, k(v)=10v, I(v)=I(x,y)=(22/3)(x3/2+y3/2), λ=1, and the boundary conditions v(0)=(1,1)T, v′(0)=(3.1,-8.1)T, v(π)=(1-π/10,-2-π/10-3cosh(2π))T, v′(π)=(2.1-cosh(2π),-8.1-2sinh(2π)), then the (ELP) equations have the form:
(16)xiv-x′′-10y′-1=0,yiv-y′′+10x′-1=0,
with the analytical solutions:
(17)x(s)=sin⁡⁡2s+cos⁡⁡2s+sin⁡⁡scosh⁡⁡2s+0.1s,y(s)=-4sin⁡⁡2s+4cos⁡⁡2s+3cos⁡⁡scosh⁡⁡2s,-4sinh⁡⁡2ssin⁡⁡s-0.1s-6.
The graph of the curve (γ) is given in Figure 2.Figure 2


### 2.2. Deformable Dynamic 2D Models

Roughly speaking, the differential fourth-order vectorial equation (11) or the differential eight-order system (12) may have many solutions, which leads to many possible energy-minimizing snakes. As we have seen in the preceding examples (Section 2.1), these solutions have a complicated form; moreover, they are often inaccessible explicitly. In order to eliminate these drawbacks, we point out in this section, two approaches which lead to a practical and more simple solution of the (ELP) equation.


#### 2.2.1. The Method of Evolution Equation [19]

Denote by(18)(γ):v=v0(s),0≤s≤1,
an initial estimate of the optimal snake and let consider a family of curves (contours)(19)(γt):v=v(t,s),v∈C2(R+×[0,1],R2),
where the parameter t≥0 describes the evolution in time of the snake and s∈[0,1] is the standard parameter of the curve. The evolution equation associated to the dynamic model is(20)∂v∂t+∂2∂s2(w2(s)∂2v∂s2)-∂∂s(w1(s)∂v∂s)-J2(∇k)∂v∂s+∇P(v)-(∇k)(J2v′)=0,
together with the initial condition:(21)v(0,s)=v0(s),0≤s≤1,
and the boundary conditions(22)v(t,0)=v0(0),v(t,1)=v0(1),∂v∂s(t,0)=v0′(0),∂v∂s(t,1)=v0′(1),t≥0.
A solution of the static problem described by (ELP) equation (11) is achieved when the solution v(t,s) becomes stable with respect to the time parameter, that is, limt→∞(∂v/∂t)(t,s)=0, uniformly with respect to the parameter s∈[0,1]; in this case, the evolution equation (20) provides a solution of the static problem (11). According to [20], we note that this approach of making the time derivative term vanish is equivalent to applying a gradient descent algorithm to find the local minimum of the energy functional E(v).


#### 2.2.2. The Method of the Lagrange Dynamics [3, 18]

A dynamic snake is represented by introducing a time-varying contour(23)v(t,s)=(x(t,s),y(t,s))T,
see (19), with a mass density μ(s) and a damping density ω(s). The Lagrange equation for a snake defined in Section 2.1 is(24)μ∂2v∂t2+ω∂v∂t+∂2∂s2(w2(s)∂2v∂s2)-∂∂s(w1(s)∂v∂s)+∇P(v)-J2(∇k)∂v∂s-(∇k)(J2v′)=0.The first two terms in the left hand side of (24) represent the inertial and damping forces, while the remaining terms, see also (11), represent the internal stretching force (the term containing ∂v/∂s), the bending (rigidity) force (the term containing ∂2v/∂s2), the external force (∇P(v)) and the balloon-type force (the last two terms). Equilibrium is achieved when these forces balance and the contour comes to rest, that is,(25)∂v∂t=∂2v∂t2=0,
which leads to the equilibrium condition (11).


### 2.3. Deformable Surfaces (3D Models)

In this section we define briefly the notion of deformable 3D model (deformable surface), both in the static and dynamic forms, and we describe a method for reducing the problem of its optimization to a 2D modelling problem.


#### 2.3.1. Energy-Minimizing Surfaces

From mathematical point of view, a 3D variational deformable model is emphasized by a family𝒜 of parameterized smooth surfaces with given boundary conditions, named admissible surfaces, and an associated energy functional.Denoting byD=[0,1]×[0,1] the unit square of ℝ2, let us consider a surface of vectorial equation:(26)(S):v=v(s,r),(s,r)∈D,
where v∈C2(D,ℝ3), v=(x,y,z)T; in this subsection we set |v|2=x2+y2+z2, vs=∂v/∂s, vss=∂2v/∂s2, vsr=∂2v/∂s∂r, vrr=∂2v/∂r2. Given the functions g∈C2(∂D,ℝ3) and h∈C1(∂D,ℝ3), where ∂D is the boundary of D, let 𝒜 be the set of admissible deformations, which consists of all functions v∈C2(D,ℝ3) satisfying the boundary conditions v(s,r)=g(s,r) and (∂v/∂n)(s,r)=h(s,r) on ∂D, where n is the normal vector with respect to the surface (S) defined by (26). Further, let us consider the following functions: the image intensity functionI∈C2(ℝ3); the potential function associated to the external forces P(v)=-λ|∇I(v)|2, λ>0; the control functions corresponding to the internal forces acting on the shape of the surface, namely, the elasticity functions w10(s;r) and w01(s;r); the rigidity functionsw20(s;r) and w02(s;r), and the twist resistance functionw11(s;r). The energy functional E:𝒜→ℝ, associated to these data, is defined as follows:(27)E(v)=∬DF(v,vs,vr,vss,vsr,vrr)dsdr,
where(28)F(v,vs,vr,vss,vsr,vrr)=w10|vs|2+w01|vr|2+w20|vss|2+2w11|vsr|2+w02|vrr|2+f(v,vs,vr),f(v,vs,vr)=P(v)+det⁡(c0v,vs,vr).
We notice that E(v) represents the sum of the internal energy (the terms of (27) excepting f(v,vs,vr)), the external energy (defined by the term containing P(v)) and the balloon energy, which is added, optionally, by the users (the term including det(c0v,vs,vr)).The triple(𝒜,I,E) is said to be a 3D deformable model, sometimes a deformable surface. The basic problem of the deformable model is to minimize its energy functional, namely, to obtain the optimal deformable surface. To this purpose, the Euler-Gauss-Ostrogradski (EGO) equation of Calculus of Variations, that is,
(29)∂F∂v-∂∂s(∂F∂vs)-∂∂r(∂F∂vr)+∂2∂s2(∂F∂vss)+∂2∂v∂r(∂F∂vsr)+∂2∂r2(∂F∂vrr)=0
is used.By simple calculation, we obtain from (28) and (29):(30)∂2∂s2(w20vss)+∂2∂r2(w02vrr)+2∂2∂s∂r(w11vsr)-∂∂s(w10vs)-∂∂r(w01vr)+12(∇f-∂∂s(∂f∂vs)-∂∂r(∂f∂vr))=0.


#### 2.3.2. Deformable Dynamic 3D Models

Similarly to the 2D model, we can suppose that a rough prior estimate of surface is accessible, namely,(31)(S0):v=v0(s,r),(s,r)∈D.
Further, this surface is refined step by step, according to (EGO) equation; so, a sequence of surfaces, which leads to the energy-minimizing surface, is provided. More exactly, let(32)(St):v=v(t,s,r),t≥0,(s,r)∈D,
be a family of surfaces, where the parameter t describes the evolution in time of the model. We associate to the previous static model (𝒜,I,E) the evolution equation(33)∂v∂t+G(v,vs,vr,vss,vsr,vrr)=0,
where G(v,vs,vr,vss,vsr,vrr) is the left hand member of (30), together with theinitial estimate (condition)(34)v(0,s,r)=v0(s,r),(s,r)∈D,
and the boundary dynamic conditions(35)v(t,s,r)=v0(s,r),(s,r)∈∂D,t≥0,∂v(t,s,r)∂n=∂v0(s,r)∂n,(s,r)∈∂D,t≥0.
A solution of the “static” problem described by (30) is achieved, when the solution v(t,s,r) becomes stable with respect to the time parameter, that is, limt→∞(∂v/∂t)(t,s,r)=0, uniformly, with respect to (s,r)∈D; in this case, the evolution equation (33) provides a solution of the static problem (30).


#### 2.3.3. The Simplified 2D Model

The problem of finding directly energy-minimizing surfaces, that is, solutions of the p.d.e. (30), is not practically possible because these solutions contain long and complicated expressions or their explicit form is inaccessible. On the other hand, by using discretized schemes for solving (33), we get a system of algebraic equations with a high computational level. These drawbacks are eliminated by passing to a 2D modeling problem, [19]. More exactly, the third component z of (S) is constrained to depend only on r, by setting z(s,r)=r. So, the surface that we seek is given as a sequence of plane curves, named slices, and the parameter r of (26) becomes the index of the corresponding slice. In this approach, the surface that we seek is viewed as a sequence of the planar curves (slices), indexed by the parameter r, so that each fixed value of r provides a closed curve, lying in a slice of the 3D-image. Consequently, let(36)(γr):v(s)=(x(s),y(s)),s∈[0,1],
be the 2D curve obtained by applying this reconstruction method, for a given r.Under the hypothesis thatwij are positive constants, the (EGO) equation (29), which corresponds to (γr), is(37)2w20d4vds4-2w10d2vds2-c0J2dvds+∇P=0,
where  J2=[01-10].If we consider in (37) c0=-0.02,w10=2.5, w20=0.4, P=r(x2+y2) and r=0.1,0.2,…,1 with boundary conditions x(0)=x(1)=1+(r2(1-r)2)/25, x′(0)=x′(1)=r(1-r)/20, y(0)=y(1)=0+r2(1-r)/25, y′(0)=y′(1)=2r(1-r)/5, we obtain the graphs of the slices and a 3D reconstruction of the surface, as we can see in Figures 3(a) and 3(b).(a)(b)In what follows we shall restrict to the study of 2D deformable models.


## 2.1. Energy-Minimizing Snakes (2D Models)

From mathematical point of view, a 2D parametric deformable model (usually known assnake) is provided by a family 𝒜 of parametrized smooth curves satisfying given boundary conditions and an associated energy-functional. More exactly, denote by C2([0,1],ℝ2) the space of all vectorial functions v=(x,y)T so that the scalar functions x=x(s) and y=y(s), 0≤s≤1 are continuous together with their derivatives up to the second-order on the standard interval [0,1], that is, x,y∈C2[0,1]; obviously, we can consider an arbitrary compact interval [a,b] of the real axis instead of [0,1]. The family 𝒜 of admissible deformations consists of all parametrized curves (snakes):(1)(γ):v(s)=(x(s),y(s))T,0≤s≤1,v∈C2([0,1],R2),
such that the values v(0),v(1),v′(0), and v′(1) are given; we adopt the notation |v|2=|x|2+|y|2.In order to find the optimal position of the snake, it is necessary to characterize its state, by means of anenergy-functional,that is associated to the class 𝒜. Let us consider the following data:(i)
theweight-functionsw1(s) and w2(s), which control the elasticity and the rigidity of (γ), respectively; generally, these are nonnegative scalar functions of class C2[0,1],(ii)
theimage intensity function I=I(x,y), which is a real function of class C2(ℝ2),(iii)
thepotential associated to the external forces, represented by a real function P(v)=P(x,y), of class C2(ℝ2). The simplest useful choice for the potential is P(v)=w3I(v), where w3 is a weight-scalar. The most used choice is P(v)=-λ|∇I(v)|, where λ>0 is a given scalar and ∇=(∂/∂x,∂/∂y)T is the Hamilton (nabla) operator; this choice will be used in this paper, too. Note that P can be defined also by P=Gσ0*I, that is, the Gaussian (variance σ0) filtered image of the input image I [10], and(iv)
the vectorial functionk(v)=(k1(v),k2(v))T of class C1(ℝ2,ℝ2) which control the local dilatation or local contraction of (γ) along its normal; usually, we take k(v)=cv, with c∈ℝ.The shape of the snake(γ) subject to the image I(v) is dictated by the energy functional:(2)E(v)=Eint⁡(v)+Eext(v)+Ebal(v),
where the terms of the right hand of (2) are defined as follows.Theinternal energy(3)Eint⁡(v)=Eels(v)+Erig(v)
is obtained by adding the elastic energy(4)Eels(v)=∫01α(s)|v′(s)|2ds,
and the rigid (bending) energy(5)Erig(v)=∫01β(s)|v′′(s)|2ds.The internal energy characterizes the deformation of a stretchy, flexible snake (contour). The values ofw1(s) and w2(s) show the extent to which the snake can stretch or bend at an arbitrary point (x(s),y(s)) of the snake.Theexternal energy,derived from the image, is given by(6)Eext(v)=∫01P(v(s))ds=-λ∫01|∇I(x,y)|2ds,
and it allows to find the edges in an image so that the snake is attracted to contour with large image gradients.Theballoon energyis an energy of constrained-type, defined as(7)Ebal(v)=-∫01det⁡(k(v),v′)ds=-∫01(k1y′-k2x′)ds.
This energy can be added, optionally, by users, in order to expand (or contract) the snake.Denoting by(8)F(s,v,v′,v′′)=w1(s)|v′2(s)|+w2(s)|v′′2(s)|+P(v(s))-det⁡(k(v),v′),P(v)=-λ|∇I(v)|2,
the following expression of energy-functional E(v) is obtained from (2)–(8):(9)E(v)=∫01F(s,v,v′,v′′)ds.By definition, the triple(𝒜,I,E) is said to be a deformable 2D model (snake).The basic goal of a deformable parametric model is to minimize its energy-functionalE(v), which leads to the energy-minimizing snake. The minimization of the snake energy gives rise to the following vectorial Euler-Lagrange-Poisson (ELP) Equation of Calculus of Variations:(10)∂F∂v-dds(∂F∂v′)+d2ds2(∂F∂v′′)=0.
Now, taking into account the relations (8) and (10), we obtain the vectorial (ELP) equation:(11)2w2(s)viv(s)+4w2′(s)v′′′(s)+2[w2′′(s)-w1(s)]v′′(s)-[2w1′(s)I2+Tr⁡(∇k)J2]v′(s)+∇P(v(s))=0,
where I2=(1001), J2=(01-10), and Tr(A) is the trace of a square matrix A.The scalar (ELP) equations, derived from (11), have the form:(12)2w2xiv+4w2′x′′′+2(w2′′-w1)x′′-2w1′x′-(∂k1∂x+∂k2∂y)y′+∂P∂x=0,2w2yiv+4w2′y′′′+2(w2′′-w1)y′′-2w1′y′+(∂k1∂x+∂k2∂y)x′+∂P∂y=0.
On the other hand, we infer from(8)(13)F(s,v,v′,v′′)=w1(x′2+y′2)+w2(x′′2+y′′2)+P(x,y)-k1y′+k2x′,
which leads to(14)∂2F∂(v′′)2=(2w1,2w2)T.
According to the Legendre conditions and the hypothesis w2>0, the relation (14) proves that any solution of the (ELP) equations (11) or (12) provides a minimum for the energy-functional E(v), namely, an energy-minimizing snake.Example 1.
If we choose in (12) w1=1, w2=0.05, and the boundary conditions v(0)=v(1)=(0,5)T, v′(0)=v′(1)=(0.5,0.5)T we obtain the general solution of the (ELP) equation:
(15)x(s)=C1e3.8042s+C2e-3.8042s+C3e2.3511s+C4e-2.3511s,y(s)=C5e3.8042s+C6e-3.8042s+C7e2.3511s+C8e-2.3511s.
Using boundary conditions, we obtain the graph of the curve (γ) in Figure 1.Figure 1Example 2.
If we choose in (12) w1=1, w2=1, k(v)=10v, I(v)=I(x,y)=(22/3)(x3/2+y3/2), λ=1, and the boundary conditions v(0)=(1,1)T, v′(0)=(3.1,-8.1)T, v(π)=(1-π/10,-2-π/10-3cosh(2π))T, v′(π)=(2.1-cosh(2π),-8.1-2sinh(2π)), then the (ELP) equations have the form:
(16)xiv-x′′-10y′-1=0,yiv-y′′+10x′-1=0,
with the analytical solutions:
(17)x(s)=sin⁡⁡2s+cos⁡⁡2s+sin⁡⁡scosh⁡⁡2s+0.1s,y(s)=-4sin⁡⁡2s+4cos⁡⁡2s+3cos⁡⁡scosh⁡⁡2s,-4sinh⁡⁡2ssin⁡⁡s-0.1s-6.
The graph of the curve (γ) is given in Figure 2.Figure 2


## 2.2. Deformable Dynamic 2D Models

Roughly speaking, the differential fourth-order vectorial equation (11) or the differential eight-order system (12) may have many solutions, which leads to many possible energy-minimizing snakes. As we have seen in the preceding examples (Section 2.1), these solutions have a complicated form; moreover, they are often inaccessible explicitly. In order to eliminate these drawbacks, we point out in this section, two approaches which lead to a practical and more simple solution of the (ELP) equation.


### 2.2.1. The Method of Evolution Equation [19]

Denote by(18)(γ):v=v0(s),0≤s≤1,
an initial estimate of the optimal snake and let consider a family of curves (contours)(19)(γt):v=v(t,s),v∈C2(R+×[0,1],R2),
where the parameter t≥0 describes the evolution in time of the snake and s∈[0,1] is the standard parameter of the curve. The evolution equation associated to the dynamic model is(20)∂v∂t+∂2∂s2(w2(s)∂2v∂s2)-∂∂s(w1(s)∂v∂s)-J2(∇k)∂v∂s+∇P(v)-(∇k)(J2v′)=0,
together with the initial condition:(21)v(0,s)=v0(s),0≤s≤1,
and the boundary conditions(22)v(t,0)=v0(0),v(t,1)=v0(1),∂v∂s(t,0)=v0′(0),∂v∂s(t,1)=v0′(1),t≥0.
A solution of the static problem described by (ELP) equation (11) is achieved when the solution v(t,s) becomes stable with respect to the time parameter, that is, limt→∞(∂v/∂t)(t,s)=0, uniformly with respect to the parameter s∈[0,1]; in this case, the evolution equation (20) provides a solution of the static problem (11). According to [20], we note that this approach of making the time derivative term vanish is equivalent to applying a gradient descent algorithm to find the local minimum of the energy functional E(v).


### 2.2.2. The Method of the Lagrange Dynamics [3, 18]

A dynamic snake is represented by introducing a time-varying contour(23)v(t,s)=(x(t,s),y(t,s))T,
see (19), with a mass density μ(s) and a damping density ω(s). The Lagrange equation for a snake defined in Section 2.1 is(24)μ∂2v∂t2+ω∂v∂t+∂2∂s2(w2(s)∂2v∂s2)-∂∂s(w1(s)∂v∂s)+∇P(v)-J2(∇k)∂v∂s-(∇k)(J2v′)=0.The first two terms in the left hand side of (24) represent the inertial and damping forces, while the remaining terms, see also (11), represent the internal stretching force (the term containing ∂v/∂s), the bending (rigidity) force (the term containing ∂2v/∂s2), the external force (∇P(v)) and the balloon-type force (the last two terms). Equilibrium is achieved when these forces balance and the contour comes to rest, that is,(25)∂v∂t=∂2v∂t2=0,
which leads to the equilibrium condition (11).


## 2.2.1. The Method of Evolution Equation [19]

Denote by(18)(γ):v=v0(s),0≤s≤1,
an initial estimate of the optimal snake and let consider a family of curves (contours)(19)(γt):v=v(t,s),v∈C2(R+×[0,1],R2),
where the parameter t≥0 describes the evolution in time of the snake and s∈[0,1] is the standard parameter of the curve. The evolution equation associated to the dynamic model is(20)∂v∂t+∂2∂s2(w2(s)∂2v∂s2)-∂∂s(w1(s)∂v∂s)-J2(∇k)∂v∂s+∇P(v)-(∇k)(J2v′)=0,
together with the initial condition:(21)v(0,s)=v0(s),0≤s≤1,
and the boundary conditions(22)v(t,0)=v0(0),v(t,1)=v0(1),∂v∂s(t,0)=v0′(0),∂v∂s(t,1)=v0′(1),t≥0.
A solution of the static problem described by (ELP) equation (11) is achieved when the solution v(t,s) becomes stable with respect to the time parameter, that is, limt→∞(∂v/∂t)(t,s)=0, uniformly with respect to the parameter s∈[0,1]; in this case, the evolution equation (20) provides a solution of the static problem (11). According to [20], we note that this approach of making the time derivative term vanish is equivalent to applying a gradient descent algorithm to find the local minimum of the energy functional E(v).


## 2.2.2. The Method of the Lagrange Dynamics [3, 18]

A dynamic snake is represented by introducing a time-varying contour(23)v(t,s)=(x(t,s),y(t,s))T,
see (19), with a mass density μ(s) and a damping density ω(s). The Lagrange equation for a snake defined in Section 2.1 is(24)μ∂2v∂t2+ω∂v∂t+∂2∂s2(w2(s)∂2v∂s2)-∂∂s(w1(s)∂v∂s)+∇P(v)-J2(∇k)∂v∂s-(∇k)(J2v′)=0.The first two terms in the left hand side of (24) represent the inertial and damping forces, while the remaining terms, see also (11), represent the internal stretching force (the term containing ∂v/∂s), the bending (rigidity) force (the term containing ∂2v/∂s2), the external force (∇P(v)) and the balloon-type force (the last two terms). Equilibrium is achieved when these forces balance and the contour comes to rest, that is,(25)∂v∂t=∂2v∂t2=0,
which leads to the equilibrium condition (11).


## 2.3. Deformable Surfaces (3D Models)

In this section we define briefly the notion of deformable 3D model (deformable surface), both in the static and dynamic forms, and we describe a method for reducing the problem of its optimization to a 2D modelling problem.


### 2.3.1. Energy-Minimizing Surfaces

From mathematical point of view, a 3D variational deformable model is emphasized by a family𝒜 of parameterized smooth surfaces with given boundary conditions, named admissible surfaces, and an associated energy functional.Denoting byD=[0,1]×[0,1] the unit square of ℝ2, let us consider a surface of vectorial equation:(26)(S):v=v(s,r),(s,r)∈D,
where v∈C2(D,ℝ3), v=(x,y,z)T; in this subsection we set |v|2=x2+y2+z2, vs=∂v/∂s, vss=∂2v/∂s2, vsr=∂2v/∂s∂r, vrr=∂2v/∂r2. Given the functions g∈C2(∂D,ℝ3) and h∈C1(∂D,ℝ3), where ∂D is the boundary of D, let 𝒜 be the set of admissible deformations, which consists of all functions v∈C2(D,ℝ3) satisfying the boundary conditions v(s,r)=g(s,r) and (∂v/∂n)(s,r)=h(s,r) on ∂D, where n is the normal vector with respect to the surface (S) defined by (26). Further, let us consider the following functions: the image intensity functionI∈C2(ℝ3); the potential function associated to the external forces P(v)=-λ|∇I(v)|2, λ>0; the control functions corresponding to the internal forces acting on the shape of the surface, namely, the elasticity functions w10(s;r) and w01(s;r); the rigidity functionsw20(s;r) and w02(s;r), and the twist resistance functionw11(s;r). The energy functional E:𝒜→ℝ, associated to these data, is defined as follows:(27)E(v)=∬DF(v,vs,vr,vss,vsr,vrr)dsdr,
where(28)F(v,vs,vr,vss,vsr,vrr)=w10|vs|2+w01|vr|2+w20|vss|2+2w11|vsr|2+w02|vrr|2+f(v,vs,vr),f(v,vs,vr)=P(v)+det⁡(c0v,vs,vr).
We notice that E(v) represents the sum of the internal energy (the terms of (27) excepting f(v,vs,vr)), the external energy (defined by the term containing P(v)) and the balloon energy, which is added, optionally, by the users (the term including det(c0v,vs,vr)).The triple(𝒜,I,E) is said to be a 3D deformable model, sometimes a deformable surface. The basic problem of the deformable model is to minimize its energy functional, namely, to obtain the optimal deformable surface. To this purpose, the Euler-Gauss-Ostrogradski (EGO) equation of Calculus of Variations, that is,
(29)∂F∂v-∂∂s(∂F∂vs)-∂∂r(∂F∂vr)+∂2∂s2(∂F∂vss)+∂2∂v∂r(∂F∂vsr)+∂2∂r2(∂F∂vrr)=0
is used.By simple calculation, we obtain from (28) and (29):(30)∂2∂s2(w20vss)+∂2∂r2(w02vrr)+2∂2∂s∂r(w11vsr)-∂∂s(w10vs)-∂∂r(w01vr)+12(∇f-∂∂s(∂f∂vs)-∂∂r(∂f∂vr))=0.


### 2.3.2. Deformable Dynamic 3D Models

Similarly to the 2D model, we can suppose that a rough prior estimate of surface is accessible, namely,(31)(S0):v=v0(s,r),(s,r)∈D.
Further, this surface is refined step by step, according to (EGO) equation; so, a sequence of surfaces, which leads to the energy-minimizing surface, is provided. More exactly, let(32)(St):v=v(t,s,r),t≥0,(s,r)∈D,
be a family of surfaces, where the parameter t describes the evolution in time of the model. We associate to the previous static model (𝒜,I,E) the evolution equation(33)∂v∂t+G(v,vs,vr,vss,vsr,vrr)=0,
where G(v,vs,vr,vss,vsr,vrr) is the left hand member of (30), together with theinitial estimate (condition)(34)v(0,s,r)=v0(s,r),(s,r)∈D,
and the boundary dynamic conditions(35)v(t,s,r)=v0(s,r),(s,r)∈∂D,t≥0,∂v(t,s,r)∂n=∂v0(s,r)∂n,(s,r)∈∂D,t≥0.
A solution of the “static” problem described by (30) is achieved, when the solution v(t,s,r) becomes stable with respect to the time parameter, that is, limt→∞(∂v/∂t)(t,s,r)=0, uniformly, with respect to (s,r)∈D; in this case, the evolution equation (33) provides a solution of the static problem (30).


### 2.3.3. The Simplified 2D Model

The problem of finding directly energy-minimizing surfaces, that is, solutions of the p.d.e. (30), is not practically possible because these solutions contain long and complicated expressions or their explicit form is inaccessible. On the other hand, by using discretized schemes for solving (33), we get a system of algebraic equations with a high computational level. These drawbacks are eliminated by passing to a 2D modeling problem, [19]. More exactly, the third component z of (S) is constrained to depend only on r, by setting z(s,r)=r. So, the surface that we seek is given as a sequence of plane curves, named slices, and the parameter r of (26) becomes the index of the corresponding slice. In this approach, the surface that we seek is viewed as a sequence of the planar curves (slices), indexed by the parameter r, so that each fixed value of r provides a closed curve, lying in a slice of the 3D-image. Consequently, let(36)(γr):v(s)=(x(s),y(s)),s∈[0,1],
be the 2D curve obtained by applying this reconstruction method, for a given r.Under the hypothesis thatwij are positive constants, the (EGO) equation (29), which corresponds to (γr), is(37)2w20d4vds4-2w10d2vds2-c0J2dvds+∇P=0,
where  J2=[01-10].If we consider in (37) c0=-0.02,w10=2.5, w20=0.4, P=r(x2+y2) and r=0.1,0.2,…,1 with boundary conditions x(0)=x(1)=1+(r2(1-r)2)/25, x′(0)=x′(1)=r(1-r)/20, y(0)=y(1)=0+r2(1-r)/25, y′(0)=y′(1)=2r(1-r)/5, we obtain the graphs of the slices and a 3D reconstruction of the surface, as we can see in Figures 3(a) and 3(b).(a)(b)In what follows we shall restrict to the study of 2D deformable models.


## 2.3.1. Energy-Minimizing Surfaces

From mathematical point of view, a 3D variational deformable model is emphasized by a family𝒜 of parameterized smooth surfaces with given boundary conditions, named admissible surfaces, and an associated energy functional.Denoting byD=[0,1]×[0,1] the unit square of ℝ2, let us consider a surface of vectorial equation:(26)(S):v=v(s,r),(s,r)∈D,
where v∈C2(D,ℝ3), v=(x,y,z)T; in this subsection we set |v|2=x2+y2+z2, vs=∂v/∂s, vss=∂2v/∂s2, vsr=∂2v/∂s∂r, vrr=∂2v/∂r2. Given the functions g∈C2(∂D,ℝ3) and h∈C1(∂D,ℝ3), where ∂D is the boundary of D, let 𝒜 be the set of admissible deformations, which consists of all functions v∈C2(D,ℝ3) satisfying the boundary conditions v(s,r)=g(s,r) and (∂v/∂n)(s,r)=h(s,r) on ∂D, where n is the normal vector with respect to the surface (S) defined by (26). Further, let us consider the following functions: the image intensity functionI∈C2(ℝ3); the potential function associated to the external forces P(v)=-λ|∇I(v)|2, λ>0; the control functions corresponding to the internal forces acting on the shape of the surface, namely, the elasticity functions w10(s;r) and w01(s;r); the rigidity functionsw20(s;r) and w02(s;r), and the twist resistance functionw11(s;r). The energy functional E:𝒜→ℝ, associated to these data, is defined as follows:(27)E(v)=∬DF(v,vs,vr,vss,vsr,vrr)dsdr,
where(28)F(v,vs,vr,vss,vsr,vrr)=w10|vs|2+w01|vr|2+w20|vss|2+2w11|vsr|2+w02|vrr|2+f(v,vs,vr),f(v,vs,vr)=P(v)+det⁡(c0v,vs,vr).
We notice that E(v) represents the sum of the internal energy (the terms of (27) excepting f(v,vs,vr)), the external energy (defined by the term containing P(v)) and the balloon energy, which is added, optionally, by the users (the term including det(c0v,vs,vr)).The triple(𝒜,I,E) is said to be a 3D deformable model, sometimes a deformable surface. The basic problem of the deformable model is to minimize its energy functional, namely, to obtain the optimal deformable surface. To this purpose, the Euler-Gauss-Ostrogradski (EGO) equation of Calculus of Variations, that is,
(29)∂F∂v-∂∂s(∂F∂vs)-∂∂r(∂F∂vr)+∂2∂s2(∂F∂vss)+∂2∂v∂r(∂F∂vsr)+∂2∂r2(∂F∂vrr)=0
is used.By simple calculation, we obtain from (28) and (29):(30)∂2∂s2(w20vss)+∂2∂r2(w02vrr)+2∂2∂s∂r(w11vsr)-∂∂s(w10vs)-∂∂r(w01vr)+12(∇f-∂∂s(∂f∂vs)-∂∂r(∂f∂vr))=0.


## 2.3.2. Deformable Dynamic 3D Models

Similarly to the 2D model, we can suppose that a rough prior estimate of surface is accessible, namely,(31)(S0):v=v0(s,r),(s,r)∈D.
Further, this surface is refined step by step, according to (EGO) equation; so, a sequence of surfaces, which leads to the energy-minimizing surface, is provided. More exactly, let(32)(St):v=v(t,s,r),t≥0,(s,r)∈D,
be a family of surfaces, where the parameter t describes the evolution in time of the model. We associate to the previous static model (𝒜,I,E) the evolution equation(33)∂v∂t+G(v,vs,vr,vss,vsr,vrr)=0,
where G(v,vs,vr,vss,vsr,vrr) is the left hand member of (30), together with theinitial estimate (condition)(34)v(0,s,r)=v0(s,r),(s,r)∈D,
and the boundary dynamic conditions(35)v(t,s,r)=v0(s,r),(s,r)∈∂D,t≥0,∂v(t,s,r)∂n=∂v0(s,r)∂n,(s,r)∈∂D,t≥0.
A solution of the “static” problem described by (30) is achieved, when the solution v(t,s,r) becomes stable with respect to the time parameter, that is, limt→∞(∂v/∂t)(t,s,r)=0, uniformly, with respect to (s,r)∈D; in this case, the evolution equation (33) provides a solution of the static problem (30).


## 2.3.3. The Simplified 2D Model

The problem of finding directly energy-minimizing surfaces, that is, solutions of the p.d.e. (30), is not practically possible because these solutions contain long and complicated expressions or their explicit form is inaccessible. On the other hand, by using discretized schemes for solving (33), we get a system of algebraic equations with a high computational level. These drawbacks are eliminated by passing to a 2D modeling problem, [19]. More exactly, the third component z of (S) is constrained to depend only on r, by setting z(s,r)=r. So, the surface that we seek is given as a sequence of plane curves, named slices, and the parameter r of (26) becomes the index of the corresponding slice. In this approach, the surface that we seek is viewed as a sequence of the planar curves (slices), indexed by the parameter r, so that each fixed value of r provides a closed curve, lying in a slice of the 3D-image. Consequently, let(36)(γr):v(s)=(x(s),y(s)),s∈[0,1],
be the 2D curve obtained by applying this reconstruction method, for a given r.Under the hypothesis thatwij are positive constants, the (EGO) equation (29), which corresponds to (γr), is(37)2w20d4vds4-2w10d2vds2-c0J2dvds+∇P=0,
where  J2=[01-10].If we consider in (37) c0=-0.02,w10=2.5, w20=0.4, P=r(x2+y2) and r=0.1,0.2,…,1 with boundary conditions x(0)=x(1)=1+(r2(1-r)2)/25, x′(0)=x′(1)=r(1-r)/20, y(0)=y(1)=0+r2(1-r)/25, y′(0)=y′(1)=2r(1-r)/5, we obtain the graphs of the slices and a 3D reconstruction of the surface, as we can see in Figures 3(a) and 3(b).(a)(b)In what follows we shall restrict to the study of 2D deformable models.


## 3. An ELP-Algorithm for Obtaining Energy-Minimizing Snakes

In this section we suppose that the following hypotheses are satisfied: the control functionsw1 and w2 are positive constants, the curves of the family (γt) given by (18) and (19) are closed for every t≥0 and k(v)=c0v, c0∈ℝ+.  Thus, the (ELP) evolution equation (20) becomes(38)2∂v∂t+2w2∂4v∂s4-2w1∂2v∂s2-2c0J2∂v∂s+∇P=0,v=v(t,s).
In order to solve numerically the partial differential equation (38), we focus on the method of finite differences, which is widely used in image processing [21]. Let δ and h be the time and the space discretization steps, respectively, and denote by ℛ={(tk,si),k≥0,0≤i≤N} the plane net of discretization, with N∈ℕ*, Nh=1, tk=kδ, and si=ih. The following notations will be used, too: vik=v(tk,si), vk=(xk,yk)T, k≥0;gk=(g1k,g2k)T, with g1k=(-1/2)((∂P/∂x)(vk)), g2k=(-1/2)((∂P/∂y)(vk)); obviously, vik=vi+Nk, i∈ℤ, because (γt),t≥0, is a closed curve. Also, we set(39)α=w1h2,β=w2h4,γ=c0h.


### 3.1. Explicit Finite Difference Scheme

We approximate the partial derivatives involved in the (ELP) evolution equation (38) as follows:(40)∂v∂t(tk,si)≈1δ(vik+1-vik);∂v∂s(tk,si)≈1h(vi+1k-vik),∂2v∂s2(tk,si)≈1h2(vi+1k-2vik+vi-1k),∂4v∂s4(tk,si)≈1h4(vi+2k-4vi+1k+6vik-4vi-1k+vi-2k).
By replacing the relations (40) in the partial differential equation (38), it result a system of algebraic equations; denoting by Vk=(Xk,Yk)T the solutions of this system (which approximate the exact values vik of (38) at the nodes of ℛ), we get the vectorial formula:(41)Vik+1-Vikδ+β(Vi+2k-4Vi+1k+6Vik-4Vi-1k+Vi-2k)-α(Vi+1k-2Vik+Vi-1k)-γJ2(Vi+1k-Vik)+12∇P=0,0≤i≤N;k≥0,
where(42)Vik=(Xik,Yik)
and α,β,γ are given by (39). The scalar equations corresponding to (41) are the following:(43)Xik+1-Xikδ+β(Xi+2k-4Xi+1k+6Xik-4Xi-1k+Xi-2k)-α(Xi+1k-2Xik+Xi-1k)-γ(Yi+1k-Yik)+12∂P∂x(Xik,Yik)=0,0≤i≤N-1;k≥0.Yik+1-Yikδ+β(Yi+2k-4Yi+1k+6Yik-4Yi-1k+Yi-2k)-α(Yi+1k-2Yik+Yi-1k)+γ(Xi+1k-Xik)+12∂P∂y(Xik,Yik)=0,0≤i≤N-1;k≥0.
Now, let K be the stiffness matrix associated to the explicit finite difference scheme, defined as the circular matrix of order N, whose first row is(44)(a1,a2,a3,0,…,0,a3,a2),
where(45)a1=2α+6β,a2=-α-4β,a3=β.
Denote by L the circular (square) matrix of order N defined by the first row (1,-1,0,0,…,0) and let IN be the identity matrix of order N. The relations (41) and (43) can be written in a matricial form as:(46)Vk+1=(IN-δK)Vk-γδL(J2Vk)+δgk,k≥0,
(47)Xk+1=(IN-δK)Xk-γδLYk+δg1k,Yk+1=(IN-δK)Yk+γδLXk+δg2k,k≥0
respectively.In what follows, the formulas (41)–(47) will be referred as (ELP) algorithm for obtaining an energy minimizing snake (in its approximating form).


### 3.2. The Residue of (ELP) algorithm

Taking into account the relation (41), the residue associated to the (ELP) algorithm is(48)Rvi=vik+1-vikδ+β(vi+2k-4vi+1k+6vik-4vi-1k+vi-2k)-α(vi+1k-2vik+vi-1k)-γJ2(vi+1k-vik)+12∇P(vik),0≤i≤N,k≥0.
By using Taylor expansions at the point (tk,si)∈ℛ we obtain(49)vik+1=vik+δ∂v∂t+δ22!∂2v∂t2+δ33!∂3v∂t3+⋯,vi±1k=vik±h∂v∂s+h22!∂2v∂s2±h33!∂3v∂s3+⋯,vi±2k=vik±2h∂v∂s+(2h)22!∂2v∂s2±(2h)33!∂3v∂s3+⋯,
where the partial derivatives ∂v/∂t and ∂lv/∂sl, l≥1 are computed at the point (tk,si)=(kδ,ih)∈ℛ.By replacing the expansions (49) in the residue’s formula (48) and using the relations (39), we derive(50)Rvi=δ(12∂2v∂t2+δ6∂3v∂t3+⋯)(tk,si)+h2w2(16∂6v∂s6+1275040h2∂8v∂s8+⋯)(tk,si)-w1h2(112∂4v∂s4+h260∂6v∂s6+⋯)(tk,si)-c0J2h(12∂2v∂v2+h224∂4v∂v4+⋯)(tk,si),0≤i≤N-1,k≥0.
If the partial derivatives of the vectorial function v are uniformly bounded on D, the relations (50) give the following estimate concerning the residue of (ELP) algorithm:(51)Rvi={O(δ)+O(h),ifc0>0,O(δ)+O(h2),ifc0=0.
Notice that the condition c0=0 means that there are not existing constrains defined by the users.


### 3.3. The Consistency of the ELP algorithm

LetTr(vi)=δRvi be the truncature errorof (ELP) algorithm at the kth iteration. Under the assumption of uniform boundedness of the partial derivatives of the vectorial function v, it follows from (51): (52)Tr⁡(vi)={O(δ2)+O(δh),ifc0>0O(δ2)+O(δh2),ifc0=0.
The relations (52) characterize the accuracy of the discretized scheme providing the (ELP)-Algorithm.On the other hand, the equality(53)lim⁡δ→0h→0Tr⁡(vi)δ=lim⁡δ→0h→0Rvi=0,
which results from (52), shows that this discretized scheme is consistent.


### 3.4. Approximation Error and the Convergence

Let us consider the approximation-errorεik at the point (tk,si)∈ℛ, namely(54)εik=vik-Vik,0≤i≤N-1,k≥0.
By replacing Vik=vik-εik from (54) into (41) and taking into account the expressions (49) and the definition (48) of Rvi, we get(55)εik+1=δRvi-βδεi+2k+δ((4α+β)I2+γJ2)εi+1k+((1-6βδ-2αδ)I2-γJ2)εik+δ(4β+α)εi-1k-βδεi-2.
Let(56)Ek=max⁡{|εi-2k|,|εi-1k|,|εik|,|εi+1k|,|εi+2k|},k≥0,
be the approximation error of (ELP) algorithm at kth iteration.The relations (55) and (56) yield:(57)Ek+1≤δ|Rvi|+(βδ+δ(4β+α)2+γ2+(1-6βδ-2αδ)2+γ2δ2+4βδ+αδ+βδδ(4β+α)2+γ2)Ek.
On the other hand, it follows from (50):(58)|Rvi|≤M1δ+|2w2-w1|M2h2+c0M3h,0≤i≤N-1,
where Mj, j≥1 are positive constants, which do not depend on δ and h.Now, the relations (57) and (58), combined with the classic inequality(59)x2+y2≤|x|+|y|,
provide the estimate:(60)Ek+1≤(10βδ+2αδ+2γδ+|1-6βδ-2αδ|)Ek+A(h,δ),
with(61)A(h,δ)=M1δ2+M2|2w2-w1|δh2+M3c0δh.
Denote by(62)ε=δh4
and let us assume that the inequality(63)6εw2(k+1)≤1
holds. It is a simple exercise to show that the relation (63) entails the inequality(64)6βδ+2αδ≤1
for N sufficiently large. Now, the relations (60) and (64) lead to:(65)Ek+1≤qEk+A(h,δ),k≥0,
where(66)E0=0,q=1+4βδ+2γδ.
Writing (65) successively for k, k-1,…,1, we get(67)Ek+1≤qk+1-1q-1A(h,δ),k≥0.
Taking into account that γ≤β (for N sufficiently large), the relations (66) and (39) imply 1+4w2ε≤q≤1+6w2ε, so that the relations (67) and (63), combined with the inequality (1+x)1/x≤e, x>0, yield(68)Ek+1≤e-14w2εA(h,δ).
Finally, we derive from (61), (62), and (68):(69)Ek+1={O(δh4)+O(h5),ifc0>0,O(δh4)+O(h6),ifc0=0.
It follows from (69) that Ek+1→0 if h→0; it is easily seen that, according to the relations (62) and (63), the hypothesis h→0 implies δ→0; consequently the following result holds.If the inequality (63) fulfills, then the (ELP) algorithm (46) is convergent and its approximation error at the(k+1)th iteration is given by the relation (69).


### 3.5. The Stability

The intuitive idea regarding the stability is that small errors in the initial conditions of a partial differential equation should cause small errors in its solution. In fact, the study of the stability is useful in connection with the theorem of Lax concerning the convergence of the discretized schemes, [21].The aim of this subsection is to examine the stability of the (ELP) algorithm (46), with c0=0. By omitting the small terms δRvi of (55), we get the relation:(70)εik+1=(1-6βδ-2αδ)εik+(αδ+4βδ)(εi+1k+εi-1k)-βδ(εi+2k+εi-2k),k≥0.
To apply the stability criterion of von Neumann, [22] we set(71)r1=αδ,r2=βδ,η1=ω1h,η2=ω2h,(72)εik=exp⁡(νkl)exp⁡(jωh)=(μkejω1ih,μkejω2ih)T,μ=exp⁡(νl),
where j and ν=ν(ω) are complex numbers, j2=-1 and ω=(ω1,ω2)T denotes the frequency.Now, we obtain from (70), (71), and (72):(73)μ=(1-2r1-6r2)+(r1+4r2)(ejη1+e-jη1)-r2(e2jη1+e-2jη1),μ=(1-2r1-6r2)+(r1+4r2)(ejη2+e-jη2)-r2(e2jη2+e-2jη2).
We choose ω1=ω2 and let η=η1/2=η2/2. The trigonometric formulas ejα+e-jα=2cosα, 1-cosα=2sin2α/2, and 1-cos2α=8sin2α/2cos2α/2, α∈ℝ together with (73) give:(74)μ=1-4r1sin⁡2η-16r2sin⁡4η.
On the other hand, according to the relation |exp(jωit)|=1, it is easy to see that the error εik of (72) does not increase in time if |μ|≤1, so that we infer from (74):(75)2r1+8r2≤1,
which represent precisely the stability criterion of von Neumann for the (ELP) algorithm (46).A combination of the relations (39), (62), (71), and (75) provides the following equivalent form of the stability condition of von Neumann:(76)2ε(4w2+w1h2)≤1.


## 3.1. Explicit Finite Difference Scheme

We approximate the partial derivatives involved in the (ELP) evolution equation (38) as follows:(40)∂v∂t(tk,si)≈1δ(vik+1-vik);∂v∂s(tk,si)≈1h(vi+1k-vik),∂2v∂s2(tk,si)≈1h2(vi+1k-2vik+vi-1k),∂4v∂s4(tk,si)≈1h4(vi+2k-4vi+1k+6vik-4vi-1k+vi-2k).
By replacing the relations (40) in the partial differential equation (38), it result a system of algebraic equations; denoting by Vk=(Xk,Yk)T the solutions of this system (which approximate the exact values vik of (38) at the nodes of ℛ), we get the vectorial formula:(41)Vik+1-Vikδ+β(Vi+2k-4Vi+1k+6Vik-4Vi-1k+Vi-2k)-α(Vi+1k-2Vik+Vi-1k)-γJ2(Vi+1k-Vik)+12∇P=0,0≤i≤N;k≥0,
where(42)Vik=(Xik,Yik)
and α,β,γ are given by (39). The scalar equations corresponding to (41) are the following:(43)Xik+1-Xikδ+β(Xi+2k-4Xi+1k+6Xik-4Xi-1k+Xi-2k)-α(Xi+1k-2Xik+Xi-1k)-γ(Yi+1k-Yik)+12∂P∂x(Xik,Yik)=0,0≤i≤N-1;k≥0.Yik+1-Yikδ+β(Yi+2k-4Yi+1k+6Yik-4Yi-1k+Yi-2k)-α(Yi+1k-2Yik+Yi-1k)+γ(Xi+1k-Xik)+12∂P∂y(Xik,Yik)=0,0≤i≤N-1;k≥0.
Now, let K be the stiffness matrix associated to the explicit finite difference scheme, defined as the circular matrix of order N, whose first row is(44)(a1,a2,a3,0,…,0,a3,a2),
where(45)a1=2α+6β,a2=-α-4β,a3=β.
Denote by L the circular (square) matrix of order N defined by the first row (1,-1,0,0,…,0) and let IN be the identity matrix of order N. The relations (41) and (43) can be written in a matricial form as:(46)Vk+1=(IN-δK)Vk-γδL(J2Vk)+δgk,k≥0,
(47)Xk+1=(IN-δK)Xk-γδLYk+δg1k,Yk+1=(IN-δK)Yk+γδLXk+δg2k,k≥0
respectively.In what follows, the formulas (41)–(47) will be referred as (ELP) algorithm for obtaining an energy minimizing snake (in its approximating form).


## 3.2. The Residue of (ELP) algorithm

Taking into account the relation (41), the residue associated to the (ELP) algorithm is(48)Rvi=vik+1-vikδ+β(vi+2k-4vi+1k+6vik-4vi-1k+vi-2k)-α(vi+1k-2vik+vi-1k)-γJ2(vi+1k-vik)+12∇P(vik),0≤i≤N,k≥0.
By using Taylor expansions at the point (tk,si)∈ℛ we obtain(49)vik+1=vik+δ∂v∂t+δ22!∂2v∂t2+δ33!∂3v∂t3+⋯,vi±1k=vik±h∂v∂s+h22!∂2v∂s2±h33!∂3v∂s3+⋯,vi±2k=vik±2h∂v∂s+(2h)22!∂2v∂s2±(2h)33!∂3v∂s3+⋯,
where the partial derivatives ∂v/∂t and ∂lv/∂sl, l≥1 are computed at the point (tk,si)=(kδ,ih)∈ℛ.By replacing the expansions (49) in the residue’s formula (48) and using the relations (39), we derive(50)Rvi=δ(12∂2v∂t2+δ6∂3v∂t3+⋯)(tk,si)+h2w2(16∂6v∂s6+1275040h2∂8v∂s8+⋯)(tk,si)-w1h2(112∂4v∂s4+h260∂6v∂s6+⋯)(tk,si)-c0J2h(12∂2v∂v2+h224∂4v∂v4+⋯)(tk,si),0≤i≤N-1,k≥0.
If the partial derivatives of the vectorial function v are uniformly bounded on D, the relations (50) give the following estimate concerning the residue of (ELP) algorithm:(51)Rvi={O(δ)+O(h),ifc0>0,O(δ)+O(h2),ifc0=0.
Notice that the condition c0=0 means that there are not existing constrains defined by the users.


## 3.3. The Consistency of the ELP algorithm

LetTr(vi)=δRvi be the truncature errorof (ELP) algorithm at the kth iteration. Under the assumption of uniform boundedness of the partial derivatives of the vectorial function v, it follows from (51): (52)Tr⁡(vi)={O(δ2)+O(δh),ifc0>0O(δ2)+O(δh2),ifc0=0.
The relations (52) characterize the accuracy of the discretized scheme providing the (ELP)-Algorithm.On the other hand, the equality(53)lim⁡δ→0h→0Tr⁡(vi)δ=lim⁡δ→0h→0Rvi=0,
which results from (52), shows that this discretized scheme is consistent.


## 3.4. Approximation Error and the Convergence

Let us consider the approximation-errorεik at the point (tk,si)∈ℛ, namely(54)εik=vik-Vik,0≤i≤N-1,k≥0.
By replacing Vik=vik-εik from (54) into (41) and taking into account the expressions (49) and the definition (48) of Rvi, we get(55)εik+1=δRvi-βδεi+2k+δ((4α+β)I2+γJ2)εi+1k+((1-6βδ-2αδ)I2-γJ2)εik+δ(4β+α)εi-1k-βδεi-2.
Let(56)Ek=max⁡{|εi-2k|,|εi-1k|,|εik|,|εi+1k|,|εi+2k|},k≥0,
be the approximation error of (ELP) algorithm at kth iteration.The relations (55) and (56) yield:(57)Ek+1≤δ|Rvi|+(βδ+δ(4β+α)2+γ2+(1-6βδ-2αδ)2+γ2δ2+4βδ+αδ+βδδ(4β+α)2+γ2)Ek.
On the other hand, it follows from (50):(58)|Rvi|≤M1δ+|2w2-w1|M2h2+c0M3h,0≤i≤N-1,
where Mj, j≥1 are positive constants, which do not depend on δ and h.Now, the relations (57) and (58), combined with the classic inequality(59)x2+y2≤|x|+|y|,
provide the estimate:(60)Ek+1≤(10βδ+2αδ+2γδ+|1-6βδ-2αδ|)Ek+A(h,δ),
with(61)A(h,δ)=M1δ2+M2|2w2-w1|δh2+M3c0δh.
Denote by(62)ε=δh4
and let us assume that the inequality(63)6εw2(k+1)≤1
holds. It is a simple exercise to show that the relation (63) entails the inequality(64)6βδ+2αδ≤1
for N sufficiently large. Now, the relations (60) and (64) lead to:(65)Ek+1≤qEk+A(h,δ),k≥0,
where(66)E0=0,q=1+4βδ+2γδ.
Writing (65) successively for k, k-1,…,1, we get(67)Ek+1≤qk+1-1q-1A(h,δ),k≥0.
Taking into account that γ≤β (for N sufficiently large), the relations (66) and (39) imply 1+4w2ε≤q≤1+6w2ε, so that the relations (67) and (63), combined with the inequality (1+x)1/x≤e, x>0, yield(68)Ek+1≤e-14w2εA(h,δ).
Finally, we derive from (61), (62), and (68):(69)Ek+1={O(δh4)+O(h5),ifc0>0,O(δh4)+O(h6),ifc0=0.
It follows from (69) that Ek+1→0 if h→0; it is easily seen that, according to the relations (62) and (63), the hypothesis h→0 implies δ→0; consequently the following result holds.If the inequality (63) fulfills, then the (ELP) algorithm (46) is convergent and its approximation error at the(k+1)th iteration is given by the relation (69).


## 3.5. The Stability

The intuitive idea regarding the stability is that small errors in the initial conditions of a partial differential equation should cause small errors in its solution. In fact, the study of the stability is useful in connection with the theorem of Lax concerning the convergence of the discretized schemes, [21].The aim of this subsection is to examine the stability of the (ELP) algorithm (46), with c0=0. By omitting the small terms δRvi of (55), we get the relation:(70)εik+1=(1-6βδ-2αδ)εik+(αδ+4βδ)(εi+1k+εi-1k)-βδ(εi+2k+εi-2k),k≥0.
To apply the stability criterion of von Neumann, [22] we set(71)r1=αδ,r2=βδ,η1=ω1h,η2=ω2h,(72)εik=exp⁡(νkl)exp⁡(jωh)=(μkejω1ih,μkejω2ih)T,μ=exp⁡(νl),
where j and ν=ν(ω) are complex numbers, j2=-1 and ω=(ω1,ω2)T denotes the frequency.Now, we obtain from (70), (71), and (72):(73)μ=(1-2r1-6r2)+(r1+4r2)(ejη1+e-jη1)-r2(e2jη1+e-2jη1),μ=(1-2r1-6r2)+(r1+4r2)(ejη2+e-jη2)-r2(e2jη2+e-2jη2).
We choose ω1=ω2 and let η=η1/2=η2/2. The trigonometric formulas ejα+e-jα=2cosα, 1-cosα=2sin2α/2, and 1-cos2α=8sin2α/2cos2α/2, α∈ℝ together with (73) give:(74)μ=1-4r1sin⁡2η-16r2sin⁡4η.
On the other hand, according to the relation |exp(jωit)|=1, it is easy to see that the error εik of (72) does not increase in time if |μ|≤1, so that we infer from (74):(75)2r1+8r2≤1,
which represent precisely the stability criterion of von Neumann for the (ELP) algorithm (46).A combination of the relations (39), (62), (71), and (75) provides the following equivalent form of the stability condition of von Neumann:(76)2ε(4w2+w1h2)≤1.


## 4. Monitoring the Behavior of Prosthetic Surgical Methods and Prosthetic Medical Materials Based on Software Implementation

In order to apply the results of the theoretical researches detailed above in the medical imaging domain, a 3D visual software environment—named MoDef—was implemented, aiming to visualize and follow up the deformation behavior of the surgical (abdominal, maxilla-facial, and orthodontic) prosthetic materials. That is performed on three distinct, but convergent, levels, as follows:(a)
3D reconstruction visual software component, aimed to tracks the evolution of the prosthetic materials, based on processing the US images of the anatomic context of a lot of surgical patients;(b)
deformable prosthetic material’s behavior forecasting software component, based on software tools which implements the above described mathematical methods;(c)
quad comparative parallel tracking software component, aimed to simultaneous supervise in time both (a) and (b) levels, in comparison with the results provided by the stochastic analysis component of the 3D visual software environment MoDef.Concerning the 3D visualizing of the prosthetic meshes by means of the MoDef software environment components, two levels of reconstruction are performed, namely(1)
on the first level, a polynomial interpolation method is applied on each slice of the US image of the prosthetic mesh, acquired based on succeeding positions of the transducer, obtained by rotating them with a constant angle in a same preestablished direction; more exactly, the curves representing the sections of the surgical mesh acquired by the transducer are extracted from the context of the US image, based on specific image processing methods, namely, contour detection methods, that are implemented at the level of the image processing operators of the MoDef environment’s image processing library. Starting with this set of basic mesh surface definition curves, extracted from the US images acquired at pre-established moments in time, a complete and consistent collection of 3D generator curve sets is obtained, by means of 3D polynomial interpolation methods, based on Lagrange, Hermite or Birkhoff operators;(2)
on the second level, the complete collection of the 3D generator curves obtained at the first level is processed based on Blended Interpolating Methods (BIM), as well as with 3D continuous representation techniques, in order to obtain “solid-view,” respectively, “wired-view” representations of the prosthetic mesh.In what follows, some preliminary experiments made in 3DS Max7, followed by some relevant results obtained with the 3D reconstruction component of MoDef 3D Visual environment are presented in Figures4 and 5.Preliminary experiments made in 3DS Max7.
(a)(b)(c)Results obtained with the 3D reconstruction component of MoDef 3D Visual environment: (a) initial 3D representation of the deformable surface of the surgical mesh, (b) curve representing a section of the surgical mesh acquired by the transducer, extracted from the context of the US image based on specific-image processing method, (c) the surface of the prosthetic mesh after the deformations produced in time due to the anatomic assimilation process, and (d) the basic set of generating curves, used to obtain the solid-view representations of the prosthetic mesh.
(a)(b)(c)(d)


## 5. Conclusions

In this paper we considered parametric (variational) deformable models and we developed an iterative method based on finite difference schemes in order to solve numerically the (ELP) equation of Calculus of Variations, which provides the energy minimizing snake. We derived estimates concerning the approximation error related to the corresponding (ELP) algorithm and we established conditions for its convergence and stability. Some considerations about the implementation of the above numerical methods where presented, too. As future targets, we intend to consider probabilistic models which offer an alternative approach by using the Bayes technique, as well as geometric deformable models which provide an efficient alternative to address some limitation of parametric deformable models.



---

*Source: 918510-2012-03-04.xml*