and fourth-order geometric partial differential equations for space curves
Computer Aided Geometric Design 31 (2014) 63–80
Contents lists available at ScienceDirect
Computer Aided Geometric Design www.elsevier.com/locate/cagd
Construction of several second- and fourth-order geometric partial differential equations for space curves ✩,✩✩ Guoliang Xu a , Xuyang Yang b a
State Key Laboratory of Scientific and Engineering Computing, Institute of Computational Mathematics, Academy of Mathematics and System Sciences, Chinese Academy of Sciences, Beijing 100190, China China Shipbuilding Industry Corporation, 750 Proving Ground, Kunming 650051, China
b
a r t i c l e
i n f o
Article history: Received 29 December 2012 Received in revised form 12 December 2013 Accepted 26 December 2013 Available online 7 January 2014 Keywords: Geometric partial differential equations Geometric flows Space curves
a b s t r a c t Geometric partial differential equations for curves and surfaces are used in many fields, such as computational geometry, image processing and computer graphics. In this paper, a few differential operators defined on space curves are introduced. Based on these operators, several second-order and fourth-order geometric flows for evolving space curves are constructed. Some properties of the changing rates of the arc-length of the evolved curves and areas swept by the curves are discussed. Short-term and long-term behaviors of the evolved curves are illustrated. © 2014 Elsevier B.V. All rights reserved.
1. Introduction In many design and construction problems, such as reverse engineering, multi-component mesh optimization and industrial geometric design from a wireframe, smooth space curves need to be generated from noisy input data. This goal can be efficiently achieved by utilizing geometric flows to smooth the input data with shape preserved. The aim of this paper is to construct such geometric flows for space curves embedded in R3 . For surfaces embedded in R3 , many geometric flows have been constructed (see Xu, 2008 for references). These include second-order geometric flows (mean-curvature flow, averaged mean-curvature flow, weighted mean-curvature flow, Gaussian curvature flow, etc.), fourth-order flows (surface diffusion flow, quasi-surface flow, Willmore flow, etc.) and sixth-order flows (minimal variation of mean-curvature flow, quasi-minimal variation of mean-curvature flow, etc.). These flows are frequently used in surface modeling, surface denoising, image processing and computer graphics. As one can imagine, all the plane curves in the xy-plane can be extended to the xyz-space by setting the z-direction as constants, obtaining three-dimensional (3D) cylindrical-type surfaces. Then, the mean curvature H of this surface is half of the signed curvature of the plane curve, and the Gaussian curvature K is zero. Then, all the geometric flows for surfaces have their counterparts for plane curves. For this reason, one may think that there is no need to study the construction of geometric flows for curves. However, for space curves in R3 , reducing the surface flows to curves is not always valid, as explained below. Since most of the curvature flows for surfaces embedded in R3 describe the surface motion in the normal direction, we take the normal as an example to explain the difference. For a single-sided surface, the surface normal can be orientated, but the normal of the space curve cannot be orientated. For instance, a straight line in 3D, as a counterpart of the plane, cannot be orientated. Furthermore, the motion of the surface in the normal direction (the tangential motion does not change ✩
This paper has been recommended for acceptance by Konrad Polthier. Supported in part by NSFC under the grants 1101401 and 81173663, NSFC key project under the grant 10990013 and Funds for Creative Research Groups of China (grant No. 11021101). ✩✩
0167-8396/$ – see front matter © 2014 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cagd.2013.12.005
64
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
the surface shape) has one degree of freedom. However, for a space curve, the motion has two degrees of freedom. Hence, the situation for the space curve is rather different from that for the plane curve. In this paper, we consider the construction of geometric flows for space curves. 1.1. Previous work on geometric flows for plane curves There has been much work on plane curve evolution (see Gage, 1983, 1984; Gage and Hamilton, 1986; Ma and Chen, 2007; Sapiro, 2001 for references). We only review a few of them. Most of the curve evolutions focus on the following curve-shortening flow (CSF)
∂x = κs N ∂t and its various generalizations. Here, κs is the signed curvature and N is the unit inward normal of a close plane curve x. CSF is also known as curvature flow or geometric heat flow, which is the counterpart of the mean-curvature flow for surfaces. In Chou (2001), some results on CSF are presented, including the properties of CSF, the anisotropic generalized CSF and the CSF in Klein geometry. In Andrews (1998), the behavior of convex curves driven by the curvature flow is considered. The long-term behavior of solutions and properties of limiting shapes are described. As shown in Gage and Hamilton (1986) and Grayson (1987), the curvature flow will shrink an embedding planar curve to a point, becoming round in the limit without developing singularities. In Dziuk (1999), an anisotropic CSF is studied, and a numerical scheme for this flow is developed. Optimal error estimates for a semi-discrete scheme are provided. A weighted curvature flow for planar curves, a generalization of the curvature flow, is discussed in Chen et al. (2007). The convexity of the curve driven by the weighted curvature flow is studied. The geometric curve evolution equations in a plane have been extended to surfaces and manifolds in Grayson (1989), Gage (1990), Spira and Kimmel (2007). Efficient numerical schemes for the implementation of the classical geodesic curve evolution equations on parametric manifolds are presented in Spira and Kimmel (2007). These flows are used on the scale space of images painted on manifolds and segmentation by an active contour model. Grayson (1989) and Gage (1990) generalized the study of curvature flow of closed curves in a plane to that in surfaces, and analyzed the long-term behavior of smooth curves immersed in a surface. 1.2. Length- or area-preserved flows for plane curves It is well known that geometric heat flow is the fastest way to shrink the length of a curve (Gage, 1986). A planar embedded curve that evolves according to this flow becomes more and more circular as t increases, and eventually shrinks to a point in finite time (Gage and Hamilton, 1986; Grayson, 1989). Such a shrinking effect is an undesirable feature in many applications, such as curve denoising, since the evolving curve will drift rapidly from the desired shape. In order to avoid the shrinking effect, Gage (1986) proposed two variations of the geometric heat flow. The first one is defined by
∂x 2π = κs − N, ∂t L
(1.1)
where L is the arc-length of the curve. The second one is defined by
∂x π x, N = κs + N + α T, ∂t A
where A is the area enclosed by the curve and α T is a tangential component that does not change the shape. The resulting flows are curve-shortening motions that preserve the area of the region enclosed by the curve. In Sapiro and Tannenbaum (1995), area-preserving multi-scale representations of planar curves are described. The representations are obtained by deforming the curve via geometric heat flows while simultaneously magnifying the plane for keeping the enclosed area constant. In Dolcetta et al. (2002), several models in phase separation theory are used to generate planar curve evolution without shrinkage. The theoretical validation of a level-set method to these models is discussed. In Ma and Zhu (2012), a length preserving curve flow
∂x 1 = κs − α (t ) N with α (t ) = ∂t 2π
κs2 ds
(1.2)
for convex curves in the plane is considered. It was shown that the area of the region bounded by the evolving curve increases, and the evolving curve converges to a circle. Pan and Yang (2008) proposed another length-preserving curve flow for convex curves in the plane of the form
∂x = ∂t
L 2π
− κs−1 N.
They proved that the convex plane curve will become more and more circular and converges to a circle.
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
65
1.3. Previous work on space curves Space curve smoothing or denoising is often required in reverse engineering, CAD modeling, animation, computer vision, image registration and mesh optimization (Leng et al., 2011). To smooth a discrete space curve (a sequence of points) with noise, B-splines were used in Guziec and Ayache (1994) to approximate the discrete data. In Dierckx (1995), noisy scattered data were smoothed by minimizing a functional that combines a term that measures the deviation of the data points from the fit with another term that measures the non-smoothness of the fit. An efficient fairing technique for a wireframe defined by two sets of space curves was presented in Liu et al. (2008). In the fairing process, the parameters are fixed and then released by an iterative gradient descent optimization method. In Cohen and Tekumalla (2005), a 3D curve-smoothing algorithm using the moving least squares (MLS) projection is proposed with the primary focus on smoothing boundaries of point clouds obtained during reverse engineering of laser scanned models. An algorithm for fairing of space curves with respect to spatial constraints based on a vector-valued curvature function is presented in Hildebrandt et al. (2005). It describes a curve-smoothing flow that satisfies strict spatial constraints and allows simultaneous control of both curvature functions. In Altschuler and Grayson (1992), the authors show how to use the curvature flow on space to define a natural continuation of the planar solution for the whole of the time. 1.4. Main contributions of this paper Our attention is focused on constructing geometric flows for space curves. First, gradient, divergence and Laplace operators on space curves are introduced. Then several second-order geometric flows and fourth-order geometric flows for space curves are constructed. We consider second- and fourth-order flows for the following reasons. First, second-order flows are easier to solve and can be used to denoise efficiently. However, the second-order flows can only achieve G 0 continuity for joining space curves. To achieve G 1 continuity, fourth-order flows are required. Second, for smoothing space curves with shape preservation, fourth-order flows are more desirable since the solution curves of the higher order flows are “softer” than those of the lower order flows. Some properties of the changing rates of the arc-length and areas are discussed. To describe the area changing rate for an evolving space curve, a notion of the area swept surface of the evolving curve is introduced. 1.5. Organization of the paper The remainder of this paper is organized as follows. In Section 2, we introduce some of the notation used and several differential operators on space curves, including the definitions of curvature, normal, gradient, divergence and Laplace operators. Some basic results for these operators are also presented. In Sections 3 and 4, we derive several geometric flows, including second-order and fourth-order flows. Some properties of these flows are discussed. Finally, we conclude the paper with a summary in Section 5. 2. Preliminary In this section, we introduce some of the notation used and several differential operators on space curves, including the definitions of curvature, normal, gradient, divergence and Laplace operators. Some basic results for these operators are also presented. 2.1. Basic differential geometry on space curves 2.1.1. Curvature and normal For a given smooth space curve C defined by
T
C = x(u ) = x(u ), y (u ), z(u ) let
: u ∈ Ω = [a, b] ,
κ be the curvature, which is given as (see Do Carmo, 1976, p. 25) x (u ) × x (u ) x (u )2 x (u )2 − (x (u )T x (u ))2 κ (u ) = = , x (u )3 x (u )3
(2.1)
where x (u ) and x (u ) denote the first- and second-order derivatives of x(u ) with respect to u. Note that κ 0. This is different from the mean curvature of a surface, which can be either positive or negative. The principal normal is given as
n(u ) =
(x (u ) × x (u )) × x (u ) x (u )x (u )2 − x (u )x (u ) T x (u ) . = (x (u ) × x (u )) × x (u ) x (u ) x (u )2 x (u )2 − (x (u )T x (u ))2
(2.2)
It should be noted that n(u ) is well defined for all u such that x (u ) and x (u ) are linearly independent. If x (u ) and x (u ) are linearly dependent, n(u ) is arbitrarily defined. We artificially define it as a zero vector. The curvature κ = 0 is
66
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
well defined. From (2.2), we can derive that n(u ) does not depend on the re-parameterization of the curve. Hence, n(u ) is uniquely defined. In addition to the curvature and normal, we define the curvature vector as
κ (u ) =
x (u )x (u )2 − x (u )x (u ) T x (u )
x (u )4
(2.3)
.
It is easy to see that
κ (u ) = κ (u ),
n(u ) =
κ (u ) if κ (u ) > 0. κ (u )
Hence, n(u ) is not continuous at the point
n (u ) = where
κ (u ) = 0. To obtain a continuous normal function, we define
κ (u ) , κ (u ) +
is a small positive number. In the following, all the stand for the same number.
2.2. Differential operators on space curves 2.2.1. Laplace operator Suppose f is a C 2 smooth function defined on the curve x(u ), then we define Laplace operator c (here the subscript c stands for that the operator is defined on curve) applying to f as
c f = divc (∇c f ),
(2.4)
where the gradient and divergence operators on a space curve are defined as
∇c f =
x f
x 2
divc v =
1
= d
x (x ) T ∇ f
x 2 T v x
x du x
,
=
[(v )T x + vT x ]x 2 − (x )T x vT x . x 4
It is easy to see that ∇c f is the projection vector of ∇ f on the tangent space of the curve. If u is arc-length and v is a tangent vector field, then
T ∇c f = x x ∇ f ,
T
divc v = v
x .
T
If we take f as the vector value function x, we have ∇c x = x(xx )2 . It is easy to see that for any tangent vector v, we have
(∇c x)v = v. Hence, ∇c x acts as an identity operator on the tangent space of the curve. From (2.4) we can derive that
c f (x) = =
f (x)x 2 − f (x)(x ) T x
x 4 (x )T ∇ 2 f (x)x x 2 + ∇ f (x)T x x 2 − ∇ f (x)T x (x )T x , x 4
(2.5) (2.6)
where f and f denote the derivatives of f (x(u )) with respect to u, the partial derivatives in ∇ f and ∇ 2 f are taken with respect to the variables x, y, z. If we take f as the vector-valued function x, then we have ∇ f = I , ∇ 2 f = 0. Therefore, from (2.3) and (2.6), we have
c x = κ .
(2.7)
It is easy to prove that the differential operators introduced above are geometric. That is, we have the following results. Theorem 2.1. Gradient operator, divergence operator and Laplace operator are independent of the specific parametric representation of the curve.
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
67
2.2.2. Green’s formula For the gradient and Laplace operators, we have the following results. Theorem 2.2. Given a curve C defined on the interval [a, b]. Let f : C → R ∈ C 2 (C), h : C → R ∈ C 1 (C), with at least one of them compactly supported. Then
hc f + ∇c h, ∇s f ds = 0.
(2.8)
C
If both f and h are C 2 , then
[hc f − f c h] ds = 0.
(2.9)
C
Proof. From the definition of the gradient operator, we have
∇c h, ∇s f ds =
C
C
(∇ h)T x (∇ f )T x ds = x 2
b =
h a
b
h (∇ f ) T x
x
a
b du = −
h a
(∇ f )T x x
[(x )T ∇ 2 f x + (∇ f )T x ]x 2 − (∇ f )T x (x )T x du = x 3
Hence (2.8) is proved. The proof of (2.9) is similar. We omit the details.
du
hc f ds.
(2.10)
C
2
Theorem 2.3. Given a space curve C defined on the interval [a, b] and given f ∈ C 2 (C), h ∈ C 1 (C), we have
hc f + ∇c f , ∇c h ds =
2
h(xi ) t(xi ), ∇c f (xi ) .
(2.11)
i =1
C
If both f and h are C 2 , then we have
[hc f − f c h] ds =
2
h(xi ) t(xi ), ∇c f (xi ) − f (xi ) t(xi ), ∇c h(xi ) ,
(2.12)
i =1
C
where t(xi ) are outward unit tangent vectors of C at the end points xi with x1 = x(a) and x2 = x(b). Proof. The proof of (2.11) is similar to the proof of (2.8), but boundary conditions need to be taken into account. In (2.10), the term
h(x(b))(∇ f (x(b))) T x (b)
x (b)
−
h(x(a))(∇ f (x(a))) T x (a)
x (a)
is added. It is easy to derive that
−
h(x(a))(∇ f (x(a))) T x (a)
x (a)
h(x(b))(∇ f (x(b))) T x (b)
x (b)
= h x(a) t x(a) , ∇c f x(a) ,
= h x(b) t x(b) , ∇c f x(b) .
Hence (2.11) is proved. Exchanging h and f in (2.11), we obtain
C
f c h + ∇c f , ∇c h ds =
2
f (xi ) t(xi ), ∇c h(xi ) .
i =1
Subtracting (2.13) from (2.11), we obtain (2.12).
2
(2.13)
68
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
2.3. Global properties of space curves Theorem 2.4. Let C be a closed smooth space curve and L be its arc-length. Then we have
x T κ ds.
L=−
(2.14)
C
Furthermore, if C is a closed smooth plane curve, we have
A=
1
x T N ds,
2
(2.15)
C
where A stands for the area of the region enclosed by the curve, N(x) is the outward normal of the curve at x. Proof. From (2.7) and the Green’s formula, we have
−
x
T
C
T
x T c x(u ), y (u ), z(u )
κ ds =
ds
C
T
∇c x, ∇c x + ∇c y , ∇c y + ∇c z, ∇c z
=
ds
C
+ = C
2
T
x(xi ) t(xi ), ∇c x(xi ) + y (xi ) t(xi ), ∇c y (xi ) + z(xi ) t(xi ), ∇c z(xi )
i =1
(∇ x)T x (x )T ∇ x + (∇ y )T x (x )T ∇ y + (∇ z)T x (x )T ∇ z x(b) T x (b) x(a) T x (a) ds + − 2 x (b) x (a) x
=
ds. C
Hence (2.14) is derived. The following facts are used in the derivation above:
x(a) = x(b),
x (a) = x (b),
t x(a) = −t x(b) .
Note that if the space curve is not closed, the boundary term does not vanish. Hence (2.14) is not true. Equality (2.15) is a well known fact (see Dierckx, 1995, p. 38), we ignore its proof. 2 Remark 2.1. The arc-length formula (2.14) depends on
κ , but is independent of the torsion of the curve.
Corollary 2.1. Let C be a closed smooth space curve. Then for any vector a ∈ R3 , we have
a T κ ds = 0.
(2.16)
C
Furthermore, if C is a closed smooth plane curve, then for any vector a ∈ R2 , we have
a T N ds = 0.
(2.17)
C
Proof. In (2.14), taking x as x + a, and noticing the fact that the translation does not change the arc-length of the curve, we obtain (2.16). (2.17) can be similarly derived from (2.15). 2 Corollary 2.2. Let C be a closed smooth space curve. We have
κ ds = 0.
(2.18)
C
Furthermore, if C is a closed smooth plane curve, we have
N ds = 0. C
(2.19)
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
Proof. In (2.16), taking a = [1, 0, 0] T , [0, 1, 0] T , [0, 0, 1] T , respectively, we obtain (2.18). (2.19) is similarly proved.
69
2
Theorem 2.5. Given a smooth curve motion flow
∂x = V (x), ∂t then for a closed space curve, we have
dL (t )
V (x) T κ (x) ds,
=−
dt
(2.20)
C(t )
and for a closed plane curve, we have
d A (t ) dt
1
=
V (x) T N(x) ds,
2
(2.21)
C(t )
where N(x) is the outward normal of the curve at x. Proof. For plane curves, the theorem is well known (see Gage, 1986, and Dierckx, 1995, p. 41). Now we prove (2.20) for space curves. From
b
L (t ) =
ds =
x du ,
a
C(t )
and Green’s formula, we have
dL (t ) dt
b = a
∂x du = ∂t x (b)
=
T
x (b)
x (b)
=
b a
T
x
∂ x du ∂t
x
T b T ∂ x(b) x (a) ∂ x(a) x ∂x − − du ∂t x (a) ∂t x ∂t a
T
x (b)
V x(b) −
x (a)
T
x (a)
V x(a) −
b a
x
x
T V (x) du .
(2.22)
If the curve considered is smooth and close, then we have
x (b)
T
V x(b) =
x (b)
x (a)
x (a)
T
V x(a) .
(2.23)
From (2.3), we have
x
x
=
x x 2 − x x T x
x 3
= κ (x) x .
Then from (2.22), we obtain (2.20). The theorem is proved.
2
If the curve considered is not closed, the equality (2.23) is not valid in general. Hence, (2.20) may not be true. However, if V (x(a)) and V (x(b)) are in the normal direction of the curve, then the terms related to the boundary vanish. Then, we obtain the following corollary. Corollary 2.3. Given a space curve normal-motion flow
∂x = v (x)n (x), ∂t we have
dL (t ) dt
=− C(t )
v (x)κ 2 (x)
κ (x) +
ds.
(2.24)
70
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
Note that N(x) and n(x) are different. N(x) is the outward normal of the plane curve, n(x) point to the center of the osculating circle of the curve at x. If the curve is convex, then N(x) = −n(x). 2.4. Area changing rate for space curves The concept of area enclosed by a close curve is valid only for a close plane curve. For a space curve, the concept of area is not well defined. In the following, we replace the concept of the enclosed area by the concept of the area of the swept surface of the geometric flow. Given a sufficiently smooth space curve C = {x(u ) ∈ R3 : u ∈ [a, b]}. Without loss of generality, we assume that the curve does not contain straight line segment. If the curve contains straight line segments, we divide the curve into several pieces at the end points of the straight line segments, such that each piece is either a straight line or contains no straight line. If a piece of the curve is a straight line, then we say this curve piece is trivial. Otherwise, we say the curve piece is non-trivial. We consider only the non-trivial curves in the following. Let u 1 < u 2 < · · · < un be all the zero points of κ (u ) of the curve C. Then around ul , there exists a vector nl ∈ R3 , such that κ (u ) and κ (u ) are in the following forms
κ (u ) = (u − ul )αl nl + O (u − ul )αl +1 ,
κ (u ) = |u − ul |αl nl + O |u − ul |αl +1 , where
(2.25) (2.26)
αl 1 is an integer. Then around ul ,
n x( u ) =
α nl κ (u ) = sign(u − ul ) l + O |u − u l | . κ (u ) nl
We define a sign for n(x) in the interval (ul , ul+1 ) as follows
sl =
sl−1 , if αl an even number, −sl−1 , if αl an odd number,
where s0 is defined as 1. Definition 2.1. If C is an open curve or a close curve (x(a) = x(b)) with sn = 1, then we say the curve is orientable. Otherwise, the curve is not orientable. If curve C is orientable, then a sign function as follows is well defined:
sign x(u ) = sl ,
if u ∈ (ul , ul+1 ),
and sign(x(ul )) = 0 for l = 1, . . . , n. Having a sign function sign(x(u )), we define an oriented normal function as
N x( u ) =
sign(x(u ))n(x(u )), n
sl−1 nl , l
κ (u ) = 0, u = u l , l = 1, . . . , n .
Then N(x) is a smooth vector-valued function and
κ (x) = sign(x)κ (x)N(x) = κs N(x), where κs = sign(x)κ (x) is signed curvature. It is a smooth function if the curve x(u ) is sufficiently smooth. Given a curve normal-motion flow for a non-trivial curve
∂x = V (x). ∂t The area of the surface swept by the curve motion is
t V (x) T N(x) ds dτ .
A (t ) = 0 C(τ )
Then the changing rate of A (t ) is
d A (t ) dt
= C(t )
V (x) T N(x) ds.
(2.27)
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
71
Fig. 3.1. Left column: The initial space curve to be evolved. Middle column: The evolved space curve family at the early stage using curvature flow (3.1). Right column: The evolved space curve family at the later stage, where the curves become circular and shrink. The figures in the first row are the space curves observed in the z-direction [0, 0, 1] T . The figures in the second row are the space curves observed in the direction [0.1, 0.9, 0.1] T .
3. Second order curve flows In this section, we propose several second-order geometric flows. Some properties of these flows are discussed. These properties are illustrated by some graphical examples. For ease of illustration, all curves used, except for Fig. 3.1, are planar. The figures presented are not schematic. They have been produced by real computation using the numerical method presented in Xu (2013). 3.1. Curvature flow Let x(u ) be a given space curve with u ∈ [a, b]. Now we define a second order flow. As for the mean-curvature flow for surface, the curvature flow for curve is defined as
∂x = κ n. ∂t
(3.1)
Remark 3.1. The counterpart in the following of the mean-curvature flow ∂∂xt = HN
∂x = κN ∂t is not correct, since the right-hand side depends on the orientation of N. Flow (3.1) can be derived as the gradient flow for minimizing the arc-length functional. For a space curve, from (2.20), we have
dL (t ) dt
κ (x)T κ (x) ds = −
=− C(t )
κ (x) 2 ds 0.
C(t )
Hence, curvature flow is arc-length shortened. For the rate of change of area enclosed by a plane curve with rotation index 1, we have
d A (t ) dt
=
1
T
κ (x)n(x) N(x) ds = −
2 C(t )
where
sign(x)κ (x) ds = −
2 C(t )
sign(x) = −n(x) T N(x) =
1
1
κs (x) ds = −π ,
2
(3.2)
C(t )
1, C is convex at x, −1, C is concave at x.
Here, the convex and concave curves are distinguished when the observer is located in the interior of the closed curve. From (3.2), we know that the curvature flow is shrinking in area. Fig. 3.1 shows an example of the curvature flow, where the left column shows the initial space curve to be evolved. The middle column shows the evolved space curve family at the early stage. The right column shows the evolved space
72
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
Fig. 3.2. Left: The initial curve to be evolved. Middle: The evolved curve family at the early stage using averaged curvature flow (3.4). Right: The steady solution of the averaged curvature flow with the initial curve.
curve family at the later stage, where the curves become circular and shrink. The figures in the first row are the space curves observed in the z-direction [0, 0, 1] T . The figures in the second row are the space curves observed in the direction [0.1, 0.9, 0.1] T . Note that we do not choose [0, 1, 0] T as the view direction, since the curves are almost in the xy-plane at the later stage. 3.2. Averaged curvature flow Let
κ0 = C
κ2 κ +
ds
C
κ ds, κ +
(3.3)
where C = {x(u ) ∈ R3 : u ∈ [a, b]}. Then the counterpart of the averaged mean curvature flow is defined as
∂x = (κ − κ0 )n . ∂t
(3.4)
Since
C(t )
κ (κ02 − κ0 κ ) ds = κ02 κ +
κ
C(t )
ds − κ0
κ +
C(t )
κ2 κ +
ds = 0,
we have
dL (t ) dt
=− C(t )
κ (κ 2 − κκ0 ) ds = − κ +
C(t )
κ (κ 2 − κκ0 ) κ (κ02 − κ0 κ ) + ds = − κ + κ +
C(t )
κ (κ − κ0 )2 ds 0. κ +
Therefore, the averaged curvature flow is arc-length shortened for a space curve. For plane curves, from (2.21), we have
d A (t ) dt
=
1
2 C(t )
κ (x)(κ (x) − κ0 (t ))n(x)T N(x) 1 ds = − κ + 2
sign(x) C(t )
κ 2 (x) − κ (x)κ0 ds κ +
= 0 if sign(x) = 1. If the curve is convex, the flow is area preserving. For a space curve, if κ > 0 everywhere, then from (2.27), we know the flow is area preserving. Fig. 3.2 shows an example of the averaged curvature flow, where an initial non-convex curve is evolved. Different from the curvature flow, the evolved curves are not shrinking. It is natural to believe that the steady solution of the flow is a circle with radius 1/κ0 . However, this is not true: the last figure shows that the steady solution of the flow is a curve consisting of piecewise circular arcs with radius 1/κ0 . The computational result shows that the total squared error C(t ) (κ − κ0 )2 ds approaches to zero. Hence, flow (3.4) is different from flow (1.1), even for the case of a plane curve. However, if = 0 and the curve considered is a close and convex plane curve with winding number 1, then from (3.3), we have
κ0 =
κ ds C
ds = C
1
κs ds =
L
2π L
,
C
and therefore, flow (3.4) is same as the flow (1.1).
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
73
Fig. 3.3. Left: The initial curve to be evolved. Middle: The evolved curve family at the early stage using averaged curvature flow (3.5). Right: The evolved curve family at the later stage.
3.3. Area-preserved averaged curvature flow Flow (3.4) is area-preserved for the curve with κ > 0. If the curvature has zero points, the area preserving property is not valid. Now we define an area-preserved averaged curvature flow. Let
1
κ0 =
κs ds.
L C
Then the flow
∂x = (κs − κ0 )N ∂t
(3.5)
is arc-length shortening and area preserving. In fact, from
κ02 ds − C(t )
κ0 κs ds = 0,
C(t )
and (2.20), we have
dL (t ) dt
=−
κs 2 − κs κ0 ds = −
C(t )
κs 2 − κs κ0 + κ02 − κ0 κs ds = −
C(t )
(κs − κ0 )2 ds 0.
C(t )
Therefore, the averaged curvature flow is arc-length shortening for a space curve. Furthermore, from (2.27), we have
d A (t )
=
dt
1
2
κs (x) − κ0 (t ) N(x)T N(x) ds = 0.
C(t )
Hence, the flow is area preserving. Fig. 3.3 shows an example of the area-preserved averaged curvature flow, where the left figure shows an initial nonconvex curve to be evolved. The middle figure shows the evolved curve family at the early stage. The right figure shows the evolved curve family at the later stage. The curve approaches a circle with radius 1/κ0 . 3.4. Arc-length preserved averaged curvature flow Let
κ1 = C
κ3 κ +
ds
C
κ2 κ +
ds.
Then we introduce another averaged curvature flow as follows
∂x = (κ − κ1 )n . ∂t Since
dL (t ) dt
=− C(t )
κ 3 − κ 2 κ1 ds = 0. κ +
(3.6)
74
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
Fig. 3.4. Left: The initial curve to be evolved. Middle: The evolved curve family at the early stage of the evolution using the averaged curvature flow (3.6). Right: The evolved curve family at the later stage.
Therefore, this averaged curvature flow is arc-length preserved for evolving space curves. If C(t ) is a convex plane curve, flow (3.6) is similar to (1.2) proposed by Ma and Zhu (2012). If C(t ) is not convex, then (3.6) is not similar to (1.2). Fig. 3.4 shows an example of this averaged curvature flow, where an initial non-convex curve is evolved. Different from the curvature flow, the evolved curves are not shrinking. Again, the steady solution of the flow is a curve consisting of piecewise circular arcs with radius 1/κ1 . The computational result shows that the total squared error C(t ) (κ − κ1 )2 ds approaches to zero. 4. Fourth order curve flows In this section, several fourth-order flows, including curve diffusion flow, quasi-curve diffusion flow, projective quasicurve diffusion flow, minimal curvature flow and minimal squared-curvature flow, are presented. Weak forms of some these flows are also given. The illustrations provided are generated using the numerical method presented in Xu (2013). 4.1. Curve diffusion flow To define a curve diffusion flow which is similar to surface diffusion flow, let us recall the surface diffusion flow (see Mullins, 1957):
∂x = −2s Hn, ∂t
(4.1)
where s , H and n denote the Laplace–Beltrami operator, mean curvature and surface normal, respectively. Note that mean curvature H depends on the orientation of the surface normal. Namely, if we take the opposite sign of the normal, then the mean curvature will change sign. Hence, the right-hand side of (4.1) is independent of the orientation of the surface normal. If we take directly the form of the surface diffusion flows as the curve diffusion flow, we obtain the following form for curve diffusion flow:
∂x = −c κ n, ∂t
(4.2)
where c and n are the Laplace operator and principal normal of the curve. However, no matter how smooth the curve x(u ) is, κ (u ) can be C 0 smooth (see (2.26)) around zero points of κ , and c κ is not well defined. Therefore, curve diffusion flow (4.2) is not well defined at the points where κ = 0. However, the flow
∂x = −c κs N ∂t
(4.3)
is well defined. From (2.20), (2.27) and Green’s formulas, we have
dL (t ) dt
c κs NT κ ds =
= C(t )
and
d A (t ) dt
c κs κs ds = −
C(t )
=− C(t )
∇c κ , ∇c κ ds 0,
C(t )
c κ NT N ds = − C(t )
c κ ds =
∇c κ , ∇c 1 ds = 0.
C(t )
Hence, the flow is arc-length shortening and area-preserving. Fig. 4.1 shows an example of the curve diffusion flow, where an initial un-convex curve is evolved. The middle figure shows the evolved curve family at the early stage using the curve diffusion flow (4.3). The right figure shows the evolved curve family at the later stage.
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
75
Fig. 4.1. Left: The initial curve to be evolved. Middle: The evolved curve family at the early stage using curve diffusion flow (4.3). Right: The evolved curve family at the later stage using the curve diffusion flow (4.3).
Fig. 4.2. Left: The initial curve to be evolved using quasi-curve diffusion flow (4.4). Middle: The evolved curve family at the early stage of the evolution. Right: The evolved curve family at the later stage of the evolution.
4.2. Quasi-curve diffusion flow For surface, we had introduced the following quasi-surface flow (see Xu, 2008)
∂x = −2s x. ∂t The flow has area shrinking property. Similar to this, we introduce the following quasi-curve diffusion flow
∂x = −c κ . ∂t
(4.4)
From (2.3), we can see that, for a sufficiently smooth curve, c κ is well defined though c κ may not be well defined. From (2.20) and Green’s formula, we have
dL (t ) dt
=
c κκ ds = −
C(t )
∇c κ : ∇c κ ds 0.
C(t )
Hence, flow (4.4) is arc-length shortening for a close space curve. From Frenet–Serret formulas, we can derive that, if
κ > 0,
c κ = κ − τ 2 − κ 2 n − 3κκ t + κτ b = c κ − τ 2 − κ 2 n − 3κ ∇c κ + κτ b. Hence, the quasi-curve diffusion flow contains the motions in the tangent and binormal directions. Note that torsion is involved in this equation. Using Green’s formula (2.12), the weak form of the quasi-curve diffusion flow is written as
C(t )
∂x φ ds = − ∂t
c xc φ ds,
∀φ ∈ C 02 (Ω),
(4.5)
C(t )
where C 02 (Ω) consists of C 2 smooth functions over Ω with support in the interior of Ω . Fig. 4.2 shows an example of the quasi-curve diffusion flow, where an initial non-convex curve is evolved. The middle figure shows the evolved curve family at the early stage using the quasi-curve diffusion flow. The right-hand figure shows the evolved curve family at the later stage. The curve will shrink with further evolution.
76
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
Fig. 4.3. Left: The initial curve to be evolved using the projective curve diffusion flow (4.6). Middle: The evolved curve family at the early stage of the evolution. Right: The evolved curve family at the later stage.
4.3. Projective quasi-diffusion flow As we have mentioned before, the quasi-curve diffusion flow contains the motions in the tangent and binormal directions. To disable the tangential motion, we construct the following form flow:
∂x = − I − ttT c κ , ∂t
(4.6)
where t is the unit tangent vector of the curve. Since κ does not depend on the orientation of the normal, the right-hand side of (4.6) is independent of the orientation of the normal. The right-hand side of (4.6) is the projection of c2 x in the normal plane, we therefore name it as projective quasi-diffusion flow. From (2.20) and Green’s formula, we have
dL (t ) dt
(c κ )T I − ttT κ ds =
= C(t )
(c κ )T κ ds = −
C(t )
∇c κ : ∇c κ ds 0,
C(t )
where A : B denotes trace of A B . Hence, flow (4.6) is arc-length shortening for evolving close space curves. Let y = κ . Then, using the Green’s formula, we can write the weak-form of (4.6) as: Find (x, y) such that T
⎧ ∂x ⎪ ⎪ φ ds = − I − tt T c yφ ds, φ ∈ C 02 (Ω), ⎪ ⎪ ∂ t ⎪ ⎨ C(t ) C(t ) ⎪ ⎪ ⎪ yψ ds = − ∇c x∇c ψ ds, ψ ∈ C 02 (Ω). ⎪ ⎪ ⎩ C(t )
(4.7)
C(t )
Fig. 4.3 shows an example of the projective curve diffusion flow (4.6), where an initial non-convex curve is evolved. The middle figure of Fig. 4.3 shows a curve family at the early stage obtained by this flow. The right-hand figure shows a curve family at the later stage. The curve will shrink continuously for further evolution. 4.4. Minimal squared-curvature flow For surfaces, Willmore flow (see Willmore, 1993) is derived from minimizing the total squared mean-curvature. For curves in R3 , we minimize the following energy
κ 2 ds.
E (C) = C
To compute its first order variation, we represent E as
E (C) =
b
2
2
κ ds = C
c x ds =
c x2 x du .
(4.8)
a
C
Assume C is an extreme point of the energy functional (4.8). Let
C(Φ, ε ) = x(u , ε ) = x(u ) + ε Φ(u ): u ∈ Ω ,
|ε | 1,
be a perturbation of the curve C, where Φ is smooth vector-valued function with compact support. Then
δ E (C), Φ =
d dε
E C(Φ, ε )
ε =0
= 0,
(4.9)
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
77
where
b
δ E (S), Φ =
2
δ c x x
b du =
a
2 c x, δ(c x) x + κ 2 δ x du .
(4.10)
a
It is not difficult to derive that
(x )T Φ δ x = , x Φ Φ (x )T x + x (Φ )T x + x (x )T Φ 2x (x )T x (x )T Φ 2x (x ) T Φ − + δ(c x) = 2 − 4 x x x 4 x 6 2x (x ) T Φ x (Φ ) T x + x (x ) T Φ 2x (x ) T x (x ) T Φ = c Φ − − + . 4 4 x x x 6 Notice that c x, x = 0. From (4.10) we have
δ E (S), Φ =
2c x, c Φ −
2x (x ) T Φ
x 4
C
=
−
2c x, c Φ ds +
b
=
−
2c x, c Φ ds + a
C
ds +
κ 2 (x )T Φ x
4x (x ) T c x
x 3
a
C
b a
b
3κ 2 x
x
+
κ 2 x x
du
, Φ du
, Φ du .
(4.11)
Using Green’s formulas, we have
δ E (S), Φ = 2
2 c x, Φ ds +
C
1
3κ 2 x
x
C
x
, Φ ds.
Since Φ is an arbitrary vector-valued function, we obtain, from (4.9), the following Euler Lagrange equation
1
2c2 x +
x
3κ 2 x
x
= 0.
To simplify the derivation, but without loss of generality, we assume the parameter of the curve is the arc-length. Then, we have
1
x
3κ 2 x
x
= 3 κ 2 x = 3κ 2 κ + 3 κ 2 x .
Since term 3(κ 2 ) x is in the tangential direction, we then construct the following geometric flow by ignoring the tangential term:
∂x = −2c κ − 3κ 2 κ . ∂t
(4.12)
We name this geometric flow as minimal squared-curvature flow. From (2.20) and Green’s formula, we have
dL (t ) dt
=
2c κκ + 3κ 4 ds =
C(t )
3κ 4 − 2∇c κ : ∇c κ ds = 3
C(t )
C(t )
κ 4 ds − 2
2 x ds,
C(t )
where the third-order derivative x is taken with respect to the arc-length parameter. If a space curve is bumpy, then the term 2 C(t ) x 2 ds is dominant. The flow is arc-length shortening. If a space curve is not bumpy, the term 3 C(t ) κ 4 may be dominant. The flow is arc-length lengthening. For instance, if the curve C is a circle, then we can derive that
C(t )
3κ 4 − 2∇c κ : ∇c κ ds =
κ 4 ds > 0.
C(t )
Hence, for a circle, flow (4.12) is arc-length lengthening.
78
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
Fig. 4.4. Left: The evolved curve family at the early stage using minimal squared-curvature flow (4.12) for the initial curve as shown in Fig. 4.2. Middle: The evolution at a stage where the curve becomes a circle. Right: The evolved curve family at the later stage. The circles become larger and larger.
Using (4.11) and noting that ∇c φ = ∇c x∇c φ , we obtain the weak-form of the minimal squared-curvature flow as
∂x φ ds = − ∂t
C(t )
2c xc φ − 3κ 2 ∇c x∇c φ ds,
∀φ ∈ C 02 (Ω).
(4.13)
C(t )
Fig. 4.4 shows an example of the minimal squared-curvature flow, where an initial non-convex curve, as shown in Fig. 4.2, is evolved. The left-hand figure shows an evolved curve family at the early stage of the evolution. It is easy to see that the curve is shortening. The middle figure shows a stage where the curve becomes a circle. The right-hand figure shows an evolved curve family at the later stage. The curve becomes a circle, and further evolution makes the circle larger and larger. 4.5. Minimal curvature flow As we have seen, the minimal squared-curvature flow is arc-length lengthening for a circle. To obtain a flow that is not arc-length lengthening for a circle, we minimize the following energy:
E (C) =
κ ds.
(4.14)
C
Assume C is an extreme point of the energy functional (4.8). Let
C(Φ, ε ) = x(u , ε ) = x(u ) + ε Φ(u ): u ∈ Ω ,
|ε | 1,
be a perturbation of the curve C, where Φ is smooth vector-valued function with compact support. Then
d δ E (C), Φ = E C(Φ, ε ) dε
ε =0
= 0,
(4.15)
where
b
δ E (S), Φ =
δ κ x du .
(4.16)
a
Without loss of generality, we may assume u is arc-length. Then it is not difficult to derive that
δ κ x = δ
x 2 x 2 − (x T x )2 x 2
= n, Φ − x x , Φ .
Assume n = 0. Then using Green’s formulas, we have
δ E (S), Φ =
C
c n + κκ + x κ , Φ ds.
(4.17)
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
79
Fig. 4.5. The left-hand figure shows the initial curve to be evolved. The remaining figures show the evolved curves using minimal curvature flow (4.18) at different stages of the evolution.
Since Φ is an arbitrary vector-valued function, then from (4.15) we obtain the following Euler Lagrange equation
c n + κκ + x κ = 0. Since term x κ is in the tangential direction, we then construct the following geometric flow by ignoring the tangential term
∂x = −[c n + κκ ]. ∂t
(4.18)
We name this flow as minimal curvature flow. From (2.21) and Green’s formula, we have
d A (t ) dt
=
1
2
κ 2 − ∇c n : ∇c n ds.
C(t )
Suppose the curve C is a plane convex curve, then from the Frenet–Serret formula we can derive that
κ 2 − ∇c n : ∇c n ds = 0.
C(t )
Hence, for a plane convex curve, flow (4.18) is area preserving. In fact, if C(t ) is convex plane curve, then C(t ) κ ds = 2π is a constant. This means that C(t ) is a minimizer of the energy functional (4.14). From the Frenet–Serret formula we can derive that, if κ > 0, (4.18) becomes
∂x = −τ b + tκ + τ 2 n, ∂t where the derivative is taken with respect to the arc-length s. From (4.17), the weak-form of the flow (4.18) can be written as
C(t )
∂x φ ds = − ∂t
c xc φ
κ
− κ ∇c x∇c φ ds,
∀φ ∈ C 02 (Ω).
(4.19)
C(t )
The right-hand side of (4.19) has a singularity when κ = 0, which reflects the fact that the normal n is not continuous. To 1 overcome this difficulty, we replace the term κ1 with κ + , where > 0 is a small number. Fig. 4.5 shows an example of the minimal curvature flow, where an initial non-convex curve is evolved. Note that the resulting curve has a sharp feature. This is understandable since the energy (4.14) is, in fact, the L 1 norm of the curvature. Remark 4.1. Note that (4.18) is derived under the assumption that n = 0. If the curve C has points where n = 0 (x and x are not linearly independent), this equation cannot be used. However, the weak form (4.19) is valid. 5. Conclusion and future work Several second-order and fourth-order geometric flows have been constructed for space curves. Some properties of the changing rates of the arc-length and areas have been discussed. These properties show that the geometric flows constructed for space curves are different in nature from their counterparts of the plane curve case and surface case. Short-term and long-term behaviors of the evolved curves using the constructed flows have been illustrated by graphical examples. These behaviors show that each of the flows has its own distinct features. As part of future work, we plan to study the applications of the geometric flows constructed. A direct application may be space curve denoising. All the flows with arc-length shortening property can be used for space curve denoising. However,
80
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
if a flow is both arc-length shortening and area-preserving, then it would be more desirable, since it will better preserve the shape of the curve. Another application can be image restoration and image denoising. For instance, the 2D version of flow (4.5) can be used to accomplish the task with sharp feature of the objects in the image to be preserved. A third application area of the geometric flows can be geometric design. In many shape design problems, shapes are designed from a wireframe consisting of space curves. Each of the smooth space curves, say B-spline curves, is constructed from interpolating a few fixed control points. The remaining degrees of freedom (DOF) of the curve need to be determined. These DOF can be determined by the geometric flows constructed. The space curves considered in this paper are expressed in parametric forms. If the curves are in the level-set form (intersections of two implicit form surfaces), how to reformulate these flows is an interesting problem that is worth further study. References Altschuler, S.J., Grayson, M.A., 1992. Shortening space curves and flow through singularities. J. Differ. Geom. 35 (2), 283–298. Andrews, B., 1998. Evolving convex curves. Calc. Var. 7, 315–371. Chen, S., Lin, M., Wu, J., 2007. A weighted curvature flow for planar curves. Int. J. Math. Comput. Simul. 1 (3), 275–278. Chou, K., 2001. The curve shortening flow: The classical approach. In: RIMS Kokyuroku, vol. 1210. Cohen, E., Tekumalla, L., 2005. Smoothing space curves with the MLS projection. In: Geometric Modeling, Vizualization and Graphics. Salt Lake City. Dierckx, P., 1995. Curve and Surface Fitting with Splines. Oxford University Press. Do Carmo, M.P., 1976. Differential Geometry of Curves and Surfaces. Prentice Hall, Englewood Cliffs, NJ. Dolcetta, I.C., Vita, S.F., March, R., 2002. Area preserving curve shortening flows: From phase transitions to image processing. Interfaces Free Bound. 4 (4), 325–343. Dziuk, G., 1999. Discrete anisotropic curve shortening flow. SIAM J. Numer. Anal. 36 (6), 1808–1830. Gage, M., 1983. An isoperimetric inequality with application to curve shortening. Duke Math. J. 50, 1225–1229. Gage, M., 1984. Curve shortening makes convex curves circular. Invent. Math. 76, 357–364. Gage, M., 1986. On an area-preserving evolution equation for plane curves. In: Contemp. Math., vol. 51, pp. 51–62. Gage, M., 1990. Curve shortening on surfaces. Ann. Sci. Éc. Norm. Super. 23, 229–256. Gage, M., Hamilton, R., 1986. The heat equation shrinking convex plane curves. J. Differ. Geom. 23, 69–96. Grayson, M., 1987. The heat equation shrinks embedded plane curves to points. J. Differ. Geom. 26 (2), 285–314. Grayson, M., 1989. Shortening embedded curves. Ann. Math. 129, 71–111. Guziec, A., Ayache, N., 1994. Smoothing and matching of 3-d space curves. Int. J. Comput. Vis. 12 (1), 79–104. Hildebrandt, K., Polthier, K., Preuss, E., 2005. Evolution of 3d curves under strict spatial constraints. In: Proceedings of the Ninth International Conference on Computer Aided Design and Computer Graphics. CAD-CG’05. IEEE Computer Society, Washington, DC, USA, pp. 40–45. Leng, J., Zhang, Y., Xu, G., 2011. A novel geometric flow-driven approach for quality improvement of segmented tetrahedral meshes. In: Proceedings of the 20th International Meshing Roundtable, pp. 347–364. Liu, Y., Lai, Y., Zang, Y., Hu, S., 2008. Fairing wireframes in industrial surface design. In: 2008 IEEE International Conference on Shape Modeling and Applications. Stony Brook, New York. ACM, pp. 29–35. Ma, L., Chen, D., 2007. Curve shortening in a riemannian manifold. Ann. Mat. Pura Appl. 186, 663–684. Ma, L., Zhu, A., 2012. On a length preserving curve flow. Monatshefte Math. 165 (1), 57–78. Mullins, W.W., 1957. Theory of thermal grooving. J. Appl. Phys. 28, 333–339. Pan, S., Yang, J., 2008. On a non-local perimeter-preserving curve evolution problem for convex plane curves. Manuscr. Math. 127, 469–484. Sapiro, G., 2001. Geometric Partial Differential Equations and Image Analysis. Cambridge University Press, Cambridge. Sapiro, G., Tannenbaum, A., 1995. Area and length preserving geometric invariant scale-spaces. IEEE Trans. Pattern Anal. Mach. Intell. 17, 67–72. Spira, A., Kimmel, R., 2007. Geometric curve flows on parametric manifolds. Chin. J. Comput. Phys. 223 (1), 235–249. Willmore, T.J., 1993. Riemannian Geometry. Clarendon Press, Oxford, England. Xu, G., 2008. Geometric Partial Differential Equation Methods in Computational Geometry. Science Press, Beijing, China. Xu, G., 2013. Numerical integration of geometric flows for space curves. Research Report No. ICM-13-01. Institute for Computational Mathematics and Scientific/Engineering Computing.
Contents lists available at ScienceDirect
Computer Aided Geometric Design www.elsevier.com/locate/cagd
Construction of several second- and fourth-order geometric partial differential equations for space curves ✩,✩✩ Guoliang Xu a , Xuyang Yang b a
State Key Laboratory of Scientific and Engineering Computing, Institute of Computational Mathematics, Academy of Mathematics and System Sciences, Chinese Academy of Sciences, Beijing 100190, China China Shipbuilding Industry Corporation, 750 Proving Ground, Kunming 650051, China
b
a r t i c l e
i n f o
Article history: Received 29 December 2012 Received in revised form 12 December 2013 Accepted 26 December 2013 Available online 7 January 2014 Keywords: Geometric partial differential equations Geometric flows Space curves
a b s t r a c t Geometric partial differential equations for curves and surfaces are used in many fields, such as computational geometry, image processing and computer graphics. In this paper, a few differential operators defined on space curves are introduced. Based on these operators, several second-order and fourth-order geometric flows for evolving space curves are constructed. Some properties of the changing rates of the arc-length of the evolved curves and areas swept by the curves are discussed. Short-term and long-term behaviors of the evolved curves are illustrated. © 2014 Elsevier B.V. All rights reserved.
1. Introduction In many design and construction problems, such as reverse engineering, multi-component mesh optimization and industrial geometric design from a wireframe, smooth space curves need to be generated from noisy input data. This goal can be efficiently achieved by utilizing geometric flows to smooth the input data with shape preserved. The aim of this paper is to construct such geometric flows for space curves embedded in R3 . For surfaces embedded in R3 , many geometric flows have been constructed (see Xu, 2008 for references). These include second-order geometric flows (mean-curvature flow, averaged mean-curvature flow, weighted mean-curvature flow, Gaussian curvature flow, etc.), fourth-order flows (surface diffusion flow, quasi-surface flow, Willmore flow, etc.) and sixth-order flows (minimal variation of mean-curvature flow, quasi-minimal variation of mean-curvature flow, etc.). These flows are frequently used in surface modeling, surface denoising, image processing and computer graphics. As one can imagine, all the plane curves in the xy-plane can be extended to the xyz-space by setting the z-direction as constants, obtaining three-dimensional (3D) cylindrical-type surfaces. Then, the mean curvature H of this surface is half of the signed curvature of the plane curve, and the Gaussian curvature K is zero. Then, all the geometric flows for surfaces have their counterparts for plane curves. For this reason, one may think that there is no need to study the construction of geometric flows for curves. However, for space curves in R3 , reducing the surface flows to curves is not always valid, as explained below. Since most of the curvature flows for surfaces embedded in R3 describe the surface motion in the normal direction, we take the normal as an example to explain the difference. For a single-sided surface, the surface normal can be orientated, but the normal of the space curve cannot be orientated. For instance, a straight line in 3D, as a counterpart of the plane, cannot be orientated. Furthermore, the motion of the surface in the normal direction (the tangential motion does not change ✩
This paper has been recommended for acceptance by Konrad Polthier. Supported in part by NSFC under the grants 1101401 and 81173663, NSFC key project under the grant 10990013 and Funds for Creative Research Groups of China (grant No. 11021101). ✩✩
0167-8396/$ – see front matter © 2014 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cagd.2013.12.005
64
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
the surface shape) has one degree of freedom. However, for a space curve, the motion has two degrees of freedom. Hence, the situation for the space curve is rather different from that for the plane curve. In this paper, we consider the construction of geometric flows for space curves. 1.1. Previous work on geometric flows for plane curves There has been much work on plane curve evolution (see Gage, 1983, 1984; Gage and Hamilton, 1986; Ma and Chen, 2007; Sapiro, 2001 for references). We only review a few of them. Most of the curve evolutions focus on the following curve-shortening flow (CSF)
∂x = κs N ∂t and its various generalizations. Here, κs is the signed curvature and N is the unit inward normal of a close plane curve x. CSF is also known as curvature flow or geometric heat flow, which is the counterpart of the mean-curvature flow for surfaces. In Chou (2001), some results on CSF are presented, including the properties of CSF, the anisotropic generalized CSF and the CSF in Klein geometry. In Andrews (1998), the behavior of convex curves driven by the curvature flow is considered. The long-term behavior of solutions and properties of limiting shapes are described. As shown in Gage and Hamilton (1986) and Grayson (1987), the curvature flow will shrink an embedding planar curve to a point, becoming round in the limit without developing singularities. In Dziuk (1999), an anisotropic CSF is studied, and a numerical scheme for this flow is developed. Optimal error estimates for a semi-discrete scheme are provided. A weighted curvature flow for planar curves, a generalization of the curvature flow, is discussed in Chen et al. (2007). The convexity of the curve driven by the weighted curvature flow is studied. The geometric curve evolution equations in a plane have been extended to surfaces and manifolds in Grayson (1989), Gage (1990), Spira and Kimmel (2007). Efficient numerical schemes for the implementation of the classical geodesic curve evolution equations on parametric manifolds are presented in Spira and Kimmel (2007). These flows are used on the scale space of images painted on manifolds and segmentation by an active contour model. Grayson (1989) and Gage (1990) generalized the study of curvature flow of closed curves in a plane to that in surfaces, and analyzed the long-term behavior of smooth curves immersed in a surface. 1.2. Length- or area-preserved flows for plane curves It is well known that geometric heat flow is the fastest way to shrink the length of a curve (Gage, 1986). A planar embedded curve that evolves according to this flow becomes more and more circular as t increases, and eventually shrinks to a point in finite time (Gage and Hamilton, 1986; Grayson, 1989). Such a shrinking effect is an undesirable feature in many applications, such as curve denoising, since the evolving curve will drift rapidly from the desired shape. In order to avoid the shrinking effect, Gage (1986) proposed two variations of the geometric heat flow. The first one is defined by
∂x 2π = κs − N, ∂t L
(1.1)
where L is the arc-length of the curve. The second one is defined by
∂x π x, N = κs + N + α T, ∂t A
where A is the area enclosed by the curve and α T is a tangential component that does not change the shape. The resulting flows are curve-shortening motions that preserve the area of the region enclosed by the curve. In Sapiro and Tannenbaum (1995), area-preserving multi-scale representations of planar curves are described. The representations are obtained by deforming the curve via geometric heat flows while simultaneously magnifying the plane for keeping the enclosed area constant. In Dolcetta et al. (2002), several models in phase separation theory are used to generate planar curve evolution without shrinkage. The theoretical validation of a level-set method to these models is discussed. In Ma and Zhu (2012), a length preserving curve flow
∂x 1 = κs − α (t ) N with α (t ) = ∂t 2π
κs2 ds
(1.2)
for convex curves in the plane is considered. It was shown that the area of the region bounded by the evolving curve increases, and the evolving curve converges to a circle. Pan and Yang (2008) proposed another length-preserving curve flow for convex curves in the plane of the form
∂x = ∂t
L 2π
− κs−1 N.
They proved that the convex plane curve will become more and more circular and converges to a circle.
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
65
1.3. Previous work on space curves Space curve smoothing or denoising is often required in reverse engineering, CAD modeling, animation, computer vision, image registration and mesh optimization (Leng et al., 2011). To smooth a discrete space curve (a sequence of points) with noise, B-splines were used in Guziec and Ayache (1994) to approximate the discrete data. In Dierckx (1995), noisy scattered data were smoothed by minimizing a functional that combines a term that measures the deviation of the data points from the fit with another term that measures the non-smoothness of the fit. An efficient fairing technique for a wireframe defined by two sets of space curves was presented in Liu et al. (2008). In the fairing process, the parameters are fixed and then released by an iterative gradient descent optimization method. In Cohen and Tekumalla (2005), a 3D curve-smoothing algorithm using the moving least squares (MLS) projection is proposed with the primary focus on smoothing boundaries of point clouds obtained during reverse engineering of laser scanned models. An algorithm for fairing of space curves with respect to spatial constraints based on a vector-valued curvature function is presented in Hildebrandt et al. (2005). It describes a curve-smoothing flow that satisfies strict spatial constraints and allows simultaneous control of both curvature functions. In Altschuler and Grayson (1992), the authors show how to use the curvature flow on space to define a natural continuation of the planar solution for the whole of the time. 1.4. Main contributions of this paper Our attention is focused on constructing geometric flows for space curves. First, gradient, divergence and Laplace operators on space curves are introduced. Then several second-order geometric flows and fourth-order geometric flows for space curves are constructed. We consider second- and fourth-order flows for the following reasons. First, second-order flows are easier to solve and can be used to denoise efficiently. However, the second-order flows can only achieve G 0 continuity for joining space curves. To achieve G 1 continuity, fourth-order flows are required. Second, for smoothing space curves with shape preservation, fourth-order flows are more desirable since the solution curves of the higher order flows are “softer” than those of the lower order flows. Some properties of the changing rates of the arc-length and areas are discussed. To describe the area changing rate for an evolving space curve, a notion of the area swept surface of the evolving curve is introduced. 1.5. Organization of the paper The remainder of this paper is organized as follows. In Section 2, we introduce some of the notation used and several differential operators on space curves, including the definitions of curvature, normal, gradient, divergence and Laplace operators. Some basic results for these operators are also presented. In Sections 3 and 4, we derive several geometric flows, including second-order and fourth-order flows. Some properties of these flows are discussed. Finally, we conclude the paper with a summary in Section 5. 2. Preliminary In this section, we introduce some of the notation used and several differential operators on space curves, including the definitions of curvature, normal, gradient, divergence and Laplace operators. Some basic results for these operators are also presented. 2.1. Basic differential geometry on space curves 2.1.1. Curvature and normal For a given smooth space curve C defined by
T
C = x(u ) = x(u ), y (u ), z(u ) let
: u ∈ Ω = [a, b] ,
κ be the curvature, which is given as (see Do Carmo, 1976, p. 25) x (u ) × x (u ) x (u )2 x (u )2 − (x (u )T x (u ))2 κ (u ) = = , x (u )3 x (u )3
(2.1)
where x (u ) and x (u ) denote the first- and second-order derivatives of x(u ) with respect to u. Note that κ 0. This is different from the mean curvature of a surface, which can be either positive or negative. The principal normal is given as
n(u ) =
(x (u ) × x (u )) × x (u ) x (u )x (u )2 − x (u )x (u ) T x (u ) . = (x (u ) × x (u )) × x (u ) x (u ) x (u )2 x (u )2 − (x (u )T x (u ))2
(2.2)
It should be noted that n(u ) is well defined for all u such that x (u ) and x (u ) are linearly independent. If x (u ) and x (u ) are linearly dependent, n(u ) is arbitrarily defined. We artificially define it as a zero vector. The curvature κ = 0 is
66
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
well defined. From (2.2), we can derive that n(u ) does not depend on the re-parameterization of the curve. Hence, n(u ) is uniquely defined. In addition to the curvature and normal, we define the curvature vector as
κ (u ) =
x (u )x (u )2 − x (u )x (u ) T x (u )
x (u )4
(2.3)
.
It is easy to see that
κ (u ) = κ (u ),
n(u ) =
κ (u ) if κ (u ) > 0. κ (u )
Hence, n(u ) is not continuous at the point
n (u ) = where
κ (u ) = 0. To obtain a continuous normal function, we define
κ (u ) , κ (u ) +
is a small positive number. In the following, all the stand for the same number.
2.2. Differential operators on space curves 2.2.1. Laplace operator Suppose f is a C 2 smooth function defined on the curve x(u ), then we define Laplace operator c (here the subscript c stands for that the operator is defined on curve) applying to f as
c f = divc (∇c f ),
(2.4)
where the gradient and divergence operators on a space curve are defined as
∇c f =
x f
x 2
divc v =
1
= d
x (x ) T ∇ f
x 2 T v x
x du x
,
=
[(v )T x + vT x ]x 2 − (x )T x vT x . x 4
It is easy to see that ∇c f is the projection vector of ∇ f on the tangent space of the curve. If u is arc-length and v is a tangent vector field, then
T ∇c f = x x ∇ f ,
T
divc v = v
x .
T
If we take f as the vector value function x, we have ∇c x = x(xx )2 . It is easy to see that for any tangent vector v, we have
(∇c x)v = v. Hence, ∇c x acts as an identity operator on the tangent space of the curve. From (2.4) we can derive that
c f (x) = =
f (x)x 2 − f (x)(x ) T x
x 4 (x )T ∇ 2 f (x)x x 2 + ∇ f (x)T x x 2 − ∇ f (x)T x (x )T x , x 4
(2.5) (2.6)
where f and f denote the derivatives of f (x(u )) with respect to u, the partial derivatives in ∇ f and ∇ 2 f are taken with respect to the variables x, y, z. If we take f as the vector-valued function x, then we have ∇ f = I , ∇ 2 f = 0. Therefore, from (2.3) and (2.6), we have
c x = κ .
(2.7)
It is easy to prove that the differential operators introduced above are geometric. That is, we have the following results. Theorem 2.1. Gradient operator, divergence operator and Laplace operator are independent of the specific parametric representation of the curve.
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
67
2.2.2. Green’s formula For the gradient and Laplace operators, we have the following results. Theorem 2.2. Given a curve C defined on the interval [a, b]. Let f : C → R ∈ C 2 (C), h : C → R ∈ C 1 (C), with at least one of them compactly supported. Then
hc f + ∇c h, ∇s f ds = 0.
(2.8)
C
If both f and h are C 2 , then
[hc f − f c h] ds = 0.
(2.9)
C
Proof. From the definition of the gradient operator, we have
∇c h, ∇s f ds =
C
C
(∇ h)T x (∇ f )T x ds = x 2
b =
h a
b
h (∇ f ) T x
x
a
b du = −
h a
(∇ f )T x x
[(x )T ∇ 2 f x + (∇ f )T x ]x 2 − (∇ f )T x (x )T x du = x 3
Hence (2.8) is proved. The proof of (2.9) is similar. We omit the details.
du
hc f ds.
(2.10)
C
2
Theorem 2.3. Given a space curve C defined on the interval [a, b] and given f ∈ C 2 (C), h ∈ C 1 (C), we have
hc f + ∇c f , ∇c h ds =
2
h(xi ) t(xi ), ∇c f (xi ) .
(2.11)
i =1
C
If both f and h are C 2 , then we have
[hc f − f c h] ds =
2
h(xi ) t(xi ), ∇c f (xi ) − f (xi ) t(xi ), ∇c h(xi ) ,
(2.12)
i =1
C
where t(xi ) are outward unit tangent vectors of C at the end points xi with x1 = x(a) and x2 = x(b). Proof. The proof of (2.11) is similar to the proof of (2.8), but boundary conditions need to be taken into account. In (2.10), the term
h(x(b))(∇ f (x(b))) T x (b)
x (b)
−
h(x(a))(∇ f (x(a))) T x (a)
x (a)
is added. It is easy to derive that
−
h(x(a))(∇ f (x(a))) T x (a)
x (a)
h(x(b))(∇ f (x(b))) T x (b)
x (b)
= h x(a) t x(a) , ∇c f x(a) ,
= h x(b) t x(b) , ∇c f x(b) .
Hence (2.11) is proved. Exchanging h and f in (2.11), we obtain
C
f c h + ∇c f , ∇c h ds =
2
f (xi ) t(xi ), ∇c h(xi ) .
i =1
Subtracting (2.13) from (2.11), we obtain (2.12).
2
(2.13)
68
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
2.3. Global properties of space curves Theorem 2.4. Let C be a closed smooth space curve and L be its arc-length. Then we have
x T κ ds.
L=−
(2.14)
C
Furthermore, if C is a closed smooth plane curve, we have
A=
1
x T N ds,
2
(2.15)
C
where A stands for the area of the region enclosed by the curve, N(x) is the outward normal of the curve at x. Proof. From (2.7) and the Green’s formula, we have
−
x
T
C
T
x T c x(u ), y (u ), z(u )
κ ds =
ds
C
T
∇c x, ∇c x + ∇c y , ∇c y + ∇c z, ∇c z
=
ds
C
+ = C
2
T
x(xi ) t(xi ), ∇c x(xi ) + y (xi ) t(xi ), ∇c y (xi ) + z(xi ) t(xi ), ∇c z(xi )
i =1
(∇ x)T x (x )T ∇ x + (∇ y )T x (x )T ∇ y + (∇ z)T x (x )T ∇ z x(b) T x (b) x(a) T x (a) ds + − 2 x (b) x (a) x
=
ds. C
Hence (2.14) is derived. The following facts are used in the derivation above:
x(a) = x(b),
x (a) = x (b),
t x(a) = −t x(b) .
Note that if the space curve is not closed, the boundary term does not vanish. Hence (2.14) is not true. Equality (2.15) is a well known fact (see Dierckx, 1995, p. 38), we ignore its proof. 2 Remark 2.1. The arc-length formula (2.14) depends on
κ , but is independent of the torsion of the curve.
Corollary 2.1. Let C be a closed smooth space curve. Then for any vector a ∈ R3 , we have
a T κ ds = 0.
(2.16)
C
Furthermore, if C is a closed smooth plane curve, then for any vector a ∈ R2 , we have
a T N ds = 0.
(2.17)
C
Proof. In (2.14), taking x as x + a, and noticing the fact that the translation does not change the arc-length of the curve, we obtain (2.16). (2.17) can be similarly derived from (2.15). 2 Corollary 2.2. Let C be a closed smooth space curve. We have
κ ds = 0.
(2.18)
C
Furthermore, if C is a closed smooth plane curve, we have
N ds = 0. C
(2.19)
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
Proof. In (2.16), taking a = [1, 0, 0] T , [0, 1, 0] T , [0, 0, 1] T , respectively, we obtain (2.18). (2.19) is similarly proved.
69
2
Theorem 2.5. Given a smooth curve motion flow
∂x = V (x), ∂t then for a closed space curve, we have
dL (t )
V (x) T κ (x) ds,
=−
dt
(2.20)
C(t )
and for a closed plane curve, we have
d A (t ) dt
1
=
V (x) T N(x) ds,
2
(2.21)
C(t )
where N(x) is the outward normal of the curve at x. Proof. For plane curves, the theorem is well known (see Gage, 1986, and Dierckx, 1995, p. 41). Now we prove (2.20) for space curves. From
b
L (t ) =
ds =
x du ,
a
C(t )
and Green’s formula, we have
dL (t ) dt
b = a
∂x du = ∂t x (b)
=
T
x (b)
x (b)
=
b a
T
x
∂ x du ∂t
x
T b T ∂ x(b) x (a) ∂ x(a) x ∂x − − du ∂t x (a) ∂t x ∂t a
T
x (b)
V x(b) −
x (a)
T
x (a)
V x(a) −
b a
x
x
T V (x) du .
(2.22)
If the curve considered is smooth and close, then we have
x (b)
T
V x(b) =
x (b)
x (a)
x (a)
T
V x(a) .
(2.23)
From (2.3), we have
x
x
=
x x 2 − x x T x
x 3
= κ (x) x .
Then from (2.22), we obtain (2.20). The theorem is proved.
2
If the curve considered is not closed, the equality (2.23) is not valid in general. Hence, (2.20) may not be true. However, if V (x(a)) and V (x(b)) are in the normal direction of the curve, then the terms related to the boundary vanish. Then, we obtain the following corollary. Corollary 2.3. Given a space curve normal-motion flow
∂x = v (x)n (x), ∂t we have
dL (t ) dt
=− C(t )
v (x)κ 2 (x)
κ (x) +
ds.
(2.24)
70
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
Note that N(x) and n(x) are different. N(x) is the outward normal of the plane curve, n(x) point to the center of the osculating circle of the curve at x. If the curve is convex, then N(x) = −n(x). 2.4. Area changing rate for space curves The concept of area enclosed by a close curve is valid only for a close plane curve. For a space curve, the concept of area is not well defined. In the following, we replace the concept of the enclosed area by the concept of the area of the swept surface of the geometric flow. Given a sufficiently smooth space curve C = {x(u ) ∈ R3 : u ∈ [a, b]}. Without loss of generality, we assume that the curve does not contain straight line segment. If the curve contains straight line segments, we divide the curve into several pieces at the end points of the straight line segments, such that each piece is either a straight line or contains no straight line. If a piece of the curve is a straight line, then we say this curve piece is trivial. Otherwise, we say the curve piece is non-trivial. We consider only the non-trivial curves in the following. Let u 1 < u 2 < · · · < un be all the zero points of κ (u ) of the curve C. Then around ul , there exists a vector nl ∈ R3 , such that κ (u ) and κ (u ) are in the following forms
κ (u ) = (u − ul )αl nl + O (u − ul )αl +1 ,
κ (u ) = |u − ul |αl nl + O |u − ul |αl +1 , where
(2.25) (2.26)
αl 1 is an integer. Then around ul ,
n x( u ) =
α nl κ (u ) = sign(u − ul ) l + O |u − u l | . κ (u ) nl
We define a sign for n(x) in the interval (ul , ul+1 ) as follows
sl =
sl−1 , if αl an even number, −sl−1 , if αl an odd number,
where s0 is defined as 1. Definition 2.1. If C is an open curve or a close curve (x(a) = x(b)) with sn = 1, then we say the curve is orientable. Otherwise, the curve is not orientable. If curve C is orientable, then a sign function as follows is well defined:
sign x(u ) = sl ,
if u ∈ (ul , ul+1 ),
and sign(x(ul )) = 0 for l = 1, . . . , n. Having a sign function sign(x(u )), we define an oriented normal function as
N x( u ) =
sign(x(u ))n(x(u )), n
sl−1 nl , l
κ (u ) = 0, u = u l , l = 1, . . . , n .
Then N(x) is a smooth vector-valued function and
κ (x) = sign(x)κ (x)N(x) = κs N(x), where κs = sign(x)κ (x) is signed curvature. It is a smooth function if the curve x(u ) is sufficiently smooth. Given a curve normal-motion flow for a non-trivial curve
∂x = V (x). ∂t The area of the surface swept by the curve motion is
t V (x) T N(x) ds dτ .
A (t ) = 0 C(τ )
Then the changing rate of A (t ) is
d A (t ) dt
= C(t )
V (x) T N(x) ds.
(2.27)
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
71
Fig. 3.1. Left column: The initial space curve to be evolved. Middle column: The evolved space curve family at the early stage using curvature flow (3.1). Right column: The evolved space curve family at the later stage, where the curves become circular and shrink. The figures in the first row are the space curves observed in the z-direction [0, 0, 1] T . The figures in the second row are the space curves observed in the direction [0.1, 0.9, 0.1] T .
3. Second order curve flows In this section, we propose several second-order geometric flows. Some properties of these flows are discussed. These properties are illustrated by some graphical examples. For ease of illustration, all curves used, except for Fig. 3.1, are planar. The figures presented are not schematic. They have been produced by real computation using the numerical method presented in Xu (2013). 3.1. Curvature flow Let x(u ) be a given space curve with u ∈ [a, b]. Now we define a second order flow. As for the mean-curvature flow for surface, the curvature flow for curve is defined as
∂x = κ n. ∂t
(3.1)
Remark 3.1. The counterpart in the following of the mean-curvature flow ∂∂xt = HN
∂x = κN ∂t is not correct, since the right-hand side depends on the orientation of N. Flow (3.1) can be derived as the gradient flow for minimizing the arc-length functional. For a space curve, from (2.20), we have
dL (t ) dt
κ (x)T κ (x) ds = −
=− C(t )
κ (x) 2 ds 0.
C(t )
Hence, curvature flow is arc-length shortened. For the rate of change of area enclosed by a plane curve with rotation index 1, we have
d A (t ) dt
=
1
T
κ (x)n(x) N(x) ds = −
2 C(t )
where
sign(x)κ (x) ds = −
2 C(t )
sign(x) = −n(x) T N(x) =
1
1
κs (x) ds = −π ,
2
(3.2)
C(t )
1, C is convex at x, −1, C is concave at x.
Here, the convex and concave curves are distinguished when the observer is located in the interior of the closed curve. From (3.2), we know that the curvature flow is shrinking in area. Fig. 3.1 shows an example of the curvature flow, where the left column shows the initial space curve to be evolved. The middle column shows the evolved space curve family at the early stage. The right column shows the evolved space
72
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
Fig. 3.2. Left: The initial curve to be evolved. Middle: The evolved curve family at the early stage using averaged curvature flow (3.4). Right: The steady solution of the averaged curvature flow with the initial curve.
curve family at the later stage, where the curves become circular and shrink. The figures in the first row are the space curves observed in the z-direction [0, 0, 1] T . The figures in the second row are the space curves observed in the direction [0.1, 0.9, 0.1] T . Note that we do not choose [0, 1, 0] T as the view direction, since the curves are almost in the xy-plane at the later stage. 3.2. Averaged curvature flow Let
κ0 = C
κ2 κ +
ds
C
κ ds, κ +
(3.3)
where C = {x(u ) ∈ R3 : u ∈ [a, b]}. Then the counterpart of the averaged mean curvature flow is defined as
∂x = (κ − κ0 )n . ∂t
(3.4)
Since
C(t )
κ (κ02 − κ0 κ ) ds = κ02 κ +
κ
C(t )
ds − κ0
κ +
C(t )
κ2 κ +
ds = 0,
we have
dL (t ) dt
=− C(t )
κ (κ 2 − κκ0 ) ds = − κ +
C(t )
κ (κ 2 − κκ0 ) κ (κ02 − κ0 κ ) + ds = − κ + κ +
C(t )
κ (κ − κ0 )2 ds 0. κ +
Therefore, the averaged curvature flow is arc-length shortened for a space curve. For plane curves, from (2.21), we have
d A (t ) dt
=
1
2 C(t )
κ (x)(κ (x) − κ0 (t ))n(x)T N(x) 1 ds = − κ + 2
sign(x) C(t )
κ 2 (x) − κ (x)κ0 ds κ +
= 0 if sign(x) = 1. If the curve is convex, the flow is area preserving. For a space curve, if κ > 0 everywhere, then from (2.27), we know the flow is area preserving. Fig. 3.2 shows an example of the averaged curvature flow, where an initial non-convex curve is evolved. Different from the curvature flow, the evolved curves are not shrinking. It is natural to believe that the steady solution of the flow is a circle with radius 1/κ0 . However, this is not true: the last figure shows that the steady solution of the flow is a curve consisting of piecewise circular arcs with radius 1/κ0 . The computational result shows that the total squared error C(t ) (κ − κ0 )2 ds approaches to zero. Hence, flow (3.4) is different from flow (1.1), even for the case of a plane curve. However, if = 0 and the curve considered is a close and convex plane curve with winding number 1, then from (3.3), we have
κ0 =
κ ds C
ds = C
1
κs ds =
L
2π L
,
C
and therefore, flow (3.4) is same as the flow (1.1).
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
73
Fig. 3.3. Left: The initial curve to be evolved. Middle: The evolved curve family at the early stage using averaged curvature flow (3.5). Right: The evolved curve family at the later stage.
3.3. Area-preserved averaged curvature flow Flow (3.4) is area-preserved for the curve with κ > 0. If the curvature has zero points, the area preserving property is not valid. Now we define an area-preserved averaged curvature flow. Let
1
κ0 =
κs ds.
L C
Then the flow
∂x = (κs − κ0 )N ∂t
(3.5)
is arc-length shortening and area preserving. In fact, from
κ02 ds − C(t )
κ0 κs ds = 0,
C(t )
and (2.20), we have
dL (t ) dt
=−
κs 2 − κs κ0 ds = −
C(t )
κs 2 − κs κ0 + κ02 − κ0 κs ds = −
C(t )
(κs − κ0 )2 ds 0.
C(t )
Therefore, the averaged curvature flow is arc-length shortening for a space curve. Furthermore, from (2.27), we have
d A (t )
=
dt
1
2
κs (x) − κ0 (t ) N(x)T N(x) ds = 0.
C(t )
Hence, the flow is area preserving. Fig. 3.3 shows an example of the area-preserved averaged curvature flow, where the left figure shows an initial nonconvex curve to be evolved. The middle figure shows the evolved curve family at the early stage. The right figure shows the evolved curve family at the later stage. The curve approaches a circle with radius 1/κ0 . 3.4. Arc-length preserved averaged curvature flow Let
κ1 = C
κ3 κ +
ds
C
κ2 κ +
ds.
Then we introduce another averaged curvature flow as follows
∂x = (κ − κ1 )n . ∂t Since
dL (t ) dt
=− C(t )
κ 3 − κ 2 κ1 ds = 0. κ +
(3.6)
74
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
Fig. 3.4. Left: The initial curve to be evolved. Middle: The evolved curve family at the early stage of the evolution using the averaged curvature flow (3.6). Right: The evolved curve family at the later stage.
Therefore, this averaged curvature flow is arc-length preserved for evolving space curves. If C(t ) is a convex plane curve, flow (3.6) is similar to (1.2) proposed by Ma and Zhu (2012). If C(t ) is not convex, then (3.6) is not similar to (1.2). Fig. 3.4 shows an example of this averaged curvature flow, where an initial non-convex curve is evolved. Different from the curvature flow, the evolved curves are not shrinking. Again, the steady solution of the flow is a curve consisting of piecewise circular arcs with radius 1/κ1 . The computational result shows that the total squared error C(t ) (κ − κ1 )2 ds approaches to zero. 4. Fourth order curve flows In this section, several fourth-order flows, including curve diffusion flow, quasi-curve diffusion flow, projective quasicurve diffusion flow, minimal curvature flow and minimal squared-curvature flow, are presented. Weak forms of some these flows are also given. The illustrations provided are generated using the numerical method presented in Xu (2013). 4.1. Curve diffusion flow To define a curve diffusion flow which is similar to surface diffusion flow, let us recall the surface diffusion flow (see Mullins, 1957):
∂x = −2s Hn, ∂t
(4.1)
where s , H and n denote the Laplace–Beltrami operator, mean curvature and surface normal, respectively. Note that mean curvature H depends on the orientation of the surface normal. Namely, if we take the opposite sign of the normal, then the mean curvature will change sign. Hence, the right-hand side of (4.1) is independent of the orientation of the surface normal. If we take directly the form of the surface diffusion flows as the curve diffusion flow, we obtain the following form for curve diffusion flow:
∂x = −c κ n, ∂t
(4.2)
where c and n are the Laplace operator and principal normal of the curve. However, no matter how smooth the curve x(u ) is, κ (u ) can be C 0 smooth (see (2.26)) around zero points of κ , and c κ is not well defined. Therefore, curve diffusion flow (4.2) is not well defined at the points where κ = 0. However, the flow
∂x = −c κs N ∂t
(4.3)
is well defined. From (2.20), (2.27) and Green’s formulas, we have
dL (t ) dt
c κs NT κ ds =
= C(t )
and
d A (t ) dt
c κs κs ds = −
C(t )
=− C(t )
∇c κ , ∇c κ ds 0,
C(t )
c κ NT N ds = − C(t )
c κ ds =
∇c κ , ∇c 1 ds = 0.
C(t )
Hence, the flow is arc-length shortening and area-preserving. Fig. 4.1 shows an example of the curve diffusion flow, where an initial un-convex curve is evolved. The middle figure shows the evolved curve family at the early stage using the curve diffusion flow (4.3). The right figure shows the evolved curve family at the later stage.
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
75
Fig. 4.1. Left: The initial curve to be evolved. Middle: The evolved curve family at the early stage using curve diffusion flow (4.3). Right: The evolved curve family at the later stage using the curve diffusion flow (4.3).
Fig. 4.2. Left: The initial curve to be evolved using quasi-curve diffusion flow (4.4). Middle: The evolved curve family at the early stage of the evolution. Right: The evolved curve family at the later stage of the evolution.
4.2. Quasi-curve diffusion flow For surface, we had introduced the following quasi-surface flow (see Xu, 2008)
∂x = −2s x. ∂t The flow has area shrinking property. Similar to this, we introduce the following quasi-curve diffusion flow
∂x = −c κ . ∂t
(4.4)
From (2.3), we can see that, for a sufficiently smooth curve, c κ is well defined though c κ may not be well defined. From (2.20) and Green’s formula, we have
dL (t ) dt
=
c κκ ds = −
C(t )
∇c κ : ∇c κ ds 0.
C(t )
Hence, flow (4.4) is arc-length shortening for a close space curve. From Frenet–Serret formulas, we can derive that, if
κ > 0,
c κ = κ − τ 2 − κ 2 n − 3κκ t + κτ b = c κ − τ 2 − κ 2 n − 3κ ∇c κ + κτ b. Hence, the quasi-curve diffusion flow contains the motions in the tangent and binormal directions. Note that torsion is involved in this equation. Using Green’s formula (2.12), the weak form of the quasi-curve diffusion flow is written as
C(t )
∂x φ ds = − ∂t
c xc φ ds,
∀φ ∈ C 02 (Ω),
(4.5)
C(t )
where C 02 (Ω) consists of C 2 smooth functions over Ω with support in the interior of Ω . Fig. 4.2 shows an example of the quasi-curve diffusion flow, where an initial non-convex curve is evolved. The middle figure shows the evolved curve family at the early stage using the quasi-curve diffusion flow. The right-hand figure shows the evolved curve family at the later stage. The curve will shrink with further evolution.
76
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
Fig. 4.3. Left: The initial curve to be evolved using the projective curve diffusion flow (4.6). Middle: The evolved curve family at the early stage of the evolution. Right: The evolved curve family at the later stage.
4.3. Projective quasi-diffusion flow As we have mentioned before, the quasi-curve diffusion flow contains the motions in the tangent and binormal directions. To disable the tangential motion, we construct the following form flow:
∂x = − I − ttT c κ , ∂t
(4.6)
where t is the unit tangent vector of the curve. Since κ does not depend on the orientation of the normal, the right-hand side of (4.6) is independent of the orientation of the normal. The right-hand side of (4.6) is the projection of c2 x in the normal plane, we therefore name it as projective quasi-diffusion flow. From (2.20) and Green’s formula, we have
dL (t ) dt
(c κ )T I − ttT κ ds =
= C(t )
(c κ )T κ ds = −
C(t )
∇c κ : ∇c κ ds 0,
C(t )
where A : B denotes trace of A B . Hence, flow (4.6) is arc-length shortening for evolving close space curves. Let y = κ . Then, using the Green’s formula, we can write the weak-form of (4.6) as: Find (x, y) such that T
⎧ ∂x ⎪ ⎪ φ ds = − I − tt T c yφ ds, φ ∈ C 02 (Ω), ⎪ ⎪ ∂ t ⎪ ⎨ C(t ) C(t ) ⎪ ⎪ ⎪ yψ ds = − ∇c x∇c ψ ds, ψ ∈ C 02 (Ω). ⎪ ⎪ ⎩ C(t )
(4.7)
C(t )
Fig. 4.3 shows an example of the projective curve diffusion flow (4.6), where an initial non-convex curve is evolved. The middle figure of Fig. 4.3 shows a curve family at the early stage obtained by this flow. The right-hand figure shows a curve family at the later stage. The curve will shrink continuously for further evolution. 4.4. Minimal squared-curvature flow For surfaces, Willmore flow (see Willmore, 1993) is derived from minimizing the total squared mean-curvature. For curves in R3 , we minimize the following energy
κ 2 ds.
E (C) = C
To compute its first order variation, we represent E as
E (C) =
b
2
2
κ ds = C
c x ds =
c x2 x du .
(4.8)
a
C
Assume C is an extreme point of the energy functional (4.8). Let
C(Φ, ε ) = x(u , ε ) = x(u ) + ε Φ(u ): u ∈ Ω ,
|ε | 1,
be a perturbation of the curve C, where Φ is smooth vector-valued function with compact support. Then
δ E (C), Φ =
d dε
E C(Φ, ε )
ε =0
= 0,
(4.9)
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
77
where
b
δ E (S), Φ =
2
δ c x x
b du =
a
2 c x, δ(c x) x + κ 2 δ x du .
(4.10)
a
It is not difficult to derive that
(x )T Φ δ x = , x Φ Φ (x )T x + x (Φ )T x + x (x )T Φ 2x (x )T x (x )T Φ 2x (x ) T Φ − + δ(c x) = 2 − 4 x x x 4 x 6 2x (x ) T Φ x (Φ ) T x + x (x ) T Φ 2x (x ) T x (x ) T Φ = c Φ − − + . 4 4 x x x 6 Notice that c x, x = 0. From (4.10) we have
δ E (S), Φ =
2c x, c Φ −
2x (x ) T Φ
x 4
C
=
−
2c x, c Φ ds +
b
=
−
2c x, c Φ ds + a
C
ds +
κ 2 (x )T Φ x
4x (x ) T c x
x 3
a
C
b a
b
3κ 2 x
x
+
κ 2 x x
du
, Φ du
, Φ du .
(4.11)
Using Green’s formulas, we have
δ E (S), Φ = 2
2 c x, Φ ds +
C
1
3κ 2 x
x
C
x
, Φ ds.
Since Φ is an arbitrary vector-valued function, we obtain, from (4.9), the following Euler Lagrange equation
1
2c2 x +
x
3κ 2 x
x
= 0.
To simplify the derivation, but without loss of generality, we assume the parameter of the curve is the arc-length. Then, we have
1
x
3κ 2 x
x
= 3 κ 2 x = 3κ 2 κ + 3 κ 2 x .
Since term 3(κ 2 ) x is in the tangential direction, we then construct the following geometric flow by ignoring the tangential term:
∂x = −2c κ − 3κ 2 κ . ∂t
(4.12)
We name this geometric flow as minimal squared-curvature flow. From (2.20) and Green’s formula, we have
dL (t ) dt
=
2c κκ + 3κ 4 ds =
C(t )
3κ 4 − 2∇c κ : ∇c κ ds = 3
C(t )
C(t )
κ 4 ds − 2
2 x ds,
C(t )
where the third-order derivative x is taken with respect to the arc-length parameter. If a space curve is bumpy, then the term 2 C(t ) x 2 ds is dominant. The flow is arc-length shortening. If a space curve is not bumpy, the term 3 C(t ) κ 4 may be dominant. The flow is arc-length lengthening. For instance, if the curve C is a circle, then we can derive that
C(t )
3κ 4 − 2∇c κ : ∇c κ ds =
κ 4 ds > 0.
C(t )
Hence, for a circle, flow (4.12) is arc-length lengthening.
78
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
Fig. 4.4. Left: The evolved curve family at the early stage using minimal squared-curvature flow (4.12) for the initial curve as shown in Fig. 4.2. Middle: The evolution at a stage where the curve becomes a circle. Right: The evolved curve family at the later stage. The circles become larger and larger.
Using (4.11) and noting that ∇c φ = ∇c x∇c φ , we obtain the weak-form of the minimal squared-curvature flow as
∂x φ ds = − ∂t
C(t )
2c xc φ − 3κ 2 ∇c x∇c φ ds,
∀φ ∈ C 02 (Ω).
(4.13)
C(t )
Fig. 4.4 shows an example of the minimal squared-curvature flow, where an initial non-convex curve, as shown in Fig. 4.2, is evolved. The left-hand figure shows an evolved curve family at the early stage of the evolution. It is easy to see that the curve is shortening. The middle figure shows a stage where the curve becomes a circle. The right-hand figure shows an evolved curve family at the later stage. The curve becomes a circle, and further evolution makes the circle larger and larger. 4.5. Minimal curvature flow As we have seen, the minimal squared-curvature flow is arc-length lengthening for a circle. To obtain a flow that is not arc-length lengthening for a circle, we minimize the following energy:
E (C) =
κ ds.
(4.14)
C
Assume C is an extreme point of the energy functional (4.8). Let
C(Φ, ε ) = x(u , ε ) = x(u ) + ε Φ(u ): u ∈ Ω ,
|ε | 1,
be a perturbation of the curve C, where Φ is smooth vector-valued function with compact support. Then
d δ E (C), Φ = E C(Φ, ε ) dε
ε =0
= 0,
(4.15)
where
b
δ E (S), Φ =
δ κ x du .
(4.16)
a
Without loss of generality, we may assume u is arc-length. Then it is not difficult to derive that
δ κ x = δ
x 2 x 2 − (x T x )2 x 2
= n, Φ − x x , Φ .
Assume n = 0. Then using Green’s formulas, we have
δ E (S), Φ =
C
c n + κκ + x κ , Φ ds.
(4.17)
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
79
Fig. 4.5. The left-hand figure shows the initial curve to be evolved. The remaining figures show the evolved curves using minimal curvature flow (4.18) at different stages of the evolution.
Since Φ is an arbitrary vector-valued function, then from (4.15) we obtain the following Euler Lagrange equation
c n + κκ + x κ = 0. Since term x κ is in the tangential direction, we then construct the following geometric flow by ignoring the tangential term
∂x = −[c n + κκ ]. ∂t
(4.18)
We name this flow as minimal curvature flow. From (2.21) and Green’s formula, we have
d A (t ) dt
=
1
2
κ 2 − ∇c n : ∇c n ds.
C(t )
Suppose the curve C is a plane convex curve, then from the Frenet–Serret formula we can derive that
κ 2 − ∇c n : ∇c n ds = 0.
C(t )
Hence, for a plane convex curve, flow (4.18) is area preserving. In fact, if C(t ) is convex plane curve, then C(t ) κ ds = 2π is a constant. This means that C(t ) is a minimizer of the energy functional (4.14). From the Frenet–Serret formula we can derive that, if κ > 0, (4.18) becomes
∂x = −τ b + tκ + τ 2 n, ∂t where the derivative is taken with respect to the arc-length s. From (4.17), the weak-form of the flow (4.18) can be written as
C(t )
∂x φ ds = − ∂t
c xc φ
κ
− κ ∇c x∇c φ ds,
∀φ ∈ C 02 (Ω).
(4.19)
C(t )
The right-hand side of (4.19) has a singularity when κ = 0, which reflects the fact that the normal n is not continuous. To 1 overcome this difficulty, we replace the term κ1 with κ + , where > 0 is a small number. Fig. 4.5 shows an example of the minimal curvature flow, where an initial non-convex curve is evolved. Note that the resulting curve has a sharp feature. This is understandable since the energy (4.14) is, in fact, the L 1 norm of the curvature. Remark 4.1. Note that (4.18) is derived under the assumption that n = 0. If the curve C has points where n = 0 (x and x are not linearly independent), this equation cannot be used. However, the weak form (4.19) is valid. 5. Conclusion and future work Several second-order and fourth-order geometric flows have been constructed for space curves. Some properties of the changing rates of the arc-length and areas have been discussed. These properties show that the geometric flows constructed for space curves are different in nature from their counterparts of the plane curve case and surface case. Short-term and long-term behaviors of the evolved curves using the constructed flows have been illustrated by graphical examples. These behaviors show that each of the flows has its own distinct features. As part of future work, we plan to study the applications of the geometric flows constructed. A direct application may be space curve denoising. All the flows with arc-length shortening property can be used for space curve denoising. However,
80
G. Xu, X. Yang / Computer Aided Geometric Design 31 (2014) 63–80
if a flow is both arc-length shortening and area-preserving, then it would be more desirable, since it will better preserve the shape of the curve. Another application can be image restoration and image denoising. For instance, the 2D version of flow (4.5) can be used to accomplish the task with sharp feature of the objects in the image to be preserved. A third application area of the geometric flows can be geometric design. In many shape design problems, shapes are designed from a wireframe consisting of space curves. Each of the smooth space curves, say B-spline curves, is constructed from interpolating a few fixed control points. The remaining degrees of freedom (DOF) of the curve need to be determined. These DOF can be determined by the geometric flows constructed. The space curves considered in this paper are expressed in parametric forms. If the curves are in the level-set form (intersections of two implicit form surfaces), how to reformulate these flows is an interesting problem that is worth further study. References Altschuler, S.J., Grayson, M.A., 1992. Shortening space curves and flow through singularities. J. Differ. Geom. 35 (2), 283–298. Andrews, B., 1998. Evolving convex curves. Calc. Var. 7, 315–371. Chen, S., Lin, M., Wu, J., 2007. A weighted curvature flow for planar curves. Int. J. Math. Comput. Simul. 1 (3), 275–278. Chou, K., 2001. The curve shortening flow: The classical approach. In: RIMS Kokyuroku, vol. 1210. Cohen, E., Tekumalla, L., 2005. Smoothing space curves with the MLS projection. In: Geometric Modeling, Vizualization and Graphics. Salt Lake City. Dierckx, P., 1995. Curve and Surface Fitting with Splines. Oxford University Press. Do Carmo, M.P., 1976. Differential Geometry of Curves and Surfaces. Prentice Hall, Englewood Cliffs, NJ. Dolcetta, I.C., Vita, S.F., March, R., 2002. Area preserving curve shortening flows: From phase transitions to image processing. Interfaces Free Bound. 4 (4), 325–343. Dziuk, G., 1999. Discrete anisotropic curve shortening flow. SIAM J. Numer. Anal. 36 (6), 1808–1830. Gage, M., 1983. An isoperimetric inequality with application to curve shortening. Duke Math. J. 50, 1225–1229. Gage, M., 1984. Curve shortening makes convex curves circular. Invent. Math. 76, 357–364. Gage, M., 1986. On an area-preserving evolution equation for plane curves. In: Contemp. Math., vol. 51, pp. 51–62. Gage, M., 1990. Curve shortening on surfaces. Ann. Sci. Éc. Norm. Super. 23, 229–256. Gage, M., Hamilton, R., 1986. The heat equation shrinking convex plane curves. J. Differ. Geom. 23, 69–96. Grayson, M., 1987. The heat equation shrinks embedded plane curves to points. J. Differ. Geom. 26 (2), 285–314. Grayson, M., 1989. Shortening embedded curves. Ann. Math. 129, 71–111. Guziec, A., Ayache, N., 1994. Smoothing and matching of 3-d space curves. Int. J. Comput. Vis. 12 (1), 79–104. Hildebrandt, K., Polthier, K., Preuss, E., 2005. Evolution of 3d curves under strict spatial constraints. In: Proceedings of the Ninth International Conference on Computer Aided Design and Computer Graphics. CAD-CG’05. IEEE Computer Society, Washington, DC, USA, pp. 40–45. Leng, J., Zhang, Y., Xu, G., 2011. A novel geometric flow-driven approach for quality improvement of segmented tetrahedral meshes. In: Proceedings of the 20th International Meshing Roundtable, pp. 347–364. Liu, Y., Lai, Y., Zang, Y., Hu, S., 2008. Fairing wireframes in industrial surface design. In: 2008 IEEE International Conference on Shape Modeling and Applications. Stony Brook, New York. ACM, pp. 29–35. Ma, L., Chen, D., 2007. Curve shortening in a riemannian manifold. Ann. Mat. Pura Appl. 186, 663–684. Ma, L., Zhu, A., 2012. On a length preserving curve flow. Monatshefte Math. 165 (1), 57–78. Mullins, W.W., 1957. Theory of thermal grooving. J. Appl. Phys. 28, 333–339. Pan, S., Yang, J., 2008. On a non-local perimeter-preserving curve evolution problem for convex plane curves. Manuscr. Math. 127, 469–484. Sapiro, G., 2001. Geometric Partial Differential Equations and Image Analysis. Cambridge University Press, Cambridge. Sapiro, G., Tannenbaum, A., 1995. Area and length preserving geometric invariant scale-spaces. IEEE Trans. Pattern Anal. Mach. Intell. 17, 67–72. Spira, A., Kimmel, R., 2007. Geometric curve flows on parametric manifolds. Chin. J. Comput. Phys. 223 (1), 235–249. Willmore, T.J., 1993. Riemannian Geometry. Clarendon Press, Oxford, England. Xu, G., 2008. Geometric Partial Differential Equation Methods in Computational Geometry. Science Press, Beijing, China. Xu, G., 2013. Numerical integration of geometric flows for space curves. Research Report No. ICM-13-01. Institute for Computational Mathematics and Scientific/Engineering Computing.
