Evolution of plane curves with a curvature adjusted tangential velocity
aa r X i v : . [ m a t h . NA ] J a n submitted to : Japan J. Indust. Appl. Math. running title : Curvature adjusted tangential velocity
Evolution of plane curves with a curvature adjusted tangential velocity ∗ Daniel ˇSEV ˇCOVI ˇC † and Shigetoshi YAZAKI ‡ Abstract.
We study evolution of a closed embedded plane curve with the normalvelocity depending on the curvature, the orientation and the position of the curve.We propose a new method of tangential redistribution of points by curvature adjustedcontrol in the tangential motion of evolving curves. The tangential velocity distributesgrid points along the curve not only uniform but also lead to a suitable concentrationand/or dispersion depending on the curvature. Our study is based on solutions tothe governing system of nonlinear parabolic equations for the position vector, tangentangle and curvature of a curve. We furthermore present a semi-implicit numericaldiscretization scheme based on the flowing finite volume method. Several numericalexamples illustrating capability of the new tangential redistribution method are alsopresented in this paper.
Key Words: curvature driven flow of a plane curve, curvature adjusted tangentialvelocity, semi-implicit scheme, flowing finite volume method, crystalline curvature flowequation.
Mathematics Subject Classification (2000):
The purpose of this paper is to study evolution of a plane curve Γ t , t ∈ [0 , T ), driven bythe normal velocity v which is assumed to be a function of the curvature k , tangentialangle ν and position vector x ∈ Γ t , v = β ( x , ν, k ) . (1.1)Geometric equations of the form (1.1) often arise from various applied problems like e.g.the material science, dynamics of phase boundaries in thermomechanics, computationalgeometry, semiconductors industry, image processing and computer vision, etc. For a ∗ The authors were/are supported by grants: VEGA 1/0657/11 (DˇS) and Grant-in-Aid for Encour-agement of Young Scientists 21740079 (SY). † Institute of Applied Mathematics, Faculty of Mathematics, Physics and Informatics, Comenius Uni-versity, 842 48 Bratislava, Slovak Republic.
E-mail : [email protected] ‡ Faculty of Engineering, University of Miyazaki, 1-1 Gakuen Kibanadai Nishi, Miyazaki 889-2192,Japan.
E-mail : [email protected] x : R / Z ⊃ [0 , → R such that Γ = Image( x ) = { x ( u ); u ∈ [0 , } and | ∂ u x | >
0. Hereafter, we denote ∂ ξ F = ∂ F /∂ξ , and | a | = √ a · a where a · b is theEuclidean inner product between vectors a and b . The unit tangent vector can be definedas t = ∂ u x / | ∂ u x | = ∂ s x , where s is the arc-length parameter and ds = | ∂ u x | du . The unitinward normal vector n is defined in such a way that det( t , n ) = 1. The signed curvaturein the direction n is denoted by k . With this orientation we have k = det( ∂ s x , ∂ s x ).Let ν be the angle of t , i.e., t = (cos ν, sin ν ) T and n = ( − sin ν, cos ν ) T . The problemof evolution of curves is stated as follows: For a given initial curve Γ = Image( x ),find a family of planar curves { Γ t } ≤ t 0) = x ( u ); u ∈ [0 , } and evolving in the normal direction according to thevelocity given by (1.1).As typical examples of β one can consider β = w ( ν ) k representing the anisotropicmean curvature flow, β = | k | γ − k (power like flow arising in image processing) and β = w ( x , ν ) k + F ( x , ν ), where w is a weight function, F is an external force and γ > v is the normal component of the following evolutionaryequation for the position vector x : ∂ t x = β n + α t , x ( · , 0) = x ( · ) . (1.2)Here α is the tangential component of the velocity vector. Notice that the motion in thetangential direction with a tangential velocity α has no effect of the shape of evolvingclosed curves [7, Proposition 2.4]. The shape is determined by the value of the normalvelocity β only. Hence the trivial setting α ≡ β = k and α = 0. However, it was documented by many authors(see e.g. [10, 14, 15, 29] and references therein) that such a choice of α may lead tovarious numerical instabilities caused by either undesirable concentration and/or extremedispersion of numerical grid points.In order to construct a stable numerical computational scheme, various choices ofa nontrivial tangential velocity have been proposed and analyzed by many authors. Wepresent a brief review of development of nontrivial tangential velocities. In [12, 13] Kimuraproposed a uniform redistribution scheme for the case when β = k by using a specialchoice of α satisfying the uniform redistribution condition (U): | ∂ u x | = L t , where L t istotal length of a curve Γ t . The author also proved convergence of a numerically discretizedsolution when assuming uniform distribution of initial grid points. In [10], Hou, Lowen-grub and Shelley (see also Hou, Klapper and Si [11]) utilized the condition (U) directly forthe case β = k . More precisely, they derived the form of a tangential velocity with ϕ ≡ ω ≡ ϕ ( k ) with ω ≡ ϕ ≡ ω ≡ β = β ( ν, k ). This result canbe considered as an improvement of that of [14] in which satisfactory results were obtainedonly for the case when β = β ( k ) is linear or sublinear function with respect to k . Next,in the series of papers [16, 17, 18], Mikula and the first author proposed a method ofasymptotically uniform redistribution. In terms of our notation, they derived (3.3) with ϕ ≡ ω ( t ) for a general class of normal velocities ofthe form β = β ( x , ν, k ). Their method was also applied to geodesic curvature flows andimage segmentation problems.Besides these uniform or asymptotically uniform redistribution methods, in the so-called crystalline curvature flow, the grid points are distributed densely (sparsely) in thosepart of the curve where the absolute value of the curvature is large (small). Although thisredistribution is far from being uniform, numerical computation is quite stable. One of thereasons for such a behavior is that polygonal curves are restricted to a class of admissiblefacet directions. In order to extract essence of the crystalline curvature flow of polygonalcurves and generalize it to a wide class of plane curve evolution, the second author showedthat the tangential velocity α = − ∂ s β/k is implicitly involved in the crystalline curvatureflow of planar curves (cf. [29]). Notice that such a tangential velocity is a local one sinceits value at some point x of a curve depends on the local properties of the curve near x .In this point, it is worthwhile to mention the paper by Deckelnick [5] who used locallydefined tangential velocity having the form α = − ∂ u ( | ∂ u x | − ) for the special case when β = k . Recently, Pauˇs et al. successfully applied curvature adjusted tangential velocityfor evolution of open curves having important applications in the dislocation dynamics(see [22, 23]).The asymptotically uniform redistribution is quite effective and applicable for a widerange of applications. However, from the approximation point of view, unless the solutioncurve is a circle, there is no reason to accept uniform redistribution for a general case ofa curve evolution. Hence the desired redistribution should take into account the shapeof evolution curves and it should also depend on the modulus of the curvature. We willfurthermore present a combination of the method of asymptotic uniform redistribution [16]and the crystalline tangential velocity discussed in [29].3s far as 3D implementation of tangential redistribution is concerned, in the recentpaper by Barrett, Garcke and N¨urnberg [1], they proposed and studied a new efficientnumerical scheme for evolution of surfaces driven by the Laplacian of the mean curvature.It turns out that a uniform redistribution of tangential velocity vectors is implicitly builtin their numerical scheme. A similar scheme with an explicit expression for the tangentialredistribution has been proposed and utilized by Morigi in [20].The organization of the paper is as follows: In Section 2, we recall the system ofPDEs for geometric equations (1.1) and (1.2). In Section 3, we will address the questionhow to control redistribution of grid points by using a nontrivial tangential velocity. Wefocus on a class of tangential velocities dependent locally on the curvature combinedwith known asymptotically uniform redistributions. The aim of Section 4 is to furthermotivate our study by constructing a curvature adjusted redistribution yielding the bestpossible approximation of the evolving family of planar curves as far as the minimizationof the error between the length of a curve and the length and enclosed of its polygonalapproximation is concerned. In Section 5, a numerical solution to the system of governingequations will be constructed by means of the so-called flowing finite volume method.In Section 6, we will present various numerical examples of applications of the curvatureadjusted tangential velocity. Finally, we conclude with some key points in Section 7. Without loss of generality, we will rewrite the normal velocity β ( x , ν, k ) as follows: β = w ( x , ν, k ) k + F ( x , ν ) . (2.1)The tangential velocity α will be defined later in Section 3. Using Fren´et formula ∂ s x = ∂ s t = k n , the following equation for the position vector follows from (1.1) and (1.2): ∂ t x = w∂ s x + α∂ s x + F n , (2.2)where the operators ∂ s and ∂ s stand for arc-length derivatives, i.e. ∂ s F = g − ∂ u F and ∂ s F = g − ∂ u ∂ s F , respectively. Here g = | ∂ u x | denotes the local length of a curve param-eterized by x . Remark 2.1 If β is a linear or superlinear function with respect to k like e.g. w ( ν ) k , w ( x , ν ) k + F ( x , ν ), | k | γ − k ( γ ≥ w appearing in (2.1) has no sin-gularity at k = 0. However, in the sublinear case, a singularity occurs at k = 0. Then aregularization at k = 0 will be necessary from both theoretical (local existence and unique-ness of solutions) as well as practical (construction of a stable numerical approximationscheme) point of view. We refer the reader to [15] for further discussion of this issue. Forexample, if w = | k | γ − , γ ∈ (0 , | k | − as follows: | e k | − = | k | − for | k | > ε and | e k | − = ε − for | k | ≤ ε , for a small ε > k = det( ∂ s x , ∂ s x ), t = ∂ s x = (cos ν, sin ν ) T , Frenet’s formulae ∂ s t = k n , ∂ s n = − k t and the relation k = ∂ s ν . The resulting system of governing PDEs togetherwith (2.2) reads as follows: ∂ t g = ( − kβ + ∂ s α ) g, (2.3) ∂ t k = ∂ s β + α∂ s k + k β, (2.4) ∂ t ν = ( ∂ k β ) ∂ s ν + ( α + ∂ ν β ) ∂ s ν + ∇ x β. t , (2.5)for u ∈ [0 , 1] and t ∈ (0 , T ), with initial conditions g ( · , 0) = g ( · ), k ( · , 0) = k ( · ), ν ( · , 0) = ν ( · ), x ( · , 0) = x ( · ), and periodic boundary conditions for g, k, x for u ∈ [0 , ⊂ R / Z and the boundary condition ν (1 , t ) = ν (0 , t ) + 2 π . The partial derivative ∇ x β is definedas ∇ x β = ( ∂ x β ( x , · , · ) , ∂ x β ( x , · , · )) T for x = ( x , x ) T . The initial functions g , k , t and x are required to satisfy the following compatibility conditions: g = | ∂ u x | > ∂ u t = g k , k = ( g ) − det( ∂ u x , ∂ u x ).In our simple and fast numerical approximation scheme proposed in Section 5 we willdiscretize (2.2) and the tangential velocity equation (3.3) with the constraint (3.4) only.With regard to [15], in the case where a more accurate numerical scheme and results arerequired, it is recommended to discretize the full system of evolutionary PDEs (2.2), (2.3),(2.4) and (2.5) for all remaining geometric quantities g , k , and ν . On the other hand, asfar as the time complexity of computation and simplicity of a numerical approximationscheme are concerned, discretization of (2.2) together with (3.3), (3.4), yields satisfactorynumerical results.However, as w may depend on k and ν , the proof of existence and uniqueness ofa smooth solution to the system of equations (2.2), (3.3), (3.4) does not seem to be astraightforward issue.Recall that in the case where α is given by the expression (3.3) with ϕ ≡ ω = 0,the short time existence and uniqueness of smooth solutions to the system of PDEs (2.2),(2.3), (2.4) and (2.5) has been shown in [15]. In the forthcoming paper [25] the authorsproved local existence and uniqueness of a classical solution to the full system of governingequations (2.2)–(2.5) and (3.3), (3.4). A key tool for construction of a suitable tangential redistribution functional α is theso-called relative local length quantity (cf. [15]). It is defined by the following ratio r : r ( u, t ) = g ( u, t ) L t , u ∈ [0 , , t ∈ [0 , T ) , L t is the total length of Γ t : L t = Z Γ t ds = Z g du, t ∈ [0 , T ) , and T > r can be explained as follows: suppose that N grid points { x ( u i , t ) } Ni =1 aredistributed on the curve Γ t for u i = i/N , i = 1 , , · · · , N . Since the arc-length s is givenby s ( u, t ) = s (0 , t ) + Z u g ( u, t ) du , we have s ( u i , t ) − s ( u i − , t ) = Z u i u i − g ( u, t ) du = L t Z u i u i − r ( u, t ) du for each i . Hence, if r ( u, t ) ≡ u at a time t , then all grid points aredistributed uniformly in the sense that s ( u i , t ) − s ( u i − , t ) = L t /N for each i at the time t .Furthermore, if the evolving family of curves fulfills the limit r ( u, t ) → u as t → T then redistribution of grid points become asymptotically uniform as t → T (cf. [16, 17]).It has been documented by many practical examples (see e.g. [16, 17]) that (asymp-totically) uniform redistribution can significantly stabilize and speed up numerical com-putation. We recall that the tangential motion has no influence on the shape of evolvedclosed planar curves. A natural question arises: Question 3.1 How to design the relative local length ratio r and, subsequently, α suchthat redistribution takes into account the shape of the limiting curve. In other words,how to densely (sparsely) redistribute grid points on those parts of a curve where themodulus of the curvature is large (small).The modulus | k | of curvature will be measured by the following shape function: Definition 3.2 Let ϕ : R → R + be a nonnegative even function ϕ ( − k ) = ϕ ( k ) > k = 0 such that ϕ ( k ) is nondecreasing for k > ϕ ( k ) ≡ ϕ ( k ) = | k | , ϕ ( k ) =1 − ε + ε | k | ( ε ∈ [0 , ϕ ( k ) = √ ε + k ( | ε | ≪ ϕ as follows: r ϕ ( u, t ) = g ( u, t ) L t ϕ ( k ( u, t )) h ϕ ( k ( · , t )) i , u ∈ [0 , , t ∈ [0 , T ) . (3.1)Here the bracket function h F i denotes the average of function F ( u, t ) along the curve Γ t : h F ( · , t ) i = 1 L t Z Γ t F ds = 1 L t Z F ( u, t ) g ( u, t ) du. r ϕ . Suppose, for a moment, that r ϕ ( u, t ) ≡ u at a time t . Then s ( u i , t ) − s ( u i − , t ) ≶ L t /N holds, if | k | satisfies ϕ ( k ) ≷ h ϕ ( k ) i on the i -th interval [ u i − , u i ]. Using this property we are in a position to construct thedesired curvature adjusted tangential redistribution which takes account of deviations ofthe shape function ϕ ( k ) from its average value h ϕ ( k ) i . Indeed, suppose that we are ableto construct a tangential velocity α in such a way that r ϕ ( u, t ) → 1, i.e., θ ( u, t ) = ln r ϕ ( u, t ) → u ∈ [0 , 1] as t → T . For t close to T , we can conclude if | k | is above/belowthe average in the sense that ϕ ( k ) ≷ h ϕ ( k ) i , then g ≶ L t holds, respectively. Distributionof grid points on corresponding subarcs is dense/sparse, respectively. Example 3.3 Consider ϕ ( k ) ≡ 1. Then θ = ln r . If one constructs α in such a way that θ ( u, t ) ≡ θ ( u, 0) for all u and t then relative local length r is preserved. On the other hand,if one takes α such that θ ( u, t ) → r → 1) for all u as t → T then redistributionbecomes asymptotically uniform. See Example 3.5.Convergence lim t → T r ϕ ( u, t ) → ω ∈ L [0 , T ) such that Z T ω ( τ ) dτ = ∞ and θ ( u, t ) = ln(( e θ ( u, − e − R t ω ( τ ) dτ + 1) for all u ∈ [0 , 1] and t ∈ [0 , T ). The previousequation can be rewritten as an ODE: ∂ t θ ( u, t ) + ω ( t )(1 − e − θ ) = 0 . (3.2) Remark 3.4 For ˜ θ = ln r ϕ one can obtain another convergence such as ˜ θ ( u, t ) = ˜ θ ( u, e − R t ω ( τ ) dτ ,then we have the equivalent ODE ∂ t ˜ θ ( u, t )+ ω ( t )˜ θ ( u, t ) = 0. The convergence rate of ˜ θ ( u, t )is the same as that of θ ( u, t ).With regard to [16, 17] one can choose the relaxation function ω ( t ) as follows: ω ( t ) = κ − κ ∂ t ln L t = κ + κ h kβ i , where κ ≥ κ ≥ ∂ t L t + R Γ t kβ ds = 0 (cf. [15]). If we formally set κ = κ = 0, thenthe function θ is constant in time t . On the other hand, if the maximal existence time isinfinite T = ∞ , we can choose κ > κ = 0. If T < ∞ and L t → t → T , wecan choose κ > 0. In both cases we have R T ω ( τ ) dτ = ∞ and so r ϕ ( u, t ) → t → T .7he equation for the tangential velocity α can be derived as follows. Differentiating θ = ln r ϕ = ln g + ln ϕ − ln R ϕg du with respect to t , using (2.3), (2.4) and the relation ∂ s ( ϕα ) = α ( ∂ s k ) ϕ ′ + ( ∂ s α ) ϕ , we obtain ϕ∂ t θ = ( − kβ + ∂ s α ) ϕ + ( ∂ s β + k β + α∂ s k ) ϕ ′ − ϕL t h ϕ i Z ∂ t ( ϕg ) du = ∂ s ( ϕα ) − f + ϕ h ϕ i h f i , where f = ϕ ( k ) kβ − ϕ ′ ( k ) ( ∂ s β + k β ), ϕ ′ ( k ) = ∂ k ϕ ( k ). Hence (3.2) holds if and only ifthe tangential velocity α satisfies the following equation: ∂ s ( ϕ ( k ) α ) ϕ ( k ) = fϕ ( k ) − h f ih ϕ ( k ) i + ω ( t ) (cid:0) r − ϕ − (cid:1) . (3.3)Equation (3.3) is written in the form ∂ s ( ϕ ( k ) α ) = F , where F := f − h f ih ϕ ( k ) i ϕ ( k ) + ω ( t ) ϕ ( k )( r − ϕ − . Since R Γ t ϕ ( k )( r − ϕ − ds = R ( L t h ϕ ( k ) i g − − ϕ ( k )) g du = 0, we obtain hF i = 0. Thusthe equation ∂ s ( ϕ ( k ) α ) = F has at least one solution α . In order to construct a uniquesolution α , we assume the following renormalization condition for α : h ϕ ( k ) α i = 0 . (3.4) Example 3.5 (uniform redistribution) In the case where ϕ ( k ) ≡ 1, the tangentialvelocity equation becomes ∂ s α = kβ − h kβ i + ω ( t )( r − − . If we set ω ≡ 0, then a solution α preserves the relative local length r (see [15]). Un-der the assumption ω κ , κ , redistribution of grid points becomesasymptotically uniform [16, 17, 18]. Example 3.6 (crystalline tangential velocity) Suppose that the evolving curve Γ t isconvex. If we consider the shape function ϕ ( k ) = | k | and ω ( t ) ≡ 0, then, with regard to(3.3) we have ∂ s ( kα ) = − ∂ s β . Taking into account renormalization constraint h ϕ ( k ) α i =0 we end up with α = − ∂ s β/k . This is exactly the same tangential velocity as it wasderived by the second author in the continuous limit of the crystalline curvature flow(see [29]). If the evolving curve Γ t has negative curvature on those parts, the shapefunction ϕ ( k ) = | k | is regularized such as, for instance, ϕ ( k ) = √ − ε + εk with ε ∈ (0 , xample 3.7 (their linear combination) In our numerical computations in Section 6,we will use the following smoothed shape function ϕ : ϕ ( k ) = 1 − ε + ε √ − ε + εk , ε ∈ (0 , . It is a linear combination of shape functions in Examples 3.5 and 3.6. Notice that ϕ ( k ) → ε → ϕ ( k ) → | k | if ε → 1, and ϕ ( k ) ≥ ϕ (0) > ε ∈ (0 , x ( l ) = ( a cos(2 πl ) , b sin(2 πl )) T be the ellipse with axes ratio a = 3 : b = 1, where l ∈ [0 , ⊂ R / Z . We will construct a new parameterization u ∈ [0 , 1] using the reparame-terization function l ( u ) such that x ( u ) = x ( l ( u )), l (0) = 0, l (1) = 1, ∂ u l ( u ) > ∂ u x ( u ) = ∂ l x ( l ( u )) ∂ u l ( u ), we have g ( u ) = g ( l ( u )) ∂ u l ( u ), where g ( l ) = | ∂ l x ( l ) | . Thenif r ϕ ≡ 1, we obtain the ODE ∂ u l ( u ) = L h ϕ ( k ) i g ( l ( u )) ϕ ( k ( l ( u )))to be solved. We have used the relation k ( u ) = k ( l ( u )). Applying the Runge-Kutta ODEsolver we calculated results depicted in Figure 3.1 (a), (b) and (c). Crystalline curvatureredistribution { x ( l i ) } is obtained from the condition t ( l i ) = ( − sin(2 πi/N ) , cos(2 πi/N )) T ,and the i -th inward normal vector of the circumscribed polygon is n ( l i ) for i = 1 , , · · · , N (Figure 3.1 (d)). The aim of this section is to further motivate the study of curvature adjusted tangentialredistribution. In what follows, we will address the question on what is the optimalredistribution of vertices { x , x , · · · , x N } belonging to a given smooth closed curve Γsuch that the discrepancy between the length/the area of Γ and that of a polygon spannedby vertices { x , x , · · · , x N } is minimal. Clearly, in the case Γ being a circle, the optimalredistribution is represented by a regular N -polygon. However, it should be obviousthat, e.g. for an oval, the optimal redistribution has to take curvature of the curveinto account. Interestingly enough, we will prove that the length/the area discrepancyminimizing redistribution X = { x , x , · · · , x N } , for N → ∞ , is closely related to thecurvature adjusted redistribution discussed in previous section with the shape functions ϕ ( k ) = | k | / for the length discrepancy ,ϕ ( k ) = | k | / for the area discrepancy . Suppose that a given smooth closed curve Γ is parameterized by a smooth function x :[0 , R . Denote by X = { x , x , · · · , x N } the set of grid points of Γ such that x i = x ( u i ) where u i = ih and h = 1 /N . 9 he length discrepancy. Let L = R Γ ds be the length of Γ and L ( X ) = N − X j =0 | x j +1 − x j | be the length of a polygon with vertices { x , x , · · · , x N } . Here we have identified x N + i = x i since Γ is assumed to be a closed curve. Clearly, L ( X ) ≤ L . Our goal is to find condi-tions under which the parameterization x ( · ) yields the minimizer X = { x , x , · · · , x N } of the problemmin X ⊂ Γ ( L − L ( X )) . (4.1)Since L ( X ) = i − X j =0 | x j +1 − x j | + | x i − x i − | + | x i +1 − x i | + N − X j = i +1 | x j +1 − x j | we obtain the following expression for the derivative of L with respect to x i in the direction y ∈ R : L ′ x i ( X ) y = − e n i . y , e n i = x i +1 − x i | x i +1 − x i | − x i − x i − | x i − x i − | Recall that the above minimization problem (4.1) is subject to the constraint x i ∈ Γ foreach i = 1 , , · · · , N . Therefore, the first order necessary condition for the constrainedminimizer X of (4.1) reads as follows: L ′ x i ( X ) y = 0 for any y ∼ t i , i = 1 , , · · · , N, where the symbol y ∼ t i means that the vector y is collinear with the unit tangent vector t i to Γ at the point x i ∈ Γ. Hence, the necessary condition for a set X = { x , · · · , x N } ∈ R × N of points belonging to Γ to be a minimizer of the functional X L − L ( X ) can berewritten as: e n i ⊥ t i , i = 1 , , · · · , N. (4.2)A graphical description of the necessary condition is depicted in Figure 4.1 (a).Next we will express the necessary condition (4.2) in terms of the parameterization x ( · ) of the curve Γ. The Taylor expansion of x ( · ) at x i = x ( ih ) yields x i ± − x i = ± ∂ u x ( u i ) h + ∂ u x ( u i ) h ± ∂ u x ( u i ) h ∂ u x ( u i ) h 24 + O ( h ) , as h → + . Derivatives ∂ ju x at u i = ih can be decomposed as follows: ∂ ju x ( u i ) = b j n i + a j t i , j = 1 , , · · · , n i and t i are, respectively, the unit inward normal and tangent vector to the curveΓ at a point x i ∈ Γ. Then a = g , b = 0 where g = g i = | ∂ u x ( u i ) | is the local lengthat the point x i ∈ Γ. Furthermore, using the Fren´et formulae g − ∂ u t = ∂ s t = k n and g − ∂ u n = ∂ s n = − k t , we obtain the recurrent relations a = g, b = 0 , a j +1 = ∂ u a j − gkb j , b j +1 = ∂ u b j + gka j . (4.3)for j ≥ k = k i is the curvature of Γ at x i ∈ Γ.Let us define Φ( h ) and Ψ( h ) as the following auxiliary functionsΦ( h ) = ∞ X j =1 a j j ! h j − , Ψ( h ) = ∞ X j =1 b j j ! h j − . Using these functions, the forward and backward difference at x i can be expressed as x i ± − x i ± h = Φ( ± h ) t i + Ψ( ± h ) n i , | x i ± − x i | h = p Φ( ± h ) + Ψ( ± h ) . Then the necessary condition (4.2) can be rewritten as: F ( h ) = 0 , where h = 1 /N . The function F is defined by F ( h ) = Φ( h ) p Φ( h ) + Ψ( h ) − Φ( − h ) p Φ( − h ) + Ψ( − h ) . Due to the symmetry F ( − h ) = − F ( h ) we have F (0) = F ′′ (0) = 0. Moreover, F ′ (0) = 0because of the fact that b = 0 and so Ψ(0) = 0. The leading order term of Taylor’sexpansion of F ( h ) is therefore given as: F ( h ) = 16 F ′′′ (0) h + O ( h )as h → + . Calculation of F ′′′ (0) can be simplified if one uses the substitution ξ ( h ) :=Ψ( h ) / Φ( h ). Notice that Φ( h ) > | h | ≪ a = g > 0. Then F ′′′ (0) = − ξ ′ (0) ξ ′′ (0). Using (4.3) we finally obtain F ′′′ (0) = 3 b a (cid:20) a − a b b (cid:21) . Since F ( h ) ≡ a − a b b = 0 , which has to be satisfied by the minimizer X in the limit h → + , i.e. for N → ∞ . Nowcalculating the corresponding terms a , a , b and b from the recurrent relations (4.3) weobtain a = g, a = ∂ u g, b = kg , b = ∂ u ( kg ) + gk∂ u g. ∂ u gg − ∂ u ( kg ) + gk∂ u gkg = 0 on Γ = { x ( u ); u ∈ [0 , } , which is clearly equivalent to the statement: | k | g / = constant , i.e. | k | / g = constant (4.4)on the curve Γ = { x ( u ); u ∈ [0 , } . It corresponds to the curvature adjusted tangentialvelocity with ϕ ( k ) = | k | / and the extended relative local length (see (3.1)) r ϕ ( u ) = g ( u ) L ϕ ( k ( u )) h ϕ ( k ( · )) i ≡ , for each u ∈ [0 , . In other words, such a curvature adjusted tangential velocity yields (in the limit of numberof grid points tending to infinity) the best possible approximation of the evolving familyof planar curves as far as the minimization of the error between the length of a curve andthe length of its polygonal approximation is concerned.Figure 3.1 (e) displays an example of the length discrepancy minimizing redistributionfor the case when ϕ ( k ) = | k | / . The values of length defect ∆ L := 1 − L ( X ) /L for several X ’s corresponding the cases inFigure 3.1 and the length optimal case are calculated as in Table 4.1. The area discrepancy. Similarly as in the case of the length functional we canask the question what is an asymptotically optimal redistribution of points on a givencurve Γ such that the discrepancy in areas enclosed by the curve Γ and its polygonalapproximation is minimal.The area A enclosed by the curve Γ is given by A = R Γ det( x , ∂ s x ) ds . The areaenclosed by the closed polygonal curve with vertices X = { x , x , · · · , x N } is given by A ( X ) = 12 N − X j =0 det( x j , x j +1 − x j ) , where x = x N . Clearly, A ( X ) = P N − j =0 det( x j , x j +1 ). The first order necessary condi-tion for a set X = { x , x , · · · , x N } of vertices to be a minimizer of the area discrepancyfunctionalmin X ⊂ Γ ( A − A ( X )) , where the vertices are constrained to belong to the curve Γ reads as follows: A ′ x i ( X ) y = 0 for any y ∼ t i , i = 1 , , · · · , N. 12e have A ′ x i ( X ) y = 12 (det( x i − , y ) + det( y , x i +1 )) = − 12 det( e t i , y ) , e t i = x i +1 − x i − for any direction y ∈ R . Hence, the necessary condition for a set X = { x , · · · , x N } ∈ R × N of points belonging to Γ to be a minimizer of the functional X ( A − A ( X )) canbe rewritten as: e t i // t i ⇔ e t i ⊥ n i , i = 1 , , · · · , N. (4.5)A graphical description of the necessary condition is depicted in Figure 4.1 (b).Now, calculating the difference e t i by means of the Taylor series and using the facts:det( t i , n i ) = 1 and det( t i , t i ) = 0 we obtain that condition (4.5) can be restated as:13 b h + O ( h ) = 0as h → + . Therefore, in the limit h = 1 /N → 0, we conclude that b = 0. Since b = ∂ u ( kg ) + gk∂ u g = kg ∂ u (ln( | k | g ) + ln g ), we obtain g | k | = constant on Γ, i.e. g | k | / = constant . It means that the area discrepancy curvature adjusted tangentialvelocity has the shape function ϕ ( k ) = | k | / . Figure 3.1 (f) displays an example of the area discrepancy minimizing redistributionfor the case when ϕ ( k ) = | k | / . The values of area defect ∆ A := 1 − A ( X ) /A for several X ’s corresponding the cases in Figure 3.1 and the area optimal case are calculated as inTable 4.1. The purpose of this section is to construct a numerical approximation scheme for solvingthe governing equation (2.2) for the position vector and the tangential velocity equation(3.3) satisfying renormalization constraint (3.4). Scheme. For a given initial N -sided polygonal curve P = S Ni =1 S i , find a family of N -sided polygonal curves {P j } j =1 , , ··· , P j = S Ni =1 S ji , where S ji = [ x ji − , x ji ] is the i -thedge with x j = x jN for j = 0 , , , · · · . The initial polygon P is an approximation ofΓ satisfying { x i } Ni =1 ⊂ P ∩ Γ , and P j is an approximation of Γ t at the time t = t j ,where t j = jτ is the j -th discrete time ( j = 0 , , , · · · ) if we use a fixed time increment τ > 0, or t j = P j − l =0 τ l ( j = 1 , , · · · ; t = 0) if we use adaptive time increments τ l > , l = 0 , · · · , j − 1. The updated curve {P j +1 } is determined from the data {P j } at theprevious time step by using discretization in space and time as follows. Discretization in space. Let P = S Ni =1 S i be an N -sided polygonal curve, where S i = [ x i − , x i ] is the i -th edge and x i is the i -th vertex ( i = 1 , , · · · , N ; x = x N ).13e denote the length of S i by r i = | x i − x i − | . The i -th unit tangent vector t i canbe defined as t i = ( x i − x i − ) /r i . Then the i -th unit tangent angle ν i is obtained from t i = (cos ν i , sin ν i ) T in the following way: firstly, from t = ( t , t ), we obtain ν =2 π − arccos( t ) if t < ν = arccos( t ) if t ≥ 0. Secondly, for i = 1 , , · · · , N wesuccessively compute ν i +1 from ν i : ν i +1 = ν i + arcsin( D ) if I > ,ν i + arccos( I ) if D > ,ν i − arccos( I ) otherwise , where D = det( t i , t i +1 ) , I = t i . t i +1 . Finally, we obtain ν = ν − ( ν N +1 − ν N ) and ν N +2 = ν N +1 + ( ν − ν ).In order to derive a discrete numerical scheme, we follow the flowing finite volumemethod adopted for curve evolutionary problems as it was proposed by Mikula et al.in [17, 19]. Let us introduce the “dual” volume S ∗ i = [ x ∗ i , x i ] ∪ [ x i , x ∗ i +1 ] of S i , where x ∗ i = ( x i − + x i ) / i = 1 , , · · · , N ; x ∗ N +1 = x ∗ ). We define the i -th unit tangent angleof S ∗ i by ν ∗ i = ( ν i + ν i +1 ) / 2. The i -th curvature k i has the constant value on S i , which isobtained from integration of k = ∂ s ν over S i with respect to s : Z S i k ds = k i Z S i ds = k i r i , Z S i k ds = Z S i ∂ s ν ds = [ ν ] x i x i − = ν ∗ i − ν ∗ i − . Hereafter, R S i F ds means R s i s i − F ds for arc-length s i satisfying x i = x ( s i , · ). Thus we have k i = ( ∂ s ν ∗ ) i , where ( ∂ s F ) i = ( F i − F i − ) /r i . The i -th curvature k ∗ i at x i can be defined as k ∗ i = ( k i +1 + k i ) / S ∗ i .Next we discretize equation (3.3) for the tangential velocity α : ∂ s ( ϕα ) = h f ih ϕ i ϕ − f + (cid:18) Lg h ϕ i − ϕ (cid:19) ω, where f = ( ∂ s β + k β ) ϕ ′ ( k ) − kβϕ ( k ). Integrating the above equation over S i yields Z S i ∂ s ( ϕα ) ds = (cid:2) ϕ ( k ) α (cid:3) x i x i − = h f ih ϕ i Z S i ϕ ( k ) ds − Z S i f ds + (cid:18) L h ϕ i Z S i g ds − Z S i ϕ ( k ) ds (cid:19) ω. Then ψ i = ϕ ( k ∗ i ) α i − ϕ ( k ∗ i − ) α i − = h f ih ϕ i ϕ ( k i ) r i − f i r i + (cid:18) Lh ϕ i N − ϕ ( k i ) r i (cid:19) ω,f i = (cid:0) ( ∂ s ( ∂ s ∗ β )) i + k i β i (cid:1) ϕ ′ ( k i ) − k i β i ϕ ( k i ) , where β i = β ( x ∗ i , k i , t i ) is constant on S i ,( ∂ s ( ∂ s ∗ β )) i = ( ∂ s ∗ β ) i − ( ∂ s ∗ β ) i − r i = 1 r i (cid:2) ∂ s β (cid:3) x i x i − , ( ∂ s ∗ F ) i = F i +1 − F i r ∗ i , ∗ i = ( r i + r i +1 ) / S ∗ i , and L = N X i =1 r i is the total length of P . The averages h f i and h ϕ i are approximated as: h f i = 1 L N X i =1 f i r i , h ϕ i = 1 L N X i =1 ϕ ( k i ) r i . To determine { α i } Ni =1 , we have to take account of the renormalization constraint h ϕα i = 0. It can be discretized as: h ϕα i = 1 L N X i =1 ϕ ( k ∗ i ) α i r ∗ i = 0 . Notice that L = P Ni =1 r i = P Ni =1 r ∗ i . We define a partial sum of { ψ i } byΨ i = i X l =2 ψ l ( i = 2 , , · · · , N ) , Ψ = 0 . With this notation we obtain ϕ ( k ∗ i ) α i r ∗ i = ϕ ( k ∗ ) α r ∗ i + Ψ i r ∗ i . Summing the above terms yields N X i =1 ϕ ( k ∗ i ) α i r ∗ i = ϕ ( k ∗ ) α L + N X i =1 Ψ i r ∗ i , L = N X i =1 r ∗ i . Hence we obtain α = − L ϕ ( k ∗ ) N X i =2 Ψ i r ∗ i , α i = 1 ϕ ( k ∗ i ) ( ϕ ( k ∗ ) α + Ψ i ) ( i = 2 , , · · · , N ) . Discretization in time. The semidiscretized (in space) evolution equation (2.2) hasthe form: ∂ t x i = w ∗ i ( ∂ s ∗ t ) i + α i ( ∂ s ∗ x ∗ ) i + F ∗ i n ∗ i , w ∗ i = w ( x i , ν ∗ i , k ∗ i ) F ∗ i = F ( x i , ν ∗ i )for i = 1 , , · · · , N , which follows from integration of ∂ t x = w∂ s t + α∂ s x + F n (see (2.2))over the volume S ∗ i . The right-hand side can be discretized as follows: ∂ t x i = w ∗ i r ∗ i (cid:18) x i +1 − x i r i +1 − x i − x i − r i (cid:19) + α i r ∗ i ( x i +1 − x i − ) + F ∗ i n ∗ i . 15n our approach, we use the semi-implicit numerical scheme for discretization in time: x j +1 i − x ji τ = a − x j +1 i − + a x j +1 i + a + x j +1 i +1 + F ∗ ji n ∗ ji ,a − = br ji − a, a = − ( a − + a + ) , a + = br ji +1 + a, a = α ji r ∗ ji , b = w ∗ ji r ∗ ji , where F ji is the i -th data of an N -sided polygonal curve P j ≈ Γ( t j ) for i = 1 , , · · · , N . Inother words, the following tridiagonal matrix has to be solved in order to obtain solutionin the new j + 1 time level: − a − τ x j +1 i − + (1 − a τ ) x j +1 i − a + τ x j +1 i +1 = x ji + F ∗ ji n ∗ ji τ ( i = 1 , , · · · , N )subject to periodic boundary conditions x j +10 = x j +1 N , x j +1 N +1 = x j +11 . We note that thetridiagonal matrix is strictly diagonally dominant provided that τ is small enough. Thetime step τ can be also chosen adaptively using the following expression τ j = r j min λ ) (cid:18) w j ∗ max r j min + | α j | max (cid:19) − , (5.1)where λ > F min = min ≤ i ≤ N F i , F max = max ≤ i ≤ N F i and | F | max = max ≤ i ≤ N | F i | . Withthis choice of τ j the solution { x j +1 i } Ni =1 exists and it is unique. This section is devoted to presentation of various numerical experiments. In all figures, forthe prescribed b τ > 0, we plot every µ b τ discrete time step using discrete points representingthe evolving curve. In every 3 µ b τ time step, we plot a polygonal curve connecting thosepoints, where µ = [[ T / b τ ] / x ] is the integer part of x ), and T is approximation of thefinal time which computed as follows. Let A t and L t be the enclosed area and the totallength of P j at time t = t j . In Figures 6.7, 6.8 and 6.9, T is the smallest discrete time t for which both conditions |A t / A t − τ − | < δ and |L t / L t − τ − | < δ are satisfied. In allother figures, T is the smallest discrete time t such that A t / A < δ . In figure captions,we provide the number of grid points N , ε ∈ (0 , 1) for ϕ ( k ) = 1 − ε + ε √ − ε + εk , κ and κ for the relaxation function ω ( t ). In all examples, the initial discretization is givenby the uniform N -division of the u -range [0 , Initial test curves. As an initial test examples we use a circle, an ellipse and thefollowing initial curves x ( u, 0) = ( x ( u ) , x ( u )) T depicted in Figure 6.1 and parameterizedby x ( u ) = cos z, x ( u ) = 0 . z + sin x + x , u ∈ [0 , , where x = sin(3 z ) sin z and z = 2 πu (left), and x ( u ) = 1 . z, x ( u ) = 1 . . z + 0 . x + 0 . x + 0 . x ) , u ∈ [0 , , x = sin(3 z ) sin z , x = 2 x , x = 3 e − x and z = 2 πu (right) (cf. [15]). Classical curvature flows. According to the classical convexification theory and theasymptotic behavior derived by M. A. Grayson, M. Gage and R. S. Hamilton, in the casewhere the geometric equation is given by β = k , any embedded curve becomes convex infinite time [9], and any convex curve shrinks to a single point. Its asymptotic shape is acircle [8]. See Figure 6.2. Affine curvature flows. Convexification also holds true in the case where the evo-lution law is the so-called affine scale space normal velocity β = k / (cf. [24]). In thiscase, convexification was proved and the asymptotic behavior was derived by G. Sapiroand A. Tannenbaum. They showed that any embedded curve shrinks to a single pointwith an ellipse as the asymptotic shape [24]. See Figure 6.3. Experimental order of convergence (EOC). In the case of the normal velocity β = k / , any ellipse shrinks to a point homothetically. It means that an ellipse is aself-similar solution to (1.1) (see Figure 6.4). By using the explicit self-similar solution,the so-called experimental order of convergence (EOC) can be obtained in the followingway.When the position vector is described by x ( u, t ) = η ( t ) z ( u ) with η (0) = 1, and whenthe initial curve is an ellipse x ( u, 0) = z ( u ) = ( a cos(2 πu ) , b sin(2 πu )) T with a, b > 0, forthe curvature we obtain k ( u, t ) = ab η ( t ) − ζ ( u ) − / , ζ ( u ) = a sin (2 πu ) + b cos (2 πu ).In the case when β = k / , ∂ t x · n = ( ab ) / η ( t ) − / ζ ( u ) − / holds. On the other hand,we have ∂ t x · n = ∂ t η ( t ) z · n = − ab ∂ t η ( t ) ζ ( u ) − / . Hence we obtain the rate η ( t ) =(1 − ( ab ) − / t ) / and the extinction time T = ( ab ) / , which is determined by η ( T ) = 0.Numerical errors at t = t j with the number of points N can be defined aserr jp ( N ) = max ≤ i ≤ N (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( x j i ) ( a η ( t j )) + ( x j i ) ( b η ( t j )) − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , if p = ∞ , N X ≤ i ≤ N (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( x j i ) ( a η ( t j )) + ( x j i ) ( b η ( t j )) − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) p ! /p , if 1 ≤ p < ∞ , where x ji = ( x j i , x j i ) T is the i -th grid point. Therefore we can define the L q (cid:0) (0 , t M ) , L p (0 , (cid:1) error such as E p,q ( N ) = max ≤ j ≤ M err jp , if q = ∞ , M X ≤ j ≤ M (cid:0) err jp (cid:1) q ! /q , if 1 ≤ q < ∞ , where t M < T . If E p,q ( N ) ≈ N − µ holds, then E p,q ( N/ ≈ µ N − µ also holds. Fromthese approximations, we may expect µ ≈ log (cid:0) E p,q ( N/ /E p,q ( N ) (cid:1) . We therefore usethe right hand side as an indicator of the numerical convergence, i.e., the EOC from the17 q (cid:0) (0 , t M ) , L p (0 , (cid:1) error is defined asEOC p,q ( N ) = log E p,q ( N/ E p,q ( N ) . Tables 6.1, 6.2, 6.3, 6.4 indicate the L q (cid:0) (0 , t M ) , L p (0 , (cid:1) errors E p,q ( N ) and EOCEOC p,q ( N ) for N = 2 , , , , , ( p, q ) ∈ { , , ∞} , and ε = 0 , . , . , . 9. Parametersare τ = 0 . N − and κ = κ = 100, and sample data are used at the time t j = 1 . j/M with j = 0 , , · · · , M = 200. Ellipses with the axes ratio 3 : 1 are the same as in Figure 6.4.In this case, the extinction time T = ( ab ) / = 1 . · · · > t M = 1 . 5. Since the governingsystem is of the parabolic type, we have used the time step τ ≈ N − . From Tables 6.1,6.2, 6.3, 6.4, we can observe EOC ≈ ε ’s. It means thatthese four tangential redistribution methods seem to be almost equally effective in viewof the EOC. The length and the area discrepancy. An ellipse is a shrinking self-similar solutionto (1.1) in the case of the normal velocity β = w ( ν ) k , as well as the affine curvature flows β = k / . Here the weight is w ( ν ) = a b T ( a sin ν + b cos ν ) , when the axes ratio of ellipse is a : b and the extinction time is T > 0. Similarly as inthe previous subsection EOC, one can obtain the position vector of a shrinking ellipse by x ( u, t ) = η ( t ) z ( u ) with the scaling function η ( t ) = p − t/T ( η (0) = 1, η ( T ) = 0) andthe initial ellipse x ( u, 0) = z ( u ) = ( a cos(2 πu ) , b sin(2 πu )) T . By using this exact solution,we will calculate the numerical test of the length and the area discrepancy in the followingway.Numerical discrepancy of the length and area defects at the time t = t j can be definedas ∆ jL = (cid:12)(cid:12)(cid:12)(cid:12) − L t j L t j (cid:12)(cid:12)(cid:12)(cid:12) , ∆ jA = (cid:12)(cid:12)(cid:12)(cid:12) − A t j A t j (cid:12)(cid:12)(cid:12)(cid:12) , respectively. Here L t = η ( t ) L and A t = η ( t ) A are the length and area of Γ t , and L t and A t are the length and area of P j at t = t j , respectively. Therefore we can define the L q (0 , t M ) numerical discrepancy of the length and area as follows∆ ∗ ,q = max ≤ j A,q for q = 1 , , ∞ . Parameters are: N = 100, τ =0 . N − and κ = κ = 100, the extinction time is T = 1, sample data are used at the18ime t j = T j/M with j = 0 , , · · · , M = 200, and the axes ratio of the initial ellipse Γ are a = 3 : b = 1. As for the shape functions we chose: ϕ ( k ) = 1 − ε + ε √ − ε + εk with ε = 0 , . , . , ϕ ( k ) = | k | / and ϕ ( k ) = | k | / . In Table 6.5, we can observe thatthe shape function ϕ ( k ) = | k | attains the minimum value in each ∆ ∗ ,q for ∗ = L, A and q = 1 , , ∞ . Weighted curvature flows. Asymptotic behavior of solutions to the weighted cur-vature flow β = w ( ν ) k is related to self-similar shrinking solutions, which need not beunique. For details we refer to the paper by Yagisita [28] and references therein. InFigure 6.5 we show evolving family of plane curves with β = w ( ν ) k . Forced curvature flows. The theory of Grayson, Gage and Hamilton was generalizedto the case of an anisotropic curvature driven flow by K. S. Chou and X. P. Zhu in [3].They considered the case where the evolution law is given by β = w ( ν ) k + F ( ν ) with w ( ν ) = σ ( ν ) + ∂ ν σ ( ν ), where σ is a given anisotropy function, σ ( ν + π ) = σ ( ν ), and F ( ν + π ) = − F ( ν ). They proved that any embedded curve becomes convex in finitetime [4], and any convex curve shrinks to a single point with the asymptotic shape beinga self-similar solution to β = w ( ν ) k , i.e., the Wulff shape of σ [3]. In Figure 6.6 weillustrate the convexification theory and the asymptotic behavior of evolving plane curvesdue to K. S. Chou and X. P. Zhu. Loss of convexity phenomenon. In the case when the normal velocity is givenby β = k + F ( x , ν ), an initially convex curve may loose its convexity in finite time fora special choice of the external force F . In Figure 6.7 we present such examples with F presented by Nakamura et al. in [21]. It is worth noting that the usual crystallinecurvature flow is unable to capture this convexity-breaking phenomena. Image segmentation by using a gradient flow. Let γ ( x ) > t . If γ is differentiable, then the L gradient flow of thefollowing energy E (Γ t ) = Z Γ t γ ( x ) ds is realized by the geometric equation (1.1) of the form β = γ ( x ) k − ∇ γ ( x ) · n . Such agradient flow can be successfully utilized in various image segmentation problems. Forexample, let an image intensity function be denoted by I : R ⊃ Ω → [0 , I ( x ) is a piecewise constant function on each pixel. Here I = 0 ( I = 1) corresponds tothe black (white) color and I ∈ (0 , 1) corresponds to a scale of gray colors. For simplicity,we assume that our target figures are drawn in white color with the black background.Then edges of the image correspond to regions where the gradient |∇ I ( x ) | is sufficientlylarge. Let us introduce an auxiliary density function γ ( x ) = f ( |∇ I ( x ) | ) where f is asmooth edge detector function such as f ( s ) = 1 / (1 + s ) or f ( s ) = e − s . Hence thesolution curve Γ t of β = γ ( x ) k − ∇ γ ( x ) · n has the tendency to minimize the energy E (Γ t ). In other words, it moves toward the edge in the image on which |∇ I ( x ) | is large.19his is a fundamental idea of image segmentation, and it has developed to a sophisticatednumerical scheme [17, 18]. An example of such a curve evolution converging to an edgein the given images is depicted in Figure 6.8. Image segmentation for sharp images. If contrast of the target image is sufficientlyhigh, a simpler scheme described in [2] can be used. We consider the geometric flow β = k + F and define the forcing term F ( x ) as follows: F ( x ) = F max − ( F max − F min ) I ( x ) ( x ∈ Ω) , where F max > F min < /F is equivalent to the minimal radius the curve canattain. Choosing an appropriate F , we can make the final shape be rounded, or we canprevent the curve from passing through the outline. A segmentation of a given sharpimage by means of a plane curve evolution is shown in Figure 6.9. In this paper, we proposed and analyzed a new class of tangential velocities by which wecan control tangential motion and grid point redistribution of plane curves that evolve withthe normal velocity depending on a general function of the curvature, tangential angle andthe position vector. The curvature adjusted tangential velocity may not only distributegrid points uniformly along the curve but also produce a desirable concentration and/ordispersion depending on the curvature. We also demonstrated that curvature adjustedtangential redistribution yields the best possible approximation of the evolving familyof planar curves as far as the minimization of the error between the length and areaof a stationary curve and that of its polygonal approximation is concerned. Numericalexperiments based on semi-implicit numerical flowing finite volume method confirmedcapability of our new method. Acknowledgments. The authors would like to express our gratitude to the anonymousreferees for their comments and suggestions. References [1] J. W. Barrett, H. Garcke and R. N¨urnberg, A parametric finite element method forfourth order geometric evolution equations , J. Comput. Phys. (2007), 441–467.[2] M. Beneˇs, M. Kimura, P. Pauˇs, D. ˇSevˇcoviˇc, T. Tsujikawa and S. Yazaki, “Appli-cation of a curvature adjusted method in image segmentation”, Bull. Inst. Math.Acad. Sin. (New Series) (4) (2008), 509–523.203] K. S. Chou and X. P. Zhu, Anisotropic flows for convex plane curves , Duke Math.J. (1999), 579–619.[4] K. S. Chou and X. P. Zhu, A convexity theorem for anisotropic flows of plane curves ,Indiana Univ. Math. J. (1999), 139–154.[5] K. Deckelnick, Weak solutions of the curve shortening flow , Calc. Var. Partial Dif-ferential Equations (1997), 489–510.[6] G. Dziuk, Convergence of a semi discrete scheme for the curve shortening flow ,Math. Models Methods Appl. Sci. (1994), 589–606.[7] C. L. Epstein and M. Gage, The curve shortening flow , Wave motion: theory,modelling, and computation (Berkeley, Calif., 1986), Math. Sci. Res. Inst. Publ. ,Springer, New York (1987), 15–59.[8] M. Gage and R. S. Hamilton, The heat equation shrinking convex plane curves , J.Diff. Geom. (1986), 69–96.[9] M. A. Grayson, The heat equation shrinks emmbedded plane curves to round points ,J. Diff. Geom. (1987), 285–314.[10] T. Y. Hou, J. S. Lowengrub and M. J. Shelley, Removing the stiffness from interfacialflows with surface tension , J. Comput. Phys. (1994), 312–338.[11] T. Y. Hou, I. Klapper and H. Si, Removing the stiffness of curvature in computing -D filaments , J. Comput. Phys. (1998), 628–664.[12] M. Kimura, Accurate numerical scheme for the flow by curvature , Appl. Math.Letters (1994), 69–73.[13] M. Kimura, Numerical analysis for moving boundary problems using the boundarytracking method , Japan J. Indust. Appl. Math. (1997), 373–398.[14] K. Mikula and D. ˇSevˇcoviˇc, Solution of nonlinearly curvature driven evolution ofplane curves , Appl. Numer. Math. (1999), 191–207.[15] K. Mikula and D. ˇSevˇcoviˇc, Evolution of plane curves driven by a nonlinear functionof curvature and anisotropy , SIAM J. Appl. Math. (2001), 1473–1501.[16] K. Mikula and D. ˇSevˇcoviˇc, A direct method for solving an anisotropic mean cur-vature flow of plane curves with an external force , Math. Methods Appl. Sci. (2004), 1545–1565.[17] K. Mikula and D. ˇSevˇcoviˇc, Computational and qualitative aspects of evolution ofcurves driven by curvature and external force , Comput. Vis. Sci. (2004), 211–225.2118] K. Mikula and D. ˇSevˇcoviˇc, Evolution of curves on a surface driven by the geodesiccurvature and external force , Appl. Anal. (2006), 345–362.[19] K. Mikula, D. ˇSevˇcoviˇc and M. Balaˇzovjech, A simple, fast and stabilized flow-ing finite volume method for solving general curve evolution equations , Commun.Comput. Phys. (1) (2010), 195–211.[20] S. Morigi, Geometric surface evolution with tangential contribution , J. Comput.Appl. Math. (5) (2010), 1277–1287.[21] K.-I. Nakamura, H. Matano, D. Hilhorst and R. Sch¨atzle, Singular limit of areaction-diffusion equation with a spatially inhomogeneous reaction term , J. Statist.Phys. (1999), 1165–1185.[22] P. Pauˇs and M. Beneˇs, Direct approach to mean-curvature flow with topologicalchanges , Kybernetika (4) (2009), 591–604.[23] P. Pauˇs, M. Beneˇs and J. Kratochv´ıl, Discrete dislocation dynamics and mean cur-vature flow , Numerical Mathematics and Advanced Applications 2009, Proceedingsof ENUMATH 2009, Springer-Berlin-Heidelberg, Part 2, (2010), 721–728.[24] G. Sapiro and A. Tannenbaum, On affine plane curve evolution , J. Funct. Anal. (1994), 79–120.[25] D. ˇSevˇcoviˇc and S. Yazaki, On a motion of plane curves with a curvature adjustedtangential velocity , submitted (arXiv:0711.2568).[26] J. A. Sethian, Level Set Methods and Fast Marching Methods: Evolving Interfaces inComputational Geometry, Fluid Mechanics, Computer Vision, and Material Science ,Cambridge University Press, New York (1999).[27] T. K. Ushijima and S. Yazaki, Convergence of a crystalline algorithm for the motionof a closed convex curve by a power of curvature V = K α , SIAM Journal onNumerical Analysis (2000), 500–522.[28] H. Yagisita, Non-uniqueness of self-similar shrinking curves for an anisotropic cur-vature flow , Calc. Var. (2006), 49–55.[29] S. Yazaki, On the tangential velocity arising in a crystalline approximation of evolv-ing plane curves , Kybernetika (2007), 913–918.22a) (b)(c) (d)(e) (f)Figure 3.1: Various redistributions of N = 12 grid points along the ellipse: (a) uniform re-distribution ( ε = 0), (b) curvature adjusted redistribution ( ε = 0 . ε = 1), (d) crystalline redistribution (corresponding to ε = 1), (e) the lengthdiscrepancy minimizing curvature adjusted redistribution with ϕ ( k ) = | k | / , and (f) the areadiscrepancy minimizing curvature adjusted redistribution with ϕ ( k ) = | k | / .(a) uniform (b) curv. adj. (c) crystalline (e) the length (f) the area X ε = 0 ε = 0 . ε = 1 optimal optimal∆ L A Table 4.1: The values of length defect ∆ L and area defect ∆ A for the cases in Figure (a),(b), (c), (e) the length optimal case ϕ ( k ) = | k | / , and (f) the area optimal case ϕ ( k ) = | k | / .Underlined values are the minimum in each defect. x i+1 x i-1 x i t~ i n Γ i x i+1 x i-1 x i t i n Γ ~ i t (a) (b)Figure 4.1: A graphical illustrations of the necessary conditions at a point x i for maximizingof the polygonal length discrepancy L ( X ) (a) and the polygonal enclosed area A ( X ) (b) of adiscrete approximation X = { x , · · · , x N } of a smooth curve Γ, respectively. Figure 6.1: Initial curves with grid points corresponding to a uniform division of u ∈ [0 , 1] with N = 100 points. (a) (b) (c) (d)Figure 6.2: Numerical examples with the normal velocity β = k . (a) and (c) are evolving curves,respectively. (b) and (d) are the corresponding final magnified curves. Numerical parameters: N = 100, ε = 0 . κ = κ = 100, τ = 0 . N − , δ = 10 τ , and b τ = 0 . Numerical examples with affine the space scale normal velocity β = k / : (a) and(c) are evolving curves, and (b) and (d) are the corresponding final magnified curves. In bothcases, parameters were chosen: N = 100, ε = 0 . κ = κ = 100, τ = 0 . N − , δ = 10 τ , and b τ = 0 . Figure 6.4: From (far left) to (far right), ε = 0 , . , . , . 9. Upper figures are evolving ellipsesstarting from the same ellipse with the axes ratio 3 : 1. Lower figures are the final magnifiedcurves. Numerical parameters: N = 128, κ = κ = 100, τ = 0 . N − , δ = 10 τ , and b τ = 0 . p E p, ( N ) EOC p, ( N ) E p, ( N ) EOC p, ( N ) E p, ∞ ( N ) EOC p, ∞ ( N )16 1 0.5380220 0.5509331 0.63363002 1.0534598 1.0704854 1.2142337 ∞ ∞ ∞ ∞ ∞ Table 6.1: ε = 0 N p E p, ( N ) EOC p, ( N ) E p, ( N ) EOC p, ( N ) E p, ∞ ( N ) EOC p, ∞ ( N )16 1 0.5290974 0.5419132 0.62550062 1.0325342 1.0494491 1.1963774 ∞ ∞ ∞ ∞ ∞ Table 6.2: ε = 0 . p E p, ( N ) EOC p, ( N ) E p, ( N ) EOC p, ( N ) E p, ∞ ( N ) EOC p, ∞ ( N )16 1 0.4042607 0.4129555 0.49102152 0.7821867 0.7914808 0.9293693 ∞ ∞ ∞ ∞ ∞ Table 6.3: ε = 0 . N p E p, ( N ) EOC p, ( N ) E p, ( N ) EOC p, ( N ) E p, ∞ ( N ) EOC p, ∞ ( N )16 1 0.4408205 0.4417817 0.47586012 0.8949908 0.8959769 0.9555951 ∞ ∞ ∞ ∞ ∞ Table 6.4: ε = 0 . ( k ) ∆ L, ∆ L, ∆ L, ∞ ∆ A, ∆ A, ∆ A, ∞ ϕ ( k ) ≡ ϕ . ( k ) 0.013191 0.039820 0.417640 0.024203 0.089520 1.000038 ϕ . ( k ) 0.005216 0.016007 0.172409 0.009032 0.033209 0.370852 ϕ . ( k ) 0.001666 0.004468 0.047258 0.002119 0.008126 0.092089 ϕ ( k ) = | k | | k | / | k | / Table 6.5: Numerical discrepancy of the length ∆ L,q and the area ∆ A,q for q = 1 , , ∞ in variouschoice of ϕ ( k ), where ϕ ε ( k ) = 1 − ε + ε √ − ε + εk . (a) (b) (c) (d)(e) (f) (g) (h)Figure 6.5: (a), (b), (c), (d) are in the case where the weight is w ( ν ) = 1 − (7 / 9) cos(3 ν ), and(e), (f), (g), (h) are those where w ( ν ) = 1 − . ν − π/ N = 100, ε = 0 . κ = κ = 100, τ = 0 . N − , δ = 10 τ , and b τ = 0 . Numerical examples for the normal velocity β = w ( ν ) k + F ( ν ), w ( ν ) = 1 − . ν ), F ( ν ) = sin( ν ). (a) and (c) are evolving curves, respectively. (b) and (d) are thecorresponding final magnified curves. In both cases, parameters were chosen as: N = 100, ε = 0 . κ = κ = 100, τ = 0 . N − , δ = 10 τ , and b τ = 0 . σ ( ν ) = 1 +(0 . / 35) cos(6 ν ) is the unique solution of w = σ + ∂ ν σ , and (e) is the locus of the boundary ofthe Wulff shape of σ , which is given by σ ( − n ) + ( ∂ ν σ ) t . (a) (b) (c) (d)Figure 6.7: (a) and (b) correspond to the external force F = − pq sin( q (4 x + x ))( − x sin ν + x cos ν ) with p = 1 . q = 3 . 0. (c) and (d) correspond to the force F = 2 pqπ cos( qπ | x | ) x · n , p = 1 . q = 1 . 15. (a) (resp. (c)) indicates evolving curves starting from the unit circle (resp.ellipse with axes ratio 1 : 2), and (b) (resp. (d)) is the final curve at T . Parameters are N = 120, ε = 0 . κ = 100, κ = 0, τ = b τ = 0 . δ = 10 − . -1.5 -1 -0.5 0 0.5 1 1.5-1.5-1-0.5 0 0.5 1 1.5 (a) (b) (c)(e) (f)Figure 6.8: (a) indicates the original bitmap image faded shadow “Venus” in Ω = [ − . , . with the 600px resolution, (b) indicates image intensity function I ( x ) in 3D, and (c) is the grayscaled auxiliary function γ ( x ). (d) indicates evolving curves starting from circle with radius 1.5,and (e) is the final curve at T . Parameters were chosen as: N = 100, ε = 0 . κ = 100, κ = 0, b τ = 0 . δ = 10 − . Adaptive time step sizing ( ) with λ = 1 has been used. (a) indicates the original bitmap image “ya” (it is the Chinese character for arrow)in Ω = [ − . , . with the 600px resolution, (b) indicates evolving curves starting from circlewith radius 2, and (c) is the final curve at T . Parameters are N = 200, ε = 0 . κ = 100, κ = 0, b τ = 0 . δ = 10 − / F max = 30 and F min = − 30. Adaptive time step sizing ( )with λ = 1 has been used.= 1 has been used.