Storage-efficient reconstruction framework for planar contours

Abstract A storage-efficient reconstruction framework for cartographic planar contours is developed. With a smaller number of control points, we aim to calculate the area and perimeter as well as to reconstruct a smooth curve. The input data forms an oriented contour, each control point of which consists of three values: the Cartesian coordinates (x, y) and tangent angle θ. Two types of interpolation methods are developed, one of which is based on an arc spline while the other one is on a cubic Hermite spline. The arc spline-based method reconstructs a G1 continuous curve, with which the exact area and perimeter can be calculated. The benefit of using the Hermite spline-based method is that it can achieve G2 continuity on most control points and can obtain the exact area, whereas the resulting perimeter is approximate. In a numerical experiment for analytically defined curves, more accurate computation of the area and perimeter was achieved with a smaller number of control points. In another experiment using a digital elevation model data, the reconstructed contours were smoother than those by a conventional method.


Introduction
A storage-efficient framework for handling cartographic planar contours is developed in this article. With a smaller number of control points for a target contour, we aim to calculate the area and perimeter of the contour as well as to draw a smooth curve. The framework can also be positioned as a contour reconstruction method.
Contour reconstruction is an important task in various engineering fields such as structural analysis, image processing, geographical analysis, and so on. A contour is a curve along which a function associated with a physical quantity has an equal value. In the usage of geography, examples of the quantity include elevation, depth, atmospheric pressure, temperature, air moisture, and so on. Although we hereafter focus on applications to cartographic systems, the same framework can also be applied to other fields.
In general, a contour is reconstructed through two interpolation procedures: we first find control points on a target curve followed by connecting the points with a curve. In most practical cases, the control points are not given directly and estimated using simulated or observed values on discrete points. The physical value of an arbitrary point is calculated by giving an interpolation function; in handling a planar contour, we set z = f (x, y) and consider a contour specified by z = c, where c is a target elevation. The most simple and commonly used function is a piecewise linear function for given elevations on a rectangular or triangular mesh, in which the resulting contour becomes a polygon and all control points can be calculated without an iterative calculation.
There are other popular functions such as bilinear and bicubic splines (Li, Zhu, and Gold 2004) to obtain a smooth curve. The resulting shape by bilinear spline is piecewise hyperbola. A disadvantage of this model is that the reconstructed contour breaks on the edges of a mesh, which is called an edge effect (Schneider 2005); the curve is G 0 continuous but not G 1 continuous. On the other hand, bicubic splines can reconstruct G 1 continuous curves but needs an iterative calculation to find control points.
The second concern in contour reconstruction is how to connect control points with a certain curve, which requires another interpolation scheme. The simplest method is to adopt a polygonal approximation, in which the ordered control points are connected with straight segments. In this method, however, the resulting contour tends to be too angular and occasionally unnatural, which becomes notable if the contour is magnified or the observation density is low.
To overcome this flaw, there would mainly be two approaches. The first one is to increase the number of control points, which can be achieved by subdividing a mesh or by using a continuous interpolation function.

OPEN ACCESS
In view of these needs, the authors have developed contour detection algorithms called the inner skin algorithm Shimakawa 2016, 2013). A feature of the algorithm is that each control point consists of three elements: the Cartesian coordinates (x,y) and tangent angle θ. The reconstructed curve forms an oriented contour, the area and perimeter of which can be calculated using Green's theorem (Stroud and Booth 2005). Goto and Shimakawa (2013) introduced a framework based on a first-order approximation, while Goto and Shimakawa (2016) managed to reduce the number of control points by a piecewise arc approximation. An oriented contour has another benefit by which the contour which is a part of a mountain or a basin can be judged; along with a visualization technique, a hill and punchbowl can be distinguished.
These works, however, do not pay attention to smoothing curves, and the reconstructed curves are G 0 continuous. Therefore, in this article, we develop a framework for reconstructing smooth contours based on arc and cubic Hermite splines. As introduced above, there are many existing researches on smoothing curves, but to the authors' best knowledge there is no research focusing on calculating geodetic quantities such as area and length. The arc spline-based method achieves G 1 continuity and exact area and perimeter can be calculated by subdividing each portion into two arcs. By contrast, the cubic Hermite spline-based method can reconstruct G 1 continuous curves and achieve G 2 continuity on some points. The exact area can be computed without subdividing each portion, but the exact perimeter is approximate. Since an accurate computation of the area and length is achieved with a small number of control points, the proposed framework is storage-efficient. In summary, the focuses of this article are threefold: • Reconstruction: draw a smooth contour; • Computation of area: calculate the area of a closed contour; • Computation of perimeter: calculate the length of a contour.

Storage format and intended use
The commonly used storage format for handling planar contours is the ordered set of two-dimensional Cartesian coordinates in the following manner: where n is the number of control points. To handle a closed contour, it is convenient if we set x n , y n = x 1 , y 1 , wherein the number of independent control points is n-1.
If we aim to draw a smooth contour or calculate feature quantities accurately, the control points are usually (1) x 1 , y 1 , x 2 , y 2 , ⋯ x n , y n For example, McCullagh (1981) handled rectangular and triangular meshes, Dyn, Levin, and Liu (1992) proposed a method by generating polyhedrons in three-dimensional space, and Gousie and Franklin (2005) adopted a spline function whose coefficients are determined using the thin-plate technique (Briggs 1974). Chai, Miyoshi, and Nakamae (1998) constructed a method to generate intermediate points by solving a partial derivative equation derived from an analogy of electric field. Furthermore, Pradhan and Ghose (2013) proposed a creative method to make a dense mesh by mimicking a spider of weaving a web.
This type of approach, however, requires recalculation whenever the target resolution is changed. Hence, a compromise approach is to calculate control points by constructing a dense mesh followed by eliminating intermediate control points that do not significantly affect the shape. This procedure is occasionally called contour compression or decimation, of which a pioneering method is called the Ramer-Douglas-Peucker (RDP) algorithm (Ramer 1972;Douglas and Peucker 1973). There are also other methods derived from the RDP algorithm (Ukasha, Saffor, and Hassan 2008;Ukasha 2011).
The second approach is to interpolate the ordered control points with a smooth curve, also known as contour smoothing. With respect to smoothing curves or surfaces for given control points, it has been a popular topic in the field of computer graphics. For example, Lawonnn et al. (2014) developed a smoothing method based on a reduction of geodesic curvature with an application example to surgical planning; Taubin (1995) introduced a method based on Gaussian smoothing that does not cause a shrinkage. The common situation in these researches is that the observed or simulated data are abundant but they include noises or errors.
By contrast, in cartographic systems, which are the main target of this article, the density of observed data is occasionally low due to technical and budget limitations. Thus, it is of use to develop a method for reconstructing a smooth curve with sparsely observed data. Adopting a spline-based framework can be popular and practical to fulfill this requirement. Li and Meek (2005) developed a method to achieve G 2 continuity for arc splines by inserting relaxed curves. Thibault and Gold (2000) developed a method called skeleton retraction by which minor ridges or valleys can be reconstructed.
In addition to visual smoothness, this article puts an emphasis on calculating the area and length of reconstructed curves with high accuracy. The aim is to apply the method to civil engineering and disaster prevention; for example the resulting area would be useful in predicting the region damaged by a flood, and length in estimating necessary resource or required time between two points or for a cycle.
taken at a small interval, whereas the data size tends to increase. Thus, it is useful to develop a framework by which the data size can be reduced.
The proposed storage format is simple; an index θ i is added to each control point x i , y i , the resulting entity of which is If the contour is closed, we set x n , y n , n = x 1 , y 1 , 1 . Here we denote the coordinate of the ith control point by x i = x i y i T . We write s i = x i+1 − x i and let this be referred to as a difference vector hereafter. Index θ i is the angle from the difference vector s i to the tangent vector on point x i in counter clockwise order. Figure 1 would help in comprehending the relations of these indices. Since the number of elements for each row increases 1.5 times compared with Equation (1), a methodology would be storage efficient if the same or better performance can be achieved with less than two-thirds of the control points. Let D and ∂D be the target closed contour and the oriented circuit along the boundary of the contour, respectively. In Figure 1 again, the area of D can be calculated by accumulating the area elements of triangle O-x i -x i+1 and the portion surrounded by segment x i -x i+1 and ∂D. We shall denote these by S 0 and S 1 hereafter. The former is equal to the value of the area by polygonal approximation, the specific value of which is: The sign is positive if x i+1 is located in the counter-clockwise side with reference to segment O-x i . The area S 1 differs from the interpolations methods, which will be examined in later sections.
The perimeter or the length of a part of ∂D can be calculated by accumulating the length elements of portion x i~xi+1 along ∂D. We denote the length of this portion by L 1 . This length clearly differs between the interpolation methods.
(2) x 1 , y 1 , 1 , x 2 , y 2 , 2 , ⋯ x n , y n , n  For reconstruction, we consider drawing a smooth contour, the primary essence of which is to determine the coordinates of the intermediate control points on portion x i~xi+1 along ∂D. To achieve scaling to any resolution, we develop a method for which the number of intermediate points can be varied arbitrarily.

Arc spline
In previous study by Goto and Shimakawa (2016), a method to partition ∂D based on roundness is proposed. The interpolation is performed by piecewise arc approximation, by which the exact area and length can be calculated. We can determine a unique arc by two control points x i , x i+1 , and the tangent angle θ i .

Calculation and reconstruction of a G 0 continuous arc
Let the center and radius of a target arc be denoted by c i and ρ i , respectively. If we put s i = x i+1 − x i and s i = |s i |, in view of Figure 2, 2 i sin i = s i follows. The area element S 1 is calculated by subtracting the area of triangle The length L 1 is calculated by If |θ i | ≃ 0, we use the following approximations because the above formulae become indeterminate: which are obtained by performing the Taylor expansion around θ i = 0.
We next draw a smooth arc by interpolation. In Figure 2, let y, c, and φ be an arbitrary interpolation point, tangent vector on endpoint x i , and angle from vector y − x i to m i in clockwise order, respectively. Angle φ follows 0 ≤ ≤ i .
The relations between the second and third terms are expressed as The values of cos i and sin i are then calculated in accordance with The angle φ i is obtained by calculating the arcsine of sin i .
If θ i . φ i < 0, we assume that the portion is a C-shaped curve and there is not an inflection point. By contrast, if θ i . φ i > 0, we assume that there is one inflection point on the portion for which an S-shaped curve is formed. We herein consider the two cases separately.

Subdivision of a C-shaped curve
To achieve G 1 continuity on the endpoints, we subdivide the portion into two arcs. Here we consider the case for θ i . φ i < 0, where the portion is a C-shaped curve. Figure 4 depicts a case of θ i > 0 and φ i < 0 with the relations of the two arcs and its pertinent indices, where we wrote θ = θ i and φ = -φ i . Let 1 , 2 and y be the radii of the arcs and intersection of the arcs, respectively. We determine ρ 1 and ρ 2 for which the tangents on y become both parallel to s i . Then the following relations hold considering the geometric properties: which can be written by then m i ∕m i represents the unit tangent vector, and can be represented by the following two forms: where R(θ) represents the rotation matrix for θ in counter-clockwise. We then obtain y as follows: Using this relation, the portion can be subdivided into any portions. One natural approach is to subdivide portion x i~xi+1 into N pieces, for which the length of each piece becomes smaller than or equal to a given target length s M . This length is given by considering the resolution of an output device. If we aim to determine N strictly, the following relation should be used: In practice, however, it is simple and would be sufficient to determine this using which we shall also use in other interpolation methods.

Assurance of G 1 continuity on the endpoints
A benefit of using an arc spline is that the feature quantities such as area and length can be computed without subdividing portions. However, the contour is not G 1 continuous on the endpoints if we use the method in the previous subsection. We hence consider assuring G 1 continuity by subdividing each portion into two arcs.
Let m i+1 and φ i be the tangent vector on x i+1 and angle from the difference vector s i = x i+1 − x i to m i+1 in counter clockwise order, respectively. With the help of Figure 3, the unit vector m i+1 ∕m i+1 on point x i+1 hold the following relations:  With respect to the tangent angles on points x i and y, the former angle is θ/2 = -θ i /2 in counter clockwise order, whereas the latter is φ/2 = φ i /2, which become the same angles in clockwise order as for θ i > 0 and φ i < 0.
In summary for the above two cases, if θ i . φ i < 0 holds, we assume that there is not an inflection point on the portion and subdivide it into two arcs. The intermediate point y is calculated using Equations (13) and (14), and the tangent angles on points x i and y are in clockwise order.

Subdivision of an S-shaped curve
The case for i ⋅ i > 0 is now examined. This case has an inflection point on the portion which we subdivide into two arcs. Letting be the inflection point, we adjust the radii of the arcs on which the tangents on are shared. Figure 5 illustrates the case of i > 0 and i > 0, where we expressed = i and = i simply. There is a degree of freedom in angle . Considering the geometric properties to vertical and horizontal directions, the following relations hold: which are expressed as Then the radii ρ 1 and ρ 2 are calculated according to Then the radii of the two arcs are obtained as and the intersection y is calculated as follows: where we replaced θ and φ with θ i and -φ i , respectively, on the last line.
To simplify this expression, we introduce the following two variables η and μ: We then obtain -φ i /2 = (1 -μ)η, and Equation (12) is rewritten as follows: If θ = θ i → 0 and φ = -φ i → 0 hold, η → 0 follows and, sin [(1 -μ)η]/ sin η becomes indeterminate. We thus derive an approximate form of Equation (14) by Taylor expansion around η = 0. The specific resulting form is With respect to the tangent angles, we first recall Figure 2, in which the tangent angle θ i on point x i is the half of the central angle of the approximated arc. Then in view of Figure 4 again, the tangent angle on point x i is θ/2 = θ i /2 and the angle on point y is φ/2 = -φ i /2 in clockwise order.
If θ i < 0 and φ i > 0, the positional relation of the portion Figure 4. If we put θ = -θ i and φ = φ i , then the geometric relation becomes the same as Equation (10) and the radii are expressed by Equation (11). On the other hand, the representation of intersection y is different from Equation (12), which results in By substituting θ and φ with -θ i and φ i respectively, we eventually obtain the same results as Equations (12) and (14). where we appended the prime symbol to to distinguish from that of Equation (18), noting the sign of the coefficients of s i is the opposite. If < i.e. i − i > 0 holds, point y is on the clockwise side with reference to segment x i ∼ x i+1 , and ′ must hold < ′ < . We thus set ′ to By substituting � = − , = − i , and = − i into Equation (24), we obtain the following simplified expression which results in the same as Equation (18). We can thus use Equation (21) to calculate y.
Regarding the tangent angles, the angles on points x i and y in counter clockwise order are which are eventually the same as Equations (22) and (23).
If > i.e. i − i < 0, point y is on the counter clockwise side of segment x i ∼ x i+1 , and we set ′ to the same as Equation (25). We can thus use Equation (21) to calculate y, and Equations (22) and (23) for the tangent angles.
In summary for the above four cases, if i ⋅ i > 0 holds, we assume there is one inflection point on the portion. The coordinate of the intermediate point y is calculated according to Equations (20) and (21), and the tangent angles on points x i and y are Equations (22) and (23), respectively.

Overall flow
The overall flow of the arc spline interpolation is confirmed. Table 1 presents a pseudo code of the proposed algorithm, the input of which is the set of Cartesian coordinates and tangent angle x i , y i , i (i = 1, ⋯ n). Algorithm 1 is designed to subdivide each portion into two arcs to make all endpoints G 1 continuous. The generic calls "Calc." are the abbreviation of "Calculate." Function 1 calculates the area and length elements for a single arc, followed subdividing the arc into multiple broken line segments. The area and length element can be calculated irrespective of the number of divisions N.
The coordinate of the inflection point y is then calculated as follows: With respect to angle δ, if < holds, point y must be on the counter-clockwise side with reference to segment x i ∼ x i+1 as shown in Figure 5. Considering the central angles of the two arcs, δ must hold < < , but there is a degree of freedom as long as δ holds this relation. One candidate is to take a certain average of θ and φ. From the results of several inspections, we propose to use the harmonic average given below: the reason for which is mentioned below.
If → 0, then → − ∕2 and → 0 follow in both Equations (13) and (20) but the signs of μ are the opposite of each other. Thus, the positions of y for Eqs. (14) and (21) are symmetric across segment x i ∼ x i+1 . This property is desirable for smaller tangent angles, and this is the reason for using the harmonic average.
In considering the tangent angles, we recall again that the tangent angle on a point is the half of the central angle. Then in view of Figure 5, the tangent angle on point x i is in clockwise order and the angle on point y is If > i.e. i − i > 0, the same geometric properties as Equations (17) and (18) hold, but point y is on the clockwise side of segment x i ∼ x i+1 , and must hold < < . If we set to the same as Equation (19), we obtain the same result as Equation (21).
In the case of i < 0 and i < 0, the positional relations of the relevant points and vectors are laterally inverted to those depicted in Figure 5. If we put = − i and = − i , the same geometric properties as Equation (17) hold, and the inflection point y is then calculated as To determine the coefficients of H, we set t = 0 and t = 1 to take into account the conditions in Equation (27): By putting, the above conditions are summarized as XHU = XI = X. Since this is followed for arbitrary X, matrix H must satisfy Accordingly, the interpolant p(t) is then represented as follows:

Cubic Hermite spline
Cubic Hermite spline, also referred to as Cubic Hermite interpolation, is widely used in computer graphics to draw a smooth curve for given control points. The curve is determined uniquely by giving two endpoints x i and x i+1 , and tangent vectors m i and m i+1 on the points.

Interpolant of a portion
An arbitrary point on the curve between portion x i ∼ x i+1 , denoted by p i , is represented in the following manner: where t is a continuous [0, 1] parameter. Functions h i (t) (i = 1, 2, 3, 4) are cubic functions of t referred to as basis function, the coefficients of which are determined by giving the following boundary conditions: For later convenience, we derive a matrix representation of Equation (26) here. First, the following symbols are introduced: and the suffix i of p i (t) will be dropped when no confusion is likely to arise. Then, p(t) and its derivative for t can be represented as follows: where N is the number of divisions introduced in Equation (8), on which the accuracy and computation time depend.

Determination of tangent vectors
In performing a cubic Hermite interpolation, it is necessary to determine the two tangent vectors m i and m i+1 in Equation (26). If the original data is given by Cartesian coordinates such as Equation (1), the tangent vectors must be determined only with x i (1 ≤ i ≤ n). Several methods have been proposed thus far, a simple and commonly used one of which is referred to as the cardinal spline (Chui 1988). We provide a tangent vector on x i by where c(0 ≤ c ≤ 1) is referred to as the shape or tension parameter. If c = 0, this method is called the Catmull-Rom spline (Catmull and Rom 1974). On the other hand, the focus of this research is to use the storage format proposed in Equation (2). This indicates that the direction of the tangent is fixed on every control point, whereas there is a degree of freedom in the length.
Here we express the tangent vector m i with two elements as follows: where m i and e i represent the length and unit tangent vector of m i , respectively. Once s i and θ i are given, e i is determined uniquely while there is a degree of freedom in determining m i (> 0). One straightforward variant is to take the length of m i in Equation (32) in the following manner: If c is close to 1, the tension is high, m i is small, and the portion x i ∼ x i+1 is approximately straight. By contrast, the portion becomes roundish if c is close to 0. The difference between using only Equation (32) and using Equations (33) and (34) lies in the direction on point x i , a demonstrative example of which will be presented in Section 5.

Smoothing the reconstructed curve
A benefit of using the cubic Hermite spline is that G 1 continuity is satisfied without subdividing a portion. On the other hand, there is a drawback in calculating the length element L 1 because Equation (31) is an approximate value and not an exact one.

Calculation of the area and length
We consider a method to calculate the area element S 1 surrounded by straight segment x i ∼ x i+1 and its interpolated curve. Noting the displacement of p(t) for an infinitesimal time dt is p � (t), S 1 is calculated in the following manner in reference to Equation (3): For simplicity, we translate x i to the origin and set x i = 0 and x i+1 = s i . Then, p(t) and its first and second derivatives are simplified into The integrand in Equation (28) is then represented by where we simply wrote u i = [Hu(t)] i and v i = [Hv(t)] i . After performing integrations of (u i v j -u j v i ) on the righthand side, we obtain Then, the area element S 1 can be calculated by The significance of this formula is that the resulting value is exact and not approximate. Once the tangent vectors are determined, the area can be calculated without subdividing the portion and does not include an approximation error. With respect to the length element L 1 , it is calculated based on Recalling Equation (29), the integrand is a square root of a quartic function of t, which is reduced to a combination of elliptic integrations and cannot be represented by elementary functions. We thus give up obtaining the exact value, and calculate an approximate one using a polygonal approximation below: Figure 6 will help in comprehending this concept, and this setting is performed for all (2 ≤ i ≤ n − 1). If the contour is closed, i.e. x 1 = x n follows, we set k 1 = k n = n−1 (1) + 1 (0) ∕2, whereas we set k 1 = 1 (0) and k n = n where This is derived by considering tiny variations of i (0) and i (1) in Equation (36). We repeat updating m until | m| < 1 is satisfied or the number of iterations reaches to a certain count, where 1 is a small positive constant.
Should one of the following conditions hold during the update, we discontinue updating and restore the initial values: where cond(J) and 2 are the condition number of matrix and a small positive constant. The first criterion is to avoid orienting the tangent vector in the opposite direction, while the second one is the condition for the existence of a solution of Equation (39). The third to fifth criteria arise from the property of cubic functions; a cubic function cannot have more than one inflection point. We thus exclude cases in which there must be at least two inflection points between two adjacent control points as shown in Figures 7 and 8. Figure 7 illustrates four cases with equal sign curvatures on the two endpoints. The upper-left case is an example with positive curvatures on points x i and x i+1 , in which e i is directed counter-clockwise with reference to s i . Figure 8 presents two cases with opposite signed curvatures, in which both To exploit another advantage, we consider smoothing the reconstructed curve for which G 2 continuity is partially achieved. Let p i (t) and i (t) (0 ≤ t ≤ 1) be the interpolant of the parametric curve and signed curvature of the ith (0 ≤ i ≤ n) portion, respectively. If i (1) = i+1 (0) holds, the curve is G 2 continuous on point i+1 . This can be achieved by adjusting the lengths of the two tangent vectors m i and m i+1 .
On a point determined by parameter t, the signed curvature i (t) is calculated based on where the sign is positive if the oriented curve moves in the counter-clockwise direction. Recalling Equation (29), the first and second derivatives of p i (t) on the two endpoints are calculated as follows: where we denoted m i = m i e i . By substituting these into Equation (35), we obtain the following curvatures on the two endpoints: These two equations can be reduced to a quartic equation for m i by eliminating m i+1 . The condition for the existence of positive solutions m i and m i+1 depends on the relations between e i × s i , s i × e i+1 , e i × e i+1 , i (0), and i (1), which is quite hard to discuss for general properties. We thus adopt a policy in which we try to smooth the curve only if there may be positive solutions m i and m i+1 that satisfy Equation (36). must be specified in addition to the same input as the arc spline interpolation. Algorithm 2 first calculates the current and target curvatures for each control point, followed by adjustment of the lengths of tangent vectors on the two endpoints. Function 2 interpolates a single portion with the tangent vectors on two endpoints specified. The function first calculates the area element, which can be obtained regardless of the number of divisions N. By contrast, the resulting length element is calculated by accumulating the lengths of the divided segments.

Numerical experiments
Two types of numerical experiments are performed. The first experiment compares the shapes of interpolated curves, while the second one inspects the relation between computation accuracy and storage efficiency. We use a PC equipped with an Intel® Core™ i5-5250U 1.60 GHz running Windows 7 64-bit. The programs are coded for MATLAB R2015b.

Comparison of the shapes
A simple closed contour featured by four control points is handled. In accordance with the storage format in Equation (2), we give the following data-set: with a constant tangent angle θ, noting the first and last points are the same to make the contour cyclic. The contour is P-shaped with both convex and concave curves. If one of the above cases arises during the update of m, we give up accomplishing G 2 continuity on the endpoints and restore the initial values of m i and m i+1 determined by Equation (37).

Overall flow
The overall algorithm of the cubic Hermite spline interpolation is summarized. Table 2 presents a pseudo code for the proposed framework. The tension parameter c Figure 7. cases with the same signed curvatures on the endpoints, in which there must be more than one inflection points. Figure 8. cases with opposite signed curvatures on the endpoints, in which there must be more than one inflection points. which is irrespective of . This yields that the direction of the arrows are different from those for cases (c)-(f).
In cases (c) and (d), four intermediate points are added between two adjacent control points by subdivision. In cases (e) and (f), the smoothing algorithm explained in Section 4.5. is applied in which a tension parameter c = 0.1 is used for an initial length.
In cases (c) and (e), all tangent vectors are directed to the inner side of the contour, while directed to the outer side in cases (d) and (f). The difference of the shapes between cases (c) and (d) is more remarkable than that of between cases (e) and (f); case (c) has a "slim" shape but case (d) is "plump. " The notable difference between cases (e) and (f) lies in the lengths of the tangent vectors. In cases where there is a steep curve and there must be an inflection point between two control points, the tangent vectors would be expanded to achieve G 2 continuity.
As the above results imply, various shapes of contours can be reconstructed by appending a tangent angle to a usual Cartesian coordinate.

Computation accuracy and storage efficiency
The relationship between the efficiency of the storage format and computation accuracy is inspected. We examine the following two closed curves: The areas of the ellipse and Cassini oval are ab and 2b 2 E(k), respectively, where k = (a/b) 2 and E(k) is the complete elliptic integral of the second kind. The perimeters of the curves are 4aE( , ;z] is the Gauss hypergeometric function. We now set (a, b) = (2, 1∕2) for the ellipse and (1, 1/0.999) for the Cassini oval, the profiles of which are depicted in Figure 10. We encode the curves in the form of Equation (2)  Marching Square (MS) (Lopes and Brodlie 2003;Maple 2003), with varying lattice intervals l = {5 ⋅ 10 −2 , 2 ⋅ 10 −2 , 1 ⋅ 10 −2 , 5 ⋅ 10 −3 , 2 ⋅ 10 −3 , 1 ⋅ 10 −4 , 5 ⋅ 10 −4 , 2 ⋅ 10 −4 } • RDP algorithm (Ramer 1972;Douglas and Peucker 1973), with lattice interval l = 2 . 10 −4 and varying tolerances d = {10 −2 , 10 −3 , 10 −4 , 10 −5 , 10 −6 , 10 −7 , 10 −8 }. MS is a contour detection algorithm for given elevations on a regular rectangular lattice; a contour is formed by connecting edges obtained for each patch number of control points increases 10-fold, the error level is decreased approximately 10 −4 -fold. The computation errors for the perimeter and area are thus estimated by O(1∕n 4 ). By contrast, the slopes of the MS and MS + RDP are approximately the same; the error is decreased approximately 10 −2 -fold if the number of control points is increased 10-fold. Thus, the computation errors are estimated by O(1∕n 2 ). These orders are the same for both the perimeter and area.
Furthermore, we measure the computation times. The time for each trial is measured using "tic" and "toc" functions, and the average time of 100 executions is adopted. Figure 13 presents the computation times for interpolating the curves including the calculation of the perimeter and area. The horizontal axis represents the number of control points n, while the vertical one gives the computation time in seconds. Each time approximately increases in proportion to n, which can be presumed in view of Tables 1 and 2.
The four values on the upper right represent the total number of interpolated points N including the original control points. For a fixed n, though N for the arc spline is two-to threefold compared with the Hermite spline, the arc spline is more lightweight than the cubic Hermite spline. This would be because the arc spline does not need a convergence calculation for adjusting the lengths of tangent vectors.
independently. The RDP algorithm is a contour compression technique by which the number of control points can be reduced; for two fixed control points, it is based on a straight line that passes the points. If the distances between intermediate points and the line are all under a given tolerance d, we remove the intermediate points. Since the RDP is designed to compress open curves, we split a closed contour at a medium control point and divide it into two series. Figures 11 and 12 compare the computation accuracies for four reconstructing methods: arc spline, cubic Hermite spline, MS, and the combination of MS and RDP (MS + RDP). The horizontal and vertical axes represent the logarithm of the computation error and the number of control points n, respectively. Figure 11 represents the errors for the Cassini oval, while Figure 12 presents those for the ellipse. The upper figures represent the errors for the perimeter, while the lower ones for the area. In the MS + RDP method, we first calculate intersections with MS for a fine lattice l = 2 . 10 −4 , followed by a reduction of the number of control points using RDP with varying tolerance d.
For a fixed error level, if a certain profile is located to the left of another one, it follows that the method is more efficient in storage than the other one. Moreover, if the negative slope of a profile is steeper than another one, the method is more beneficial in terms of space complexity than the other one.
In view of Figures 11 and 12, the arc splines significantly excel the other three methods in both storage efficiency and space complexity. With respect to the cubic Hermite spline, the errors for the Cassini oval are close to the MS + RDP method, while those for the ellipse are beneficial to some extent. The MS + RDP method excels the MS in most cases, which would because the control points are generated from the MS with the smallest lattice interval. In addition, the errors of the MS + RDP converge on a level close to the MS in n > 10 4 .
In terms of space complexity, the arc spline yields the best performance among the four methods. If the  rectangular lattice are recorded, and adopt a bicubic spline function (Li, Zhu, and Gold 2004) for an elevation function.
As a study area, we choose an isolated island called Yakushima, located in the southwest part of Japan, an overview of which is shown in Figure 14. There is the highest peak called Mt. Miyanoura, which is close to the center of the island. The peak is 1936(m) in altitude and located on 30 °20′9″N, 130°30′17″E. We use a DEM data, the lattice interval of which is 2.25 × 1.5 (s) in longitude and latitude directions. This corresponds to approximately 60.1 × 46.2 (m) around the peak. We detect and encode contours using a piecewise arc algorithm (Goto and Shimakawa 2016). Now the existing and proposed reconstruction methods are compared by investigating a peninsula on the northern part and a small valley located in northeast. These targets are pointed to by black arrows in Figure 14. The outermost thick contours represent coastlines z = 0 and the internal ones are those for every five meters. We extracted coastlines by assuming the elevation of the sea level is z = −1m, the reproducibility of which is thus not necessarily good. As long as the shapes of the reconstructed curves are compared, there seem few differences between the arc spline-and cubic Hermite spline-based methods. Contrary to the experiments in the previous section, which spline is used may not notably affect the visual shapes.
In view of Figure 15, the contours reconstructed by the MS are angular, and the contours break only on the edges of patches. Since a linear interpolant is used in the MS, the contour interval within a single edge is constant. Moreover, the peak of a local hill appears only on an intersection of four adjacent patches. These features are unprobable in real terrains and a flaw compared

Application to geographic information systems
Practical application examples to geographic information systems are presented. We use a digital elevation model (DEM) data in which the elevations on a regular    Figure 16 depicts the reproduced contours around the valley in the northeast. In this zone, a river flows from the lower left toward the upper right, which is not depicted in the figures. In the proposed methods, several closed chain lines are observed. These indicate that there is a concave inside the contour. Since there are double chain contours in the lower left part, the depths of the concaves are at least 5 m. This is a benefit of using oriented contours, and such features cannot be identified by the conventional methods.
The proposed methods would be useful for both personal and engineering purposes. In contrast to using contour image data, the proposed approach can detect contours for arbitrary elevations, as well as reconstruct curves for an arbitrary resolution. Since the input data is vector-based, the storage space can be reduced in mapping applications. In the use of civil engineering, the G 1 and G 2 continuities are an important feature in designing a walkway or cycle road. For another application, the volume of a pit or hill can be estimated more accurately, which is achieved by accumulating the areas of with the proposed methods. Another disadvantage is that ambiguity arises in connecting points on edges where the elevations of two opposite corners of a patch are greater or smaller than those of the remaining two elevations. We examine a hill on the upper right as an example. Focusing on the patch in the left result pointed to by a dashed arrow, there seems a saddle which is not observed in the middle and right results. This would be an adverse effect of the ambiguity. According to a real topographical map, no saddle exists around there and there is only a single hill.
The salient benefit of the proposed methods is that detailed computations and analyses are achievable even if only a coarse data-set is available. For example, because of the good reproducilibity of local convexities or concavities, small local peaks or pits can be located more precisely. Ridges and valleys would also be detected more easily. If we aim to find a route to climb up to the center hill from a coast, the path indicated by the dashed arrow in the right figure can be a candidate. It is not easy to find this from the left figure.  closed contours to vertical direction. Such usage is useful for example in predicting the area of water immersion caused by torrential rain or a landslide dam.

Conclusions
This research has developed a storage-efficient reconstruction framework for vector-based planar contours. The storage format is Cartesian coordinates with tangent angles. Two types of interpolation methods were developed; one is based on the arc spline, while the other one is on the cubic Hermite spline. In a numerical experiment for analytically defined curves, the arc spline-based method achieved better performance in computation accuracy than the cubic Hermite-based method. The shapes of reconstructed curves by the two methods would be notably different where there is a steep curve. By contrast, in another experiment using DEM data, there were little differences in the visual shapes. Which spline should be adopted depends on its objective; a benefit of using the arc spline is that an exact perimeter can be calculated, while with the cubic Hermite spline G 2 continuity is achieved in most control points. Although this research has examined using overland DEM data, the same framework is also applicable to undersea data. The calculation of a blackish-water region and estimation of total reserves are probable applications.