{'fraction_non_alphanumeric': 0.042853259022633046, 'fraction_numerical': 0.02086590090396248, 'mean_word_length': 4.393511730205279, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fairing of Discrete Planar Curves by Integrable Discrete Analogue of Euler's Elasticae

We construct a method to fair a given discrete planar curve by using the integrable discrete analogue of Euler's elastica, which is a discrete version of the approximation algorithm presented by D. Brander, et al. We first give a brief review of the integrable discrete analogue of Euler's elastica proposed by A. I. Bobenko and Yu. B. Suris, then we present a detailed account of the fairing algorithm, and we apply this method to an architectural problem of characterizing the keylines of Japanese handmade pantiles.


Introduction
The Euler's elastica (elastic curve) is a class of planar curves characterized as the solutions to the variational problem of minimizing the elastic energy under a certain boundary condition. It has been regarded as one of the most important class of planar curves because it is endowed with rich mathematical structure: exact solutions, integrability, geometry of elliptic curves, and so on, while it serves as a simple but realistic model of thin inextensible elastic rod (see, for example, [19,29]). Brander et al. [5] have proposed an algorithm to fair a given planar curve segment by an Euler's elastica, motivated mainly by the development of the robotic hot-blade cutting technology. In this work, motivated by a problem of architecture to characterize the keylines of Japanese handmade pantiles, where the curve data is obtained in the form of discrete point data, we aim to construct a fairing method of discrete planar curves by using the integrable discrete analogue of the Euler's elastica proposed by Bobenko and Suris [3], which is referred to as the discrete elastica in this paper.
This paper is organized as follows. We give a brief review of the Euler's elastica and the discrete elastica in Sections 2 and 3, collecting the information on variational formulations, exact solutions and continuum limits with proofs. We present a detailed account of the fairing of a given discrete planar curve by the discrete elastica in Section 4. Finally, application to the characterization of the keylines of Japanese handmade pantiles is discussed in Section 5. For various formulas of the Jacobi elliptic functions used in this paper, the readers may refer to [24], for example.

Euler's elastica
Let ( ) ∈ R 2 ( ∈ R) be an arc length parameterized planar curve. By definition, it holds that with respect to variations with fixed endpoints and tangent vectors at the endpoints, under the condition of preserving the total length. The Euler-Lagrange equation yields the following differential equations for the curvature and the angle function: where ∈ R is a constant. Moreover, the angle function satisfies where > 0 is a constant.
Derivation of equation (2.6) is given in various literatures such as [29]. Here, we show a concise derivation using the variation of the tangent vector. Consider the functional where the first term is the elastic energy, and the second and third terms correspond to the preservation and ( ) − (0), respectively, with = ( ) ∈ R and ∈ R 2 being the Lagrange multipliers. The variation of is calculated by using the Frenet formula (2.4) as , . (2.9) The first term is the boundary term that vanishes due to the boundary condition, and the second term gives the Euler-Lagrange equation, Taking the scalar product with and , we have respectively. Multiplying to both sides of the second equation of equation (2.11), we find that it is integrated to give where is a constant of integration. Then differentiating equation (2.7) and using equation (2.15) yields which is nothing but equation (2.6).
(2) Equations (2.6) and (2.7) can be seen as travelling-wave reductions of the (focusing) modified KdV equation and the sine-Gordon equation, respectively, where the former describes the integrable deformation of planar curves [8,13].
It is known that the differential equations (2.6) and (2.7) can be solved in terms of the Jacobi elliptic functions. This is well-known, but we present the solutions for the readers' convenience. In the literature, the solutions are often constructed from the the first integral of (2.6) given by where is a conserved quantity (constant). Here we present those solutions and verify them by the differential equations for the Jacobi elliptic functions.

Proposition 2.3.
The curvature and the angle function can be expressed in terms of the Jacobi elliptic functions as follows: = dn( , ),  [17,24]. Thus, there exists an analytic continuation for all Jacobi elliptic functions in the range ≥ 0. In this way, cases (i) and (ii) can be regarded as one. As it is known, the case (i) (resp. (ii)) yields the elastica without (resp. with) inflection points; see Figure 4. 3 Integrable discrete Euler's elastica

Basic framework of discrete planar curves
We first introduce the basic framework for discrete planar curves [12,14,20]. Let ∈ R 2 ( ∈ Z) be a discrete planar curve with | +1 − | = ℎ, where ℎ > 0 is a constant. We also assume det +1 − , − −1 ≠ 0, i.e., not three consecutive points are collinear. Then is called a discrete planar curve with segment length ℎ. We define the discrete tangent and normal vectors by respectively, where the discrete angle function Θ is the angle of measured from the horizontal axis in the counterclockwise direction. The discrete Frenet frame is defined by (2), (3.2) and the discrete Frenet formula by where = Θ − Θ −1 is the angle between two adjacent tangent vectors. Equation (3.3) is the discrete version of the Frenet formula (2.4); see Figure 2. The discrete curvature can be defined as the reciprocal of the radius of the osculating circle touching two adjacent segments at their midpoints [12] (see Figure 3), We note that the discrete Frenet formula (3.3) can be written in terms of as

Discrete Euler's elastica
The discrete Euler's elastica [3,4,11] is defined as a discrete planar curve with segment length ℎ that is a critical point of the functional with respect to variation with fixed endpoints and end edges. Note that means that the two functionals yield the same critical points. As mentioned in [4], can be regarded as a discrete analogue of the elastic energy (2.5), in a sense that 2/ℎ log(1 + ℎ 2 2 /4) is the potential energy of the bending force proportional to the discrete curvature at each vertex.
Taking into account the preservation of − −1 and − 0 by introducing the Lagrange multipliers ∈ R and ∈ R 2 , respectively, consider the functional Proposition 3.1. The Euler-Lagrange equation for the functional (3.7) is given by (3.8) or, equivalently, in terms of the discrete curvature by where ∈ R is a constant.
Proof. The variation of is calculated as follows: where is the boundary term given by which vanishes due to the boundary condition. Setting = 0 gives the Euler-Lagrange equation, which proves the first part. For the second part we use the discrete Frenet formula (3.5) written explicitly in terms of the tangent and normal vectors, From the discrete Frenet formula (3.3) and equation (3.4) we have that thus taking the inner product of both hand sides of equation (3.8) with gives − +1 = ℎ , . (3.14) Then, taking the inner product of the first equation in the discrete Frenet formula (3.12) with gives which implies that exists a constant ∈ R such that Finally taking the inner product of the second equation in the discrete Frenet formula (3.12) with , and using equations (3.14) and (3.16), yields which is exactly equation (3.9) with = 2 + ℎ 2 .
By using a technique similar to the one shown in [31], we construct explicit solutions to equation (3.9) from its discrete first integral, where ∈ R is a constant (conserved quantity). Here, we avoid the long computation required and simply present the solutions corresponding to equations (2.20) and (2.22), and verify them by using the addition formulas for the Jacobi elliptic functions.

Proposition 3.2.
Let ∈ R be a constant. Then the following functions satisfy equation (3.9). .

Remark 3.3.
(1) Comparing equations (2.20) and (3.19), we see that there exists a constant Ω such that = (Ω ). Indeed, we have (3.25) There is also a similar relationship between equations (2.22) and (3.20). This implies that the discrete curvature is an "exact discretization" of the smooth curvature ( ).
(3) Equation (3.9) is also known as the McMillan map, which is a special case of the Quispel-Roberts-Thompson (QRT) map solved by elliptic functions [26]. It can also be regarded as an autonomous version of a discrete Painlevé II equation [15,16,27] (4) It is known that position vectors of both smooth and discrete elasticae admit explicit formulas in terms of the elliptic theta functions [21,23].

Discrete Euler's elastica in terms of a potential function
Following [13], we say that is a potential function if it is such that In this context, the discrete curvature is written as

Proposition 3.4. Suppose that satisfies
where ∈ R is a constant. Then we have: (1) It holds that where Λ ∈ R is a constant.
(2) The discrete curvature satisfies equation Proof. The first statement is shown as follows. Multiplying sin( +1 − −1 4 ) to equation (3.30), using the product-to-sum formula, and rearranging terms gives: which implies (3.31). In order to show the second statement, we introduce for simplicity in the notation. We have that We expand equation ( where we used equation (3.31). From equation (3.36) and the sum-to-product formula, we have We use the following two expressions: On the one hand, from equation (3.30), On the other hand, where we used (3.38), which is rewritten as From the definition of discrete curvature, equation (3.41) is equivalent to equation (3.9) with = 2(1 − 2 )/Λ 2 , which proves the second statement.
Remark 3.5. Proposition 3.4 provides discrete analogues for equations (2.7) and (2.15). In fact, by putting We next present explicit solutions for equation (3.30). Part of Proposition 3.6 and 3.8 can be found in [30], in a slightly different context: in that work, the function is regarded as the angle function (here denoted as Θ ) instead of a potential function.

Remark 3.7.
(1) In case (i), the parameter in equation (3.9) is given by

46)
which implies that in equation (3.43) corresponds to in equation (3.19). In case (ii), so that in equation (3.44) corresponds to in equation (3.20). These correspondences can be verified directly by computing from equations (3.43) and (3.44), respectively.
We finally present the variational formulation for equation (3.30).
with respect to variations of the potential angle with fixed endpoints. (1) The discrete sine-Gordon equation [2,10,14], where , are lattice intervals. In fact, assuming that depends only on = + , equation (2) The discrete potential modified KdV equation [9], which describes the isoperimetric and equidistant deformation of discrete planar curves [13,20], is transformed to the discrete sine-Gordon equation (3.50) by → (−1) . In this sense, equation (3.30) can also be regarded as a reduction of the discrete potential modified KdV equation.

Approximation of discrete curves
In this section, we construct an algorithm to approximate a given discrete planar curve to a discrete elastica. Among the many possible discretizations for the elastica, the advantages of using the one shown in this work can be described as follows: First, the discrete elasticae are endowed with the same integrable structure as in their smooth counterpart, i.e., they possess several conserved quantities, can be obtained via a variational principle, and their explicit solutions are expressed in terms of Jacobi elliptic functions. Moreover, it is known that variational integrators have controlled error in their solutions [7,18,25]. In particular, the explicit expression for the discrete curvature is an "exact discretization" of the smooth curvature ( ) as discussed in Remark 3.3, and the potential function has the same functional shape as the smooth angle function ( ). From these observations, we expect this discretization to have good numerical properties.

General discrete Euler's elastica segment
To describe a general curve segment in the plane, we include the freedom of rotation in the equations. We do this by shifting the angle function and discrete angle function by a constant ∈ R in all the expressions. In particular, equation (3.30) goes to where we put = ℎ 2 /4, with > 0 a constant. From Proposition 3.6, we have where > 0, , ∈ R are constants. The parameter determines the shape of the elastica, the initial point, and is related with the length and point aggregation of the curve segment. Finally, from its starting point 0 ∈ R 2 , a discrete elastica segment is calculated recursively by Note that we can expand the sine and cosine in equation (4.4) and make use of equation (4.3) to obtain an explicit expression in terms of the Jacobi elliptic functions. We conclude that a general discrete elastica segment can be characterized by seven parameters: where 0 , 0 ∈ R are the two components of the initial point 0 . We write as ( ) to the discrete elastica with parameters .

Fairing process
Given a general discrete curve segment ∈ R 2 ( = 0, . . . , − 1), we look for a discrete elastica ( ) ∈ R 2 that is the closest, in a 2 -distance sense, to . Namely, we seek to find * such that * = arg min We solve this non-convex problem using the Interior Point Optimizer (IPOPT) package, that for our purpose can be seen as a gradient descent-like method for nonlinear optimizations [32]. For its implementation we need to compute the gradient and the Hessian of the objective function For the gradient, we have which is computed recursively from equation (4.4), Then, the only non-trivial derivatives we need to compute are 11) or equivalently, where we denoted ≡ sin − 2 and ≡ cos − 2 . Finally, we use equations (4.2), (4.3) and the derivatives of the Jacobi elliptic functions with respect to their argument and module to complete the computation. For the Hessian, we use a numerical quasi-Newton approximation, internally computed by the package.
The IPOPT method needs a starting pointˆ, that we refer as the initial guess. In the following subsection we describe the algorithm that we use to obtain the initial guess, which is a discrete analogue of the one provided in [5].

Initial parameters
The initial guess, that starts the IPOPT method, can be obtained in a numerically stable manner thanks to two geometric properties of the discrete elastica: Proposition 4.1 and Corollary 4.1.
Remarkably, these are geometrically equivalent to the same properties for the smooth elastica [5]. Let and define the projection of the curve onto I as 14) and the angle measured from I as Ψ = 2 + Θ − or, equivalently, such that  Note that, by putting 1 = cos and 2 = sin , equation (4.16) can be expressed as where , ∈ R are the two components of .

21)
where ∈ R is a constant.
Proof. From the definition of and Ψ in equations (4.14) and (4.15), respectively, we obtain Then, putting this into the discrete Frenet formula ( − −1 )/ℎ = ( + −1 )/2 gives where we used (4.16). After expanding the right hand side of the previous expression and then adding ± 2 , we conclude that there exists a constant ∈ R such that, for all , To estimate ( , , , ) we use some results from the smooth elastica to avoid unnecessary complexity in the discrete case. We use the following approximations: From equation We obtain an approximation of the parameter from these expressions. We see from equation  Hence, from equation (4.28) we obtain Then, we obtain the initial guessˆby solving equations (4.25) and (4.26) in the least-square sense and using several of the equations mentioned above. We proceed as follows: • (parameter ) From equation (4.25), compute (ˆ1,ˆ2,ˆ) = arg min • (parameter and ) For simplicity, let = + . Define ∈ N as the number of segments in which the function is monotone. We counted manually, although it could also be estimated by, for example, where is the complete elliptic integral of the first kind, and the term in brackets (+1) is added only if both 0 and −1 are simultaneously increasing or decreasing. Now we can simply invert the Jacobi elliptic function at the endpoints = 0 and = − 1 to obtain and . From equations (4.27) and (4.28), we have the following: which can be rewritten as if is even.
Here we used the fact that sn −1 = • arcsin, cn −1 = • arccos, and is the elliptic integral of the first kind. Finally,ˆ= 0 ,ˆ= • (parameters 0 and 0 ) From the previous steps, using all the recovered parameters, construct a discrete elastica segment that starts at the origin,ˆ. Then, (4.55) Figure 5 illustrates typical examples of the fairing by the discrete elasticae obtained by using the above algorithm.

Background and outline
Sangawara (Japanese pantiles) are the most common type of roof tiles in Japan and are thought to be unique to the country. Although the number of buildings with sangawara roofs are decreasing, the landscapes with the remaining constructions having sangawara roofs are consider by the community to be one of the most beautiful scenes and culturally Japanese. Traditionally, sangawara were handmade from local clay by placing a clay plate on a wooden mold, beating it with a board Figure 5: Typical examples of the fairing by the discrete elasticae. Black squares: input curve, green triangles: initial guess, red circles: output elastica.
called tataki, and stroking it with a board called nadeita ( Figure 6). In recent times, they are mass-produced by metal mold presses in limited areas. The mold shapes are thought to be based on the shape of the sangawara in the handmade era, but companies keep their designs a trade secret, so it is not clear. We consider that it is important to characterize aesthetically pleasing curves like sangawara with mathematical formulas to be used in architectural design. Because of this, and the fact that the process involves bending the clay plate, we thought that the shape of sangawara could possibly be approximated by elasticae. In Section 5.2 we explain how the handmade sangawara (simply referred to as pantiles) were collected, in Section 5.3 we obtain the keyline of each pantile, and in Section 5.4 we approximate those keylines to discrete elasticae.

Generation of 3D keyline data of pantiles
The pantiles we measured had been used in a house built around 1900 in Settsu City, Osaka Prefecture (Figure 7). From the characteristic shapes of the pantiles, they were likely used from the original construction or replaced before the revision of the urban building law in 1924 after the Great Kanto Earthquake. According to the owners, most of the tiles were blown away when the 2nd Muroto Typhoon hit in 1961, so they were collected and re-roofed. After that, only a few of the pantiles were replaced before the house was demolished in March 2017. In our survey before the demolition, we found that the roofs of this house were covered by four different sizes of pantiles that ranged from 240 to 280 mm in working width. Prior to dismantling the building, we preserved six rows of pantiles (A to F in Figure 8) that covered those four sizes. The pantiles varied in shape due to their handmade nature, so we preserved six rows instead of only four individual pantiles. We measured 37 pantiles of the C and F with a working width of 270 mm (the most commonly used on this house), excluding the eave pantiles (C01, F01). We used the NextEngine's Ultra HD 3D laser scanner and placed it and the pantiles as shown in Figure 9 for the measurements. The mesh data obtained from the 3D laser scanner was read by 3D Systems' RapidWorks 64 4.1.0 reverse modeling software. We used this software to synthesize and decimate its polygons finer than the scanner's measurement accuracy (0.3 mm) and then automatically heal incorrect data and fill holes in the meshes.
The two front edges are the 3D keylines of the pantile that could be observed when pitched. The lower edge could be generated by extracting the outer boundary curve of the mesh with RapidWorks. The upper edge could be generated by extracting the curve network from the mesh data with RapidWorks. Figure 10 shows an example of the two 3D keylines (upper and lower edges) generated. Each of these 3D keylines form an open polygon.

Generation and conversion of 2D keylines
Plane fitting of 3D keyline by principal component analysis In order to generate 2D keyline data from the 3D keylines obtained from the scanning data, we need plane fitting of the 3D data. This can be done by projecting the points in R 3 on the plane which minimizes the sum of squared distances from the points. As is well-known, such plane is constructed by applying the principal component analysis of the so-called the covariance matrix. More concretely, let 1 , 2 , . . . , ∈ R 3 be three-dimensional column vectors that represent spatial We consider the covariance matrix defined by which is real, symmetric positive semi-definite matrix, and let 1 , 2 , 3 be orthonormal eigenvectors that correspond to eigenvalues 1 , 2 , 3 ( 1 ≥ 2 ≥ 3 ≥ 0) of . Then it is known that the plane that includes¯and is parallel to 1 and 2 (i.e., with 3 as the normal vector) minimizes the sum of the squared distances from 1 , 2 , . . . , (see, for example, [1]). Therefore, plane fitting is obtained by projecting 1 , 2 , . . . , onto . Following this idea, the -coordinates of the point sequence consisting of all the vertices of the open polygon that constitutes one 3D keyline are converted into the -coordinate system, where the center of gravity of the point sequence is the origin and the orientations of 1 , 2 , and 3 are the , and axes, respectively. The orientations of 1 , 2 , and 3 are determined to be those of the , and axes shown in Figure 9, respectively. When converting the coordinate in the -coordinate system into the coordinate in the -coordinate system,

Equilateral open polygon approximation of 2D keyline
The 2D keylines obtained by projecting the 3D keylines onto the planes in the previous paragraph (simply referred to as the 2D keylines before approximation) can be approximated by equilateral open polygons using the following procedure.
First, let the end point on the negative side of the axis of the 2D keyline before approximation (simply referred to as the left end) be the starting point 1 after approximation, and let the point on the keyline whose distance from 1 is (= 1.0 mm) be 2 . Next, let the point on the same keyline where the distance from 2 is and does not return to the left end ( 1 ) be 3 . In the same way, repeat the operation to let the point whose distance from is and does not return to −1 be +1 (Figure 11). However, if there are multiple points that satisfy this condition, let the point on Figure 10: Example of 3D keylines generated (C07).
the keyline closer to the right end (the end point on the positive side of the axis) be +1 ( Figure  12). When the only point that has a distance from is the point that returns to −1 , is the stop point of the 2D keyline after approximation ( Figure 13). All the points from the start point 1

Coordinate conversion of 2D keylines after approximation by principal component analysis
When measuring the pantiles as shown in Figure 9, we place them one by one on the table by hand. Therefore, the position and rotation angle of each pantile is slightly different. To compare the curves of the keylines of the measured pantiles, we have to eliminate the effects of the positions and rotation angles and convert them into a coordinate system determined by only the keyline shape. This effect is already eliminated by conversion from the -coordinate system into the -coordinate system, where the center of gravity of the point sequence consisting of all the vertices of the 2D keyline is , the orientation of its first principal component is the axis, and the orientation of its second principal component is the axis. However, the number and position of the points that compose the point sequence have been changed due to the equilateral open polygon approximation. Therefore, the center of gravity of the point sequence consisting of all vertices of the approximated 2D keyline does not generally coincide with , nor does the orientation of its first principal component coincide with the axis, nor does the orientation of its second principal component coincide with the axis. Because the relationship between the shape of the keyline and the coordinate system has been lost, the coordinate system must be converted again into one determined solely by the shape of the keyline. Therefore, by performing the principal component analysis again on the point sequence consisting of all the vertices of the 2D keyline after approximation, we convert the coordinates into the -coordinate system with the center of gravity of the point sequence as and the first principal component as the axis. In the -coordinate system, let 1 , 2 , . . . , be the two-dimensional column vectors that represent the plane coordinates of the points of the 2D keyline after approximation. Then their center of gravitȳ is given by¯= 1 ∑︁
Let 1 , 2 be normalized eigenvectors that correspond to the eigenvalues 1 , 2 ( 1 ≥ 2 ) of covariance matrix defined by Here, the directions of 1 and 2 are determined so that the diagonal components of ( 1 2 ) are all positive. When converting the coordinate in the -coordinate system into the coordinate in the -coordinate system, is satisfied because ( 1 2 ) is an orthogonal matrix. In the following, the -coordinates obtained by this coordinate conversion will be redefined as -coordinates.

Estimation of inflection point of approximate curve of 2D keyline
The approximate curve of the 2D keyline of the pantile has one global inflection point and is clearly asymmetric across the inflection point. In general, it is difficult to formulate such a curve in terms of a single elastica. Therefore, we estimate the global inflection point of the approximate curve, then we divide the curve on both sides of the inflection point and approximate them with different elastic curves. We estimate the inflection point using a method inspired by the Ramer-Douglas-Peucker (RDP) algorithm [6,28], which was devised to simplify an open polygon with many vertices by thinning out the vertices (Figure 14). We call this the RDP method:  (4) When points + 1 and + 2 become adjacent vertices on the open polygon, the calculation is terminated and the inflection point is estimated to be + 2.
Because the RDP algorithm is intended to thin out the vertices, its calculation is terminated when an error falls within a certain range. However, the RDP method in this study is meant to estimate the inflection points. Therefore, we repeat the calculation only in the section where the inflection point is estimated to exist until it agrees with the original open polygon, and we estimate the converged point to be the inflection point.

Approximation of 2D keyline by discrete elastica
For the right (valley) sides of the inflection points of the 37 lower keylines (Figure 15) of the 37 pantiles, we tried to approximate the discrete elastica by the method already described. Among = ( 0 , 0 , ℎ, , , , ) shown in (4.5), ( 0 , 0 ) were set as the inflection point to allow the discrete elastica to pass through the inflection point after approximation, and the remaining five parameters were optimized. An example of approximation of a keyline of a pantile (C07 lower keyline) by discrete elastica is shown in Figure 16, the calculation results of the 37 approximated discrete elasticae are shown in Figure 17, and the calculation results of the parameters ℎ, , , , are shown in Table 1 As Table 1 shows, the 35 key lines (except for F16 and F24) show very little variation in the calculated results of the parameters, and the discrete elasticae are similar in both shape and rotation angle. The values are close to 0.3, and the shapes are close to sine curves. However, the values of , , and for F16 and F24 are extremely different from the others. The values exceed 10, and the shapes are close to arcs. Figure 18 shows the 37 translated discrete elasticae shown in Figure 17 with the inflection point ( 0 , 0 ) at the origin. The variations of the discrete elasticae are larger than that in Figure 17, suggesting that a certain number of variations and errors may be included in the estimations of the Figure 14: Example of inflection point estimation by RDP method (C07 lower keyline). In the first three curves we apply steps (1), (2) and (3), respectively. In the fourth and fifth curves, following step (3), we switch point 4 to point 5, and point 7 to point 11, respectively. Finally, in the augmented image at the bottom, following step (4), the inflection point is estimated to be point 11.
inflection points. In F16 and F24 in particular, the estimated inflection points are located closer to the right (valley) sides. Therefore, we can assume that the areas around the inflection points where the curvature is small are largely omitted, and the discrete elasticae are approximated to be close to the arcs. The effect the local unevenness of the keyline has on the estimation of the inflection point should be examined. In conclusion, we found that 35 of the 37 lower keylines of the handmade pantiles could be approximated by discrete elasticae with a very small variation on the right (valley)