Geometric Computation of Curvature Driven Plane Curve Evolutions.
We present a new numerical scheme for planar curve evolution with a normal velocity equal to $F(\kappa)$, where $\kappa$ is the curvature and F is a nondecreasing function such that F(0)=0 and either $x\mapsto F(x^3)$ is Lipschitz with Lipschitz constant less than or equal to 1 or $F(x)=x^\gamma$ for $\gamma\geq 1/3$. The scheme is completely geometrical and avoids some drawbacks of finite difference schemes. In particular, no special parameterization is needed and the scheme is monotone (that is, if a curve initially surrounds another one, then this remains true during their evolution), which guarantees numerical stability. We prove consistency and convergence of this scheme in a weak sense. Finally, we display some numerical experiments on synthetic and real data.
Introduction
In this paper, we investigate the evolution of a closed smooth plane curve, when each point of
the curve moves with a normal velocity depending on the curvature of the curve at this point.
More precisely, we study evolution of a curve C obeying the equation
@t
where N(s; t) is the inner normal vector to the curve at the point with parameter s at evolution
time t. Equations of type (1) model phenomena in physics or material science. They also play
an important role in digital image analysis. Indeed, it was proved in [1] that any image analysis
process satisfying some reasonable properties and invariance (essentially causality, stability, invariance
with respect to isometries and contrast changes) is described by an equation of type (1),
or more precisely by the corresponding grey level evolution described by the scalar equation
@t
which has to be considered in the viscosity sense [6]. In this equation
x
z e-mail: moisan@cmla.ens-cachan.fr
x Mail adress: CMLA, ENS-CACHAN, 61 Avenue du Pr'esident Wilson, 94 235 Cachan Cedex FRANCE.
Phone:
is the curvature of the level line passing through x. Equation (2) means that the level lines of u
move with respect to Equation (1). The case F (with the convention F
for ff 2 R and x 2 R) is particularly interesting since it is the only contrast invariant equation
that also commutes with affine transformations preserving the area (this set is named the special
affine group, and is composed of the mappings of the type x 7! Ax
is a linear transformation such that det A = 1). The case F
has extensively been studied ([8, 9]), and existence, regularity, and vanishing in finite time have
been proved. The affine invariant case F x 1=3 has also been studied in [2, 22, 23]. For
results have been exposed in [25]. The equivalence between curve motion
and contrast invariant smoothing was proved in [7].
In this paper, we focus on a possible numerical algorithm solving (1). According to [7], if u
is a continuous real-valued function in R 2 , we can solve Equation (2) by applying this algorithm
to all level lines of u. Of course, this implies that these level lines must not intersect during their
evolution. We thus require that, like Equation (1) the algorithm satisfies an inclusion principle,
meaning that the order (with respect to inclusion) is respected. Among the first attempts to
solve Equation (1), Osher and Sethian ([19, 24]) solve (2) by introducing the signed distance
function to the curve at time Unfortunately, this algorithm satisfies neither inclusion
principle nor rotation invariance. In addition, the evolution a priori depends on the chosen
distance function, since the scalar algorithm is not contrast invariant. Moreover, the data (and
thus the CPU time) becomes rapidly huge when a high precision is requested. Some attempts
on curves by finite difference schemes have also been made ([15, 16]) with interesting results,
but the nongeometric nature of the scheme still prevents inclusion principle from being satisfied.
Moreover, for evolution driven by a power a the curvature larger than 1, the discretization of the
curves are enclined to becoming sparser around points with high curvature, thus preventing a
good accuracy. On the contrary, a completely different scheme has been implemented in [17, 18]
for the affine invariant case. It is fully geometrical and also satisfies inclusion principle (implying
numerical stability). In [11], a theoretical algorithm for moving hypersurface by a power of Gauss
Curvature has also been studied. We generalize and implement this algorithm in the plane for
nonconvex curves and for more general functions of the curvature. We just mention that a
numerical approach on curves has the advantage that the resolution is not limited by the pixel
size, which allows a very high precision. Moreover, the computation time is far shorter than
in a scalar approach. On the other hand, topology changes (e.g. a single curve breaking into
two connected components) are automatically handled in a scalar approach. This should not
be a real problem for plane curves evolution since it is likely that no topological changes occur.
This has been proved at least in [9] for the Mean Curvature Flow and in [2] for Affine Curve
evolution. On the contrary, it is known that topology changes do occur in higher dimension. In
addition, it seems difficult (at least not trivial) to generalize our algorithm in higher dimensions
even for hypersurfaces in a three dimensional space.
The paper is as follows. We first give in section x2 some preliminary definitions and introduce
some operators on sets. This provides operators on curves which are boundaries of sets. These
operators are consistent with curvature dependent differential operators. They satisfy some
monotonicity and continuity properties, allowing to extend them to real valued functions. We
prove some consistency results and give the proof of convergence without entering too much
into details since we prefer to focus on numerical applications. In x3, we adapt the operator
previously defined to the special case of evolution driven by a power of the curvature and in x4,
we give an algorithm with the same scaling covariance properties as the curve Scale Space. We
then show in x5 some numerical experiments in the convex case as well as in the nonconvex case.
2 Definition and Properties
We first give some notations (previously used in [18]). Let C a semi-closed curve in R 2 (that is
an oriented simple curve dividing the plane in exactly two connected components) and K the
interior of C (which is the bounded component of R 2 nC when C is closed). We suppose that C
is oriented such that K lies on "the left" when C is positively described. More rigorously, if we
assume that the plane is counterclockwise oriented, the inner normal is such that the tangent
vector and the inner normal form an orthonormal direct basis. We assume that a smooth
parameterization is defined on C (for example piecewise C 1 ).
A chord is a segment of the form ]C(s); C(t)[ that does not intersect any point of C with
parameter between s and t. A chord set C s;t is the connected set enclosed by a chord ]C(s); C(t)[
and the curve C(]s; t[). We say that C s;t is a oe-chord set and (s; t) is a oe-chord if L 2 (C s;t
and if the area of any chord set strictly included in C s;t is strictly less than oe. We denote by
K oe (C) the set of oe-chord sets.
a oe-chord of K and C s;t the associated oe-chord set. We call chord arc
distance of C (or of C s;t ) the number ffi(C([s; t]); [C(s); C(t)]) where ffi is the Hausdorff semi-
distance (in particular it is not commutative) defined by
For x in R 2 , we also denote by ffi s;t (x) the signed distance from x to the oriented line ]C(s); C(t)[,
that is
are in R 2 , we denote by [x; y] the determinant of the 2 \Theta 2 matrix with columns x and
y). We also denote by K
oe (C) the sets of positive oe-chords i.e the chord-sets C s;t satisfying
In the same way, we can define
oe the set of negative chord sets. Remark that for positive
chord sets, the chord-arc distance is nothing but sup ffi s;t (x) for x 2 C([s; t]) and for negative
chord sets, the chord-arc distance is
To finish with notations, we set
We now define a mapping on the set of plane sets that we shall assume piecewise regular for
sake of simplicity. Let G be a nondecreasing 1-Lipschitz function defined in R+ and such that
K be a set in R 2 whose oriented boundary
assume that L 2 (K) ? oe. For C s;t a chord set of K with chord-arc distance
h, we write
oe (C s;t
ae
!oe 2=3
'oe
In the sequel, we shall briefly say that  oe (C s;t ) is a modified chord set. We remark that the
right-hand term of the inequality above is nonnegative because of the Lipschitz assumption on
G. Hence, a modified oe-chord set is always included in its associated oe-chord set. On Figure 1,
we represent a oe-chord and its modified chord. The modified oe-chord set is filled.
PSfrag replacements
oe
!oe 2=3 G( h
s

Figure

1: A oe-chord set before and after transform.
K be the interior of a piecewise C 1 semi-closed curve. We define
oe (A): (4)
We will refer to E oe as an erosion operator and to E oe (K) as the eroded of K at scale oe.
Remark 2. The algorithm in [18] corresponds to
Remark 3. We can generalize Definition 1 to sets with several connected components by applying
the erosion to each component.
Lemma 4 Let K be a convex set of R 2 . Then E oe (K) is also convex.
Proof. We can write
Kn oe (A);
which proves that E oe K is convex as an intersection of convex sets (each of them is the intersection
between K and a half plane).
Lemma 5 If K is a smooth compact set, then E oe (K) is a compact set.
Proof. This is obvious since K is an intersection of compact sets.
The proposition below is crucial in a theoretical point of view as well in a numerical one since
it is necessary to obtain a stable numerical algorithm.
Proposition 6 (Inclusion Principle) Let K 1 ae K 2 . Assume that G is nondecreasing and
1-Lipschitz. Then
Proof. Assume that x 2 K 2 and x 62 E oe (K 2 ). We prove that x 62 E oe (K 1 ). If x 62 K 1 then
Assume now that x 2 K 1 . By assumption, there exists a
oe 0 -chord of K 2 (that we denote by C) with oe 0  oe such that x belongs to the modified chord-set.
The Lipschitz condition on G implies that x also belongs to the oe 0 -chord set. The same oe 0 -chord
delimits in K 1 a unique chord-set containing x. Let oe 00 be its area : we have oe 00  oe 0 and thus
oe 00  oe. It then suffices to prove that this chord excludes x from E oe (K 1 ). Consider the situation
illustrated in Figure 2.
are the chord-arc distances of C in K 1 and K 2 .
are the chord-arc distances of the associated modified chords.
(iii) l is the difference of length between K 1 and K 2 in the direction that is orthogonal to
the chords, i.e.
It is enough to prove that l + l 1  l 2 . But, we know that l
oe 2=3 !G( h 2
is 1-Lipschitz, we conclude.
PSfrag replacements
l
l 1
l 2
x
oe 00

Figure

2: Inclusion Principle.
Proposition 7 Assume that is of class C 2 . Let M 2 C such that the curvature of C
at M is not equal to 0. Then
lim
!oe 2=3
is the positive part of , defined by
Proof. Assume that C is concave at M (that is the curvature at M is strictly negative in
a neighborhood of M ). Then for oe small enough, any oe-chord (s; t) such that M 2 C([s; t]) is
a (strictly) negative oe-chord. Hence M 2 E oe (K) and the proposition follows from
Assume now that C is strictly convex at M (thus in a neighborhood of M ). Then for any chord
parallel to the tangent at M enclosing a oe-chord set containing M , the chord-arc distance h
satisfies
which can be easily established for a parabola, then for any regular curve by approximation (see
[18] for example). Then the result follows from 3 and 4, which imply
!oe 2=3
and the fact that G is continuous.
The following property is a continuity property allowing to extend E oe to any plane set, and
then to define an operator acting on functions with real values.
Proposition 8 (Continuity) Let K n a sequence of compact smooth sets. Let
Assume that K is also a smooth set. Then
Proof. Since K ae K n for any n, by monotonicity we have also
the first part of the equality. In order to prove the reverse inclusion, we can assume that the
family K n is nonincreasing. Without loss of generality, we also suppose that (K n ) converges to
K for the Hausdorff distance between compact sets. Assume that x 62 E oe (K). By definition,
there is a chord (s; t) with area not more than oe such that the modified chord excludes x. Since
oe (K) is closed, its complementary is an open set; thus we can assume that the area of C s;t is
strictly less than oe. In K n it also defines a chord and for n large enough, the area of this chord
is also less than oe (by using convergence of measures). Moreover, as K n tends to K for the
Hausdorff distance, the chord-arc distance also converges. Since G is continuous, this implies
that the chord excludes x in K n for n large enough. Hence x 62 "E oe (K n ) and this ends the
proof.
Remark 9. The compactness assumption in the proposition above is far from necessary. It suffices
for example that the boundary of the sets is locally convex or concave. This ensures that the
erosion is local when oe is small. We can then conclude by the same kind of arguments. Note
also that the continuity property allows to define the erosion on any closed set by approximating
closed sets by smooth closed sets.
We can now extend E oe to real-valued functions in R 2 . First, if u R, we define as usual
the level set with value , the set
Applying a theorem by Matheron (see [14]) yields the following
F a set of real valued functions in R 2 such that the level sets of the elements
of F are compact and smooth. Then, we can extend E oe to elements of F by setting
It is equivalent to define E oe (u) by its level sets,
This uniquely defines a monotone, translation invariant operator commuting with nondecreasing
continuous functions.
We also define a dual operator D oe (called dilation operator) by
the subscript designing the complementary set in R 2 . This operator satisfies the same properties
as E oe except the consistency result where the positive part of the curvature has to be replaced
by the negative part.
By standard arguments ([4], [10], [20]), we can derive the following consistency result on the
operator acting on functions.
Proposition 11 Let u be a C 3 function and x a point such that Du(x) 6= 0 and (u)(x) 6= 0.
Then
Proof. The whole proof is not very difficult but a bit technical and long. Thus, we do
not enter into all details and we shall skip some points. The aim is to prove that E oe (u) only
depends on local features of u. Choose ff such that oe
oe tends to 0 (note that there is no incompatibility; choose r = oe 1=4 for instance). By using
translation invariance and contrast invariance, we assume that
of the proof is the following. If r is small enough, the curvature of the level lines of u has a
strict sign in D(0; r). As a consequence, for any oe-chord, we can estimate the chord-arc distance
by Equation (5). Moreover the same kind of approximation (made on a circle or a parabola)
shows that the length of the chord is of order ( oe
) 1=3 . Hence, any oe-chord intersecting D(0; r
asymptotically included in D(0; r(1 (because of the choice of r). Assume
first that (u)(0) ? 0. Define u+ and u \Gamma by
and u+ elsewhere (we can replace the infinite value by very large
numbers). The global inequalities
enough, the level lines of u are uniformly strictly concave in D(0; r). Thus, there is no positive
oe-chord of level sets of u \Gamma and u+ intersecting D(0; r
Hence
r
Assume now that  ! 0. Define
in D(0; r) and
in D(0; r) and elsewhere. The constant k is chosen such that for oe small enough,
we have
PSfrag replacements

Figure

3: Case  ? 0. On the left, a level set of u \Gamma . The oriented boundary is the bold line. The
level set is the bounded connected component delimited by the curve. The dashed line is the
boundary of the eroded set. On the right, a level set of u+ : the boundary is the bold oriented
line and the level set is the unbounded component. In both cases, there is no positive oe-chord
set intersecting D(0; r
Hence, the erosion has no effect.
This is possible since we assumed that u is C 3 . Monotonicity yields
As v and w have trivial level sets out of D(0; r), it is quite easy to estimate their image by E oe .
We now use the consistency result (Proposition 7). The only trick is that the level lines of u
and v are not parabola in the canonical form. Nevertheless, with some few arguments we can
be led back to this situation ([4, 10, 20]). The computation of the eroded level sets of v and w
is drawn on Figure 4. From this, it is no longer difficult to prove that
and
are constants depending on D 2 u(0). We can apply the same result to the
dilation operator D oe to obtain the second part of the proposition.
For the next proposition, we first extend G to make it odd, that is if x ! 0, we set
\GammaG(\Gammax).
Proposition 12 Let u be a C 3 function. Suppose that Du(x) 6= 0 and (u)(x) 6= 0. Then
Proof. This follows from the fact that : 1) E oe and D oe are monotone and commute with
addition of constants, 2) near a point with gradient and curvature different from zero, the
arguments developed in the previous proposition are uniform.
In order to be complete, we just give without proof (it is simple when introducing the affine
erosion operator ([18])) a lemma controlling the behavior of E oe at critical points.
PSfrag replacements

Figure

4: Case  ! 0. On the left, a level line of v (oriented bold line). The eroded line is the
dashed line. On the right, the same thing with w.
lim
oe 2=3
the limit being taken when y and oe tend to 0 independently.
This and the consistency result above allow to deduce the convergence result we now enounce
(see [21]). We denote T
Theorem 14 Let u 0 bounded and uniformly continuous in R 2 . For oe ? 0, let define
R 2 \Theta R! R by
Then, when oe tends to 0, u oe converges locally uniformly to the unique viscosity solution of the
equation
@t
The usual definition of viscosity solution ([6]) is no longer valid here because of the possible
singularity of the operator at critical points. An extended notion of solution is defined in [3]
and the solution still exists and is unique.
3 Curvature Power
3.1 Approximation of power functions
In this section, we show that the previous study can be adapted to the particular case of power
functions F not 1-Lipschitz on the whole real line,
we define
where
is the largest positive number at which the power function x 3fl has a derivative less than 1. As
G is 1-Lipschitz, we can then apply E oe to this function. We could think that this scheme is not
consistent with motion by curvature power. Indeed, if the curvature is too large for fixed oe then
it may happen that the chord arc distance is also very large and the erosion is then given by the
linear part of G. Nevertheless, we shall see that by an adequate scaling, G is not evaluated in
its linear part. From now on, we slightly change the notations in the definition of E oe . This will
simplify the statements in the case of power functions. For a oe-set C s;t , we now set
oe (C s;t
and we still define the erosion operator by
oe (A): (23)
When the scale tends to 0, h also tends to 0 and G is not taken in its linear part. The fact that
we get an operator consistent with a power function is due to the homogeneity properties of
power functions (this implies that except for power functions, this new definition of the erosion
operator makes no sense). We can adapt the proof of the inclusion principle 6 to prove that
the modified erosion operator E oe still satisfies this inclusion principle. Consistency on curves
(circles is enough !) is easy to establish. Continuity is not a problem as well. Thus, by using
Matheron's Theorem, we can extend this erosion operator to an operator acting on functions.
The following proposition asserts that this operator is consistent with a power of the curvature.
Proposition 15 Let u be a C 3 function. Let x be such that Du(x) 6= 0 and (u)(0) 6= 0. Then
Moreover, consistency is locally uniform.
Proof. As usual, we assume that In a first time, assume that (u)(0) ? 0.
We use the same locality argument as in the first proof of consistency in this chapter. Since the
curvature of u in a small ball with radius r ? 0 is bounded from below by a positive constant,
say the level sets of u have no positive oe-chord in D(0; r). We set
We use the same locality argument as in the first
proof of consistency as above: For oe small enough, the level sets of u \Gamma have no positive oe-chord
in D(0; r). Thus the erosion has no effect upon u \Gamma in D(0; r). By using monotonicity, we have
and the result is proved in the case  ? 0.
Let us now come to the most difficult case: Assume that
cxy be the Taylor expansion of u at the origin with
parameter and let
If r is chosen small enough, we have v(x)  u(x) in D(0; r). By extending v by \Gamma1 out of
D(0; r) this remains true everywhere. Moreover, we can assume that the curvature of the level
lines of v is still strictly negative. Indeed, its value is  is the value at the
origin. We want to approximate E oe (v)(0). Let now j ? 0 be also small, such that the curvatures
of the level lines of v is larger than that in this part of the proof the curvatures
are all negative, hence a circle with a small curvature will also have a small radius). We can
again invoke the same locality property of the erosion in the case of a curve with a strictly
negative curvature: We know that the chord-arc distance is equal to O(oe 2=3
the length of the chord is a We deduce from this that the oe-chord sets
containing 0 must be included in a ball with radius O(oe 1=3 ). A small computation shows that
the modified corresponding oe-chord sets are included in a ball with radius O(oe fl ). The constant
in these terms are clearly uniform because the curvature is bounded from above by a negative
constant. In particular, they do not depend on " and j. Let now x be in a D(0; r). We call
C \Gamma"\Gammaj (x) the disk of curvature that is tangent to the level line of v at x and that is
in the same side as  v(x) (x). Because of the comparison of the curvatures, we can still assume
that r is small enough such that we have the inclusion
Now, since the erosion operator is local, we also have
Assume that E oe (v)(0)  . By definition, this means that 0 62 E oe (  (u)). By using inclusion
principle, we deduce that 0 62 E oe (C \Gamma"\Gammaj (x)) for any point such that . On the level line
of v with the same value , we can find a unique point x  such that x  and the normal at x
are colinear (this is due to the strict convexity of the level lines and the theorem of intermediary
values). This point is also characterized by the fact that the distance between the origin and
the tangent to the level line is minimal. Hence, the modified oe-chord set of C \Gamma"\Gammaj
also contains the origin. Let be the coordinates of x  . A simple calculation gives
Since x  and Du(x  ) are colinear, we have
By using consistency on disks, we have jx
From this, we finally deduce
where we have approximated jDv(x  )j by p up to a O(oe 2fl ) term. This analysis can be performed
since eventhough the constant were not explicited, we already stress
that they do not depend upon " and j. We then deduce that
Let now search an upper bound to E oe (u)(0). We do not repeat all the arguments since there
will be some similarity with the research of a lower bound. We approximate u by its Taylor
expansion and define 2cxy such that u  w is a small ball
of radius r with
enough, the curvature of the level lines of w is
smaller than which can also be chosen negative if r is small enough. The locality of the
oe-chords still holds. We now define C +"+j (x) as above; its radius is equal
r is small enough, for any x the level set  w(x) (w) is included in C +"+j (x) inside D(0; r). The
rest of the proof is still an application of comparison principle and the asymptotic behavior on
erosion on disks. Assume that E oe (w)(0)  . This means that In particular
is a above. This implies that x   ! 3fl
We use the characterization (24) of x  and deduce that we must have
Since, this is true for any " ? 0 and j ? 0, we also obtain
Again, we can pass to the limit since the o(1) term is uniform in " and j, since all the curvatures
may be taken bounded by above by a strictly negative constant. Thus, if  ! 0, we have
With the case   0, this gives the result.
The case of the dilation can be deduced by the relation D oe (\Gammau). The uniform
consistency follows from the fact that the oe-chords are uniformly bounded in some ball with
radius O(oe 2=3 ) and the constants of these terms are bounded as soon as the curvature have an
absolute value strictly more that a positive constant.
Uniform consistency yields consistency for the alternate operator.
Corollary
At critical points, we need to describe the behavior of the erosion in the following manner. Let
R be a C 2 function such that
tends to 0.
Lemma any sequence of points x n tending to 0, we have
lim
oe 2fl
where the limit is taken as n and oe tends to their respective limit independently and T h designs
either E oe or D oe or D oe
This lemma can be easily established by finding some estimates on the radius of the circle after
erosion. We can then prove the following convergence theorem.
Theorem
3 . Let T the alternate dilation-erosion for the curvature
power function x fl . Let u 0 in BUC(R 2 ) and define u oe by
Then, when oe tends to 0, u oe tends locally uniformly to the unique viscosity solution of the
equation
@t
with initial value u 0 .
As soon as the power fl is more than 1, then the usual notion of viscosity solution is not
appropriate since the elliptic operator jDuj(curv u) fl is singular at critical points. Ishii and
Souganidis proved in [12] that existence and uniqueness where still true if test functions were
restricted to a class of functions with flat critical point. This flatness is given by the same
conditions in the previous lemma. This point apart, the scheme of the proof is standard so we
do not explicit it. The only new point is the previous lemma in the case where test functions
are stationnary. In this case, the lemma directly gives the solution and we leave the rest of the
proof to the reader.
3.2 Scale Covariance
We denote by H  the dilation with ratio , that is H  t be the evolution
semi-group of
@t
that is, is the curve evolving according to Equation (25) above (the solution exists
and is unique at least for short times for smooth initial data (see [25]). The semi-group S t
satisfies the relation
Indeed, let C 1 (t) the evolving curve defined by
Then, it is simple to check that C 1 satisfies Equation (25) with initial condition H  C(0). This
is exactly what Equation (26) asserts. The erosion operator E oe does not satisfy the same
covariance property. We thus define a modified operator
O
a 2
where a and oe are positive parameters depending upon t (and possibly on C). We want O t to
satisfy covariance Property (26), like S t . This will be true as soon as
and oe(H  C;
We also want O t to be consistent with S t , that is, for any convex set K with C 3 boundary,
O t
the term o(t) being measured with the Hausdorff distance (we also use an abusive notation by
denoting S t (K) the set whose boundary is S t (@K)). By using consistency result above, we see
that a, oe and t must be linked by the relation
Assume that oe ? 0 is fixed. If a is chosen large enough such that, for any oe chord of K with
chord-arc distance h the inequality
a
ff fl (30)
holds, then the modified chord-arc distance is then given by G near the origin, that is Equation
(21) is not involved. Precisely, for M 2 @K (with smooth boundary) consistency writes
a
a
a
Notice that it is interesting to take the smallest possible value of a (given by the case of equality
in (30)) in order to get the largest possible scale step t from (29). We can summarize these
results in the following
Proposition 19 Let h(A) denote the chord-arc distance of a chord A. Then, the operator O t
defined by (27) with
sup
t! \Gamma3fl a 3fl \Gamma1 \Delta 1=2fl (32)
is consistent with (25) and satisfy the same scaling property as S t in (26).
Remark 20. (Error Analysis) In order to obtain consistency, it is not necessary that h
a
ff fl for
all (oe; a); for instance, if we fix a and define O t by (27) and (29), we still have consistency since
inequality (30) holds for oe small enough. However if a and oe do not satisfy (30), the difference
between O t C and S t C is of the same magnitude as t (since Equation(21) is involved). On the
contrary, this difference is O(t 2 ) if a and oe satisfy (30) (because of Consistency 5).
4 Algorithm
4.1 General method
Each iteration of the operator defined above involves three parameters: the scale step t (that can
be viewed as a time step), the erosion area oe and the saturation length a. These three quantities
have to satisfy (32), which leaves only one degree of freedom. A usual numerical scheme would
consider t as the free parameter (the time step, related to the required precision), and then
define a and oe from t. In the present case, this would not be a good choice for two reasons.
First, a is defined as an explicit function of oe but as an implicit function of t, which suggests
that oe may be a better (or at least simpler) free parameter than t. Second, t has no geometrical
interpretation in the scheme we defined, contrary to oe which corresponds in some way to "the
scale at which we look at the curve". In particular, oe is constrained by the numerical precision at
which the curve is known: roughly speaking, if C is approximated by a polygon with a precision
" (corresponding, for example, to the Hausdorff distance between the both of them), then we
must have oe AE " 2 in order that the effect of the erosion at each iteration overcomes the effect
of the spatial quantization. For all these reasons, we choose to fix oe as the free parameter, and
then compute t and a using (32). If the scale step t obtained this way is too large, we can simply
adjust it by reducing oe while keeping the same value of a. We propose the following algorithm
for the evolution of a convex set K at final scale T with area precision oe.
1. Let
2. While
- For each oe-set of K t , compute the chord-arc distance.
- Set a to the maximal value of these distances.
decrease oe in order to keep the previous equality.
Apply operator O t to K t , yielding K t+ffit .
Increment t by ffi t.
In practice, it is of course impossible to deal with all the oe-chords. In fact, the curve is a
polygonal line and we take the chords with an end point equal to a vertex of the polygon.
4.2 Computation of the erosion
The boundary of O t (K) is included in the envelop of all the modified chords. To obtain an
approximation of this set, we explicitly determine the position of the unique point of each chord
belonging to the envelop. This result is a generalization of the middle point property exposed
in [18],
Lemma 21 (Middle point Property) Let K be a strictly convex set. Let A 1 be a oe-set of
K with oe-chord C 1 . Let A 2 be another oe-set, and let C 2 be its oe-chord. Then, when dH
tends to 0, the intersection point of C 1 and C 2 tends to the middle point of C 1 .
Proof. (quoted from [18]). Let ' be the geometrical angle between C 1 and C 2 . If A 1 and A 2
are close enough, C 1 and C 2 intersect at a point that we call I('). We also call r 1 (') and r 2 (')
the length of the part of the chord C 2 on each side of I('). Since C 1 and C 2 are oe-chords, we
PSfrag replacements

Figure

5: The middle point property.
This implies that lim(r 1 tends to 0.
The boundary of O t (K) is included in the envelop of the "modified" oe-chords. In the case
these chords are the oe-chords themselves. For all other values of fl, this is no longer
true and we have to compute the position of the intersection of closer and closer chords. This is
the purpose of
Proposition 22 Let K be a strictly convex set, and a oe-chord of K. Consider P the
farthest point of C([s; t]) from C. If (L; h) are the coordinates of P in the direct orthonormal
referential whose origin is the middle point of C and whose first axis is directed by C(t)-C(s),
then the contribution of the modified chord arising from C is either void or the unique point with
coordinates
a
a
in the same referential.
Proof. Examine the situation on Figure 6. Let C ' the oe-chord making an angle ' with C. We
search the coordinates of the intersection point of the modified chord of C and C ' when ' tends
to 0. We set x(') the abscissa of the point . By the middle point property, we know that
' the modified chord of C ' . The distance between these two chords
is
a
is the chord-arc distance of C ' . Let (L('); H(0)) be the
coordinates of the common point of C 0 and C 0
' . Elementary but a bit fastidious geometry proves
PSfrag replacements'

Figure

The modified middle point property.
that
a
implying that the limit point we are looking for has coordinates given by (33).
Remark 23. In the case of the general function of the curvature and without introducing the
scaling preventing saturation phenomenon, the coordinates of the point are
!oe 2=3
!oe 2=3
As precised in the proposition, the limit point may not belong to the boundary of O t (K). Indeed
@O t (K) is in general strictly included in the envelop of the modified oe-chords. In general, this
envelop is not even the boundary of a convex set! Nevertheless, if we know that C is a convex
curve then it is simple to decide whether a point has to be kept or not by comparing its position
with adjacent modified chord. Hence we can remove the bad points and obtain a convex set.
For example, on Figure 7, we display the envelop of the modified oe-chords of a square. If the set
is a "corner", the explicit computation can be made and shows the same behavior in the corner
as in

Figure

7. The eroded set is obtained by removing the parts with cusps in the corners of
the square.
5 Numerical Experiments
To finish, we display numerical experiments, first in the case of convex sets. By a change of
scale variable, we implemented an approximation of the equation
@t
For this rescaling, scale and space are homogeneous (precisely, if T t maps C to C(t) by this
equation, we check that T t

Figure

7: Envelop of the modified oe-chords for a square 2). The result is not convex. The
bold line is what is to be kept.
5.1 Closed convex curves
The first example is the the case of circles. The radius is explicitly computable for the scale
space since
Remark that the extinction scale for a circle with initial radius R(0) is R(0) for any fl. On

Figure

8, we display the evolution of a circle with radius 10, for 9.

Figure

8: Evolution of circles for displayed at scale 0, 1,., 9 (fast computation).
In

Table

1, we give the theoretical and computed radius for several values of fl and for a circle
with initial radius equal to 10. We performed two sets of experiments with different precisions,
and for each of them we give the CPU time on a Pentium II 366MHz, the number of performed
iterations and the final obtained radius. In Figures 9 and 10, we display the evolution of convex
closed polygons. On each figure, the display scales are the same for all the different values of fl.
5.2 Unclosed curves evolution
Until now, we have only studied the evolution of compact convex sets. The boundary of such
sets is a closed convex curve. It is possible to make nonclosed convex curves evolve by fixing their
end points as in [5]. This is equivalent to symmetrize and periodize the curve. The steady state
to this evolution is a segment if the two ends are disjoint. If they are equal, then a singularity
fast computation slow computation
R theo R comp # iter CPU
2.0 6.47 6.39 112 1.22 6.45 421 27
3.0 7.65 7.56 108 1.38 7.64 400 28
10.0 9.66 9.62 100 1.19 9.66 214 15

Table

1: The shortening at scale 9 of an initial circle with radius 10 is considered. For different
values of fl, we give the theoretical value of the final radius (R theo ), the corresponding values
obtained with different precisions (fast or slow computation), the number of performed iterations
and the used CPU time.

Figure

9: Evolution of a triangle. Up-left:
occurs in finite time. Such an evolution is displayed on Figure 11 for a nonclosed convex curve
for several values of fl.
As can be derived from [8, 9, 2], if C moves by Equation 25 and is locally the graph of a
function, y then y satisfies the equation
We use this equation to determine whether a convex curve with distinct fixed end points becomes
a straight line in finite time.
Proposition 24 Let u 0 be a strictly convex function on [\Gamma1; 1] with u

Figure

10: Evolution of a pentagon. Up-left:
be the solution of
Then, if identically zero in finite time. If fl  1, then the steady state is
attained for infinite time.
Proof. From Equation (37), we deduce that u is subsolution of the following equation
On the other hand, we can derive from Equation (37) an equation for u 0 . This equation is also
parabolic. By maximum principle, the supremum of u 0 is attained at time
that u is supersolution of
For both Equations (38), (39), we can compute separable solutions of the type g(t)f(x). Both
functions f and g then satisfy an ordinary differential equation that is explicitely solvable for
g and can be expressed in terms of elliptic functions for f . We can then check that the time
components becomes null in finite time if and only if fl ! 1. In order to conclude, it suffices to
bound u 0 fromabove if 1 and from below if fl  1 by a adequate f (the spatial part of the
separable solution) and to apply maximum principle.

Figure

11: Evolution of an angle for different powers of the curvature
2:0). The displayed curves correspond to scales that are integer multiples of
a same fixed value.
Let us interpret this experiment in terms of image processing. Imagine that we process the
filtering of a grey level image by applying (25) to its level lines, as done in [13] for
This will have smoothing effects, and in particular one can expect to remove pixellization effects.
The periodic structure corresponding to the initial state of Figure 5.2 is a "staircase" line that
corresponding to the discretization of a perfect straight line (oriented at 45 ffi ). As the evolution
scale increases, this infinite staircase is smoothed and eventually becomes a perfect straight line
in finite time if fl ! 1. A natural question now arises: how choose fl in order to smooth these
staircase effects with the smallest possible damaging effects on the image? In other terms, what
power of the curvature regularizes discrete lines in the shortest time? The previous proposition
asserts that only powers smaller than 1 can straigthen a staircase in finite time. Experiments of

Figure

5.2 corroborate this result and indicate that the straightening time increases with fl. But
this result has to be counterbalanced in the following. By smoothing an image, we would like
to remove small undesirable details while keeping the rest of the image unchanged. The results
and experiments on circles (Figure 8) tend to prove that large powers can do this in a better
way than small powers (i.e. when a circle with initial radius 9 disappears at scale
with radius 10 is less changed for large values of fl).
5.3 Generalization to nonconvex sets
An algorithm for nonconvex curves has been proposed in [18] for the affine erosion
apply the same method in the general case. It consists in splitting a curve into its convex and
concave components. This decomposition is unique and well defined. By this method, we do
not find inflexion points but inflexion segments and we define an inflexion point as the middle
point of an inflexion segment. We then apply the erosion operator to all the convex components
by fixing the end points (which are the inflexion points defined above). Once this is done we
gather all the parts to form a new curve. We reapply the decomposition and the erosion to
this new curve. Notice that near an inflexion point, the ending segments of the joining convex
components do not stay parallel in general; thus, inflexion points have have no reason to stay
still. Practically, this is what we observe. Moreover, as fl increases, they seem to move more
and more slowly, which also seems logical. On Figures 12 and 13, we display the result of the
algorithm on nonconvex curves. Nevertheless, we do not have any strict justification for the
convex component decomposition and the displacement of inflexion points should be studied
more carefully.

Acknowledgements

. The authors would like to thanks Jean-Michel Morel for all valuable
conversations and advice.

Figure

12: Scale space of a "T" shape: the curves are displayed at the same evolution scale.
From left to right and top to bottom: fl =0.4, 1, 2, 3. CPU times are respectively 5, 9, 23 and

Figure

13: Scale space of a hand curve: the curves are displayed at the same evolution scale.
From left to right and top to bottom: initial curve,fl =0.4, 1, 2. CPU times are respectively 6,



--R

Axioms and fundamental equations of image processing.
On the affine heat equation for nonconvex curves.
Convergence of approximation schemes for fully nonlinear second order equations.
Partial differential equations and mathematical morphology.
Is scale space possible
User's guide to viscosity solution of second order partial differential equations.
Motion of level sets by mean curvature.
The heat equation shrinking convex plane curves.
The heat equation shrinks embedded plane curves to round points.
Partial Differential Equations and Image Iterative Filtering.
An approximation scheme for gauss curvature flow and its convergence.
Generalized motion of noncompact hypersurfaces with velocity having arbitrary growth on the curvature tensor.
Geometric multiscale representation of numerical images.
Random Sets and Integral Geometry.
Solution of nonlinear curvature driven evolution of plane convex curves.

Traitement num'erique d'images et de films
Affine plane curve evolution: A fully consistent scheme.
Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulation
Approximation de propagation de fronts avec ou sans termes locaux.
Approximation of viscosity solution by morphological filters.
Affine invariant scale space.
On affine plane curve evolution.
Curvature and the evolution of fronts.
The generalized curve shortening problem.
--TR

--CTR
L. Alvarez , A.-P. Blanc , L. Mazorra , F. Santana, Geometric Invariant Shape Representations Using Morphological Multiscale Analysis, Journal of Mathematical Imaging and Vision, v.18 n.2, p.145-168, March
curve evolution;image processing;level set methods;viscosity solutions;numerical approximation