Contour tracing for geographical digital data

: Our purpose is to trace a contour in the form of a polygon. In this research, we use a bicubic spline function for interpolation of the elevation data on a grid covering the area of concern. We construct the polygon as a data consisting of ordered contour points on sides of the grid. The contour enters a cell at an entry point and goes out at an exit point on its sides. The polygon is formed connecting these points. A problem occurs as to which two points should be connected when a cell of the grid has more than three contour points on its sides. As for existing methods of the differential geometry such as discretization using tangential increments, it is difficult to predetermine a suitable step size to arrive at a next contour point correctly if several contour components wind closely to each other within a cell. As a solution, we take an algebraic approach exploiting a simple fact that a bicubic function is viewed as a univariate cubic function with a parameter. From this perspective, we identify the exit point examining the behavior of the real roots of the cubic equation for the contour in terms of the numerical order. Our method enables us to faithfully trace the contour of bicubic spline functions which provide smoother and better fitting curves than bilinear spline functions used by the other authors. Computation time is exhibited in the numerical experiment for an island in Japan.


PUBLIC INTEREST STATEMENT
This research paper is intended for researchers and practitioners who may be interested in making a contour map from geographical digital data. The data-set consists of altitudes at vertices of a rectangular grid. We interpolate these data by a bicubic spline function. We then calculate the locations of contour points on the sides of the grid. These points are connected to form a polygon which represents a contour of the altitude.

Introduction
We develop a method to trace a contour for the surface interpolated by a bicubic spline function in the form of a polygon. The principal features of this method are: • We use the elevation data given on a rectangular grid.
• We interpolate the elevation data by a bicubic spline function which provides fitting and smooth approximation with moderate computational load.
• We detect all the contour points on the sides of the grid and connect them to trace a contour.
• In connecting contour points on the sides of a cell, contour points can be linked correctly even if there are more than three points.
• We can obtain the same polygon no matter which point we choose on the contour as a starting point, different from iterative methods using discretization.
We store contours in the form of polygon data which is easy to handle when applying to various engineering needs such as calculation of the volume of a water reservoir.
The algorithm of tracing and drawing contour is sometimes referred to as "contour tracing" and "contouring." The pioneering work is Cottafava and Moli (1969). Their method is to compute contour points on a rectangular grid and connect them together to form a contour. Since they use the intermediate value theorem to detect a contour point, the number of contour points on a side of a cell is always presumed to be no more than one. McLain (1971) uses the weighted least square method to approximate the elevation function. This method cannot express the contour winding intricately since the direction of increment vector is presumed to change at an angle of zero or 45 degrees to the previous one. Lopes and Brodlie (1998) use bilinear splines to interpolate the elevation function. By incorporating one or two contour points in the interior of a cell, they enhance the preciseness of contour tracing for the bilinear hyperbola. Maple (2003) uses the method called marching squares which approximates the contour in a rectangular cell by 16 possible patterns.
Recently, in geographical modeling and analysis, the importance of application using discrete altitude data such as geographic information system is widely recognized. In particular, the Triangular Irregular Network (TIN) expressing a terrain surface in a set of triangles is frequently used. The TIN allows the arbitrarily scattered points. For terrain representation, they use only the key points which are less than those aligned regularly in Digital Elevation Model (DEM). Gold, Charters, and Ramsden (1977) propose the method of contour trace using TIN, and Dobkin, Wilks, Levy, and Thurston (1990) extend it to mappings between higher dimensional spaces. Saux, Thibaud, Li, and Kim (2004) extract topographic properties on scattered elevation data model by TIN. Goto and Shimakawa (2015) propose the iterative method of discretization using gradient vectors of an elevation function which is constructed by interpolation with a bicubic spline function.
The preceding works mentioned above do not faithfully represent contours winding intricately in a rectangular cell. In this case, McCullagh (1981) subdivides a grid cell into a series of subcells enabling contours to be laced smoothly through the cell. In our research, however, we presume a fixed rectangular grid and still maintain the accuracy. The outline of our process is illustrated as follows. The altitude data on a rectangular grid are given in Figure 1. To reproduce the surface, we interpolate the altitude data on the grid by bicubic spline in Figure 2. We then calculate the contour points which are intersections between the rectangular grid and the contour. We finally produce a polygon of the contour for an altitude as in Figure 3. Although the process is same as that of Lopes and Brodlie (1998), we adopt bicubic spline functions for interpolating elevation data. Taking full advantage of their accuracy, we achieve representation of complicated contours such as having multiple points on a side of a cell shown in Figure 4 connecting A to B, C to D, and E to F. Our key idea which makes it possible is to connect the calculated contour points using the order of real roots for cubic equations.

Elevation function
We construct an elevation function over an m × n rectangular grid with interval δ by interpolating data of altitudes with a bicubic spline function. Although we suppose that a cell of the grid is a square δ × δ for simplicity, the discussion below can be applied easily to a rectangle δ x × δ y . There are various types of spline polynomial functions with different degrees. However, the most commonly used are cubic splines because of their precision and efficiency in computation. In particular, cubic B-splines are preferred. A base function of one-dimensional B-spline is given by For bicubic spline functions, k is set to 4 in expression (1). As for the definition of right-hand side of (1), readers should refer to Piegl and Tiller (2012). The parameters ξ i , ξ i+1 , …, ξ i+k are called knots and used for controlling the smoothness of the base function. In our case, they coincide with consecutive nodes of the grid. If there are no identical knots in expression (1), then the base function is C 2 -class at each knot. We generate base functions of bicubic spline by taking a tensor product of two base functions of one-dimensional B-spline. Thus, the base functions of two-dimensional B-spline are M 4,i (x) M 4,j (y) (i = 1, …,m; j = 1, …, n). Put (x i , y j ) = ((i − 1)δ, (j − 1)δ) for the coordinates of (i, j)-th grid point, set ξ i+l = x i+l (l = 0, …, k) for M 4,i (x), and ξ j+l = y j+l (l = 0, …, k) for M 4,j (y). The form of M 4,i (x) is given by is illustrated in Figure 5. Taking a suitable linear sum, we obtain a bicubic spline function: (1) (2)

Roots of the cubic equation
The contour of an altitude h 0 is given by the equation h(x, y) = h 0 which is a cubic equation in one variable while fixing the other. We use some aspects of cubic equations to draw contours. Thus in this section, we consider a general cubic equation Putting we obtain Where Let the cubic and quartic roots of 1 with the positive imaginary part be denoted by ω and i, respectively. According to the Cardano's method, we put The roots of (6) are given as follows.
In the case of D < 0, the equation has one real root and two imaginary roots.
In the case of D = 0, we have U = V = −q/2 from expression (9). If q < 0, x 3 is a simple root and x 1 , x 2 are a multiple root as illustrated in Figure 7. If q > 0, x 1 is a simple root and x 2 , x 3 are a multiple root as illustrated in Figure 8.

Behavior of the roots
We here generally consider a parameterized cubic equation Each coefficient is a continuous function of t. Solving Equation (12) for s, we obtain the roots s 1 (t), s 2 (t), s 3 (t) corresponding to x 1 , x 2 , x 3 for (4), respectively. The counterparts of q and D are parameterized as q(t) and D(t). The behavior of the roots is classified as follows.
(i) In the case of D(t) < 0, one root is real and the other roots are imaginary.
(iii) In the case of D(t) = 0 for t = t d , the Equation (12) have a multiple root. We then have is a simple root and s 2 (t d ), s 3 (t d ) are confluent as a multiple root.

Proposed algorithm
We trace the contour points on the sides of a grid cell based on the behavior of the real roots of the cubic equations. To make a polygon, we propose an algorithm composed of the following four phases: (i) Interpolation with bicubic spline functions.
(ii) Detection of the grid cells with relevant contour points on their sides.
(iii) Computation of the coordinates for the contour points.
(iv) Selection of the next contour point on the sides of a grid cell.
More precisely, (i) we interpolate the altitudes for the grid points with a bicubic spline function. (ii) we count the number of contour points on each side of a grid cell by Sturm's method (Bourdon, 1853) and detect the grid cells which the contour traverses. (iii) we compute the coordinate values of contour points on the sides of the grid cells detected above using a numerical method. (iv) we determine the exit point corresponding to the entry point and append it to the polygon. In the case that the cell has more than three contour points on the sides, we select the exit point according to the behavior of the real roots of the cubic equations. We elaborate on the first and fourth phases below.

Interpolation
We describe below an algorithm for determining coefficients c ij in expression (3). We give matrix representations A and B of the system of base functions for one-dimensional B-spline in x and y, respectively, as follows.
where (x i , y j ) (i = 1, …, m; j = 1, …, n) represent locations of grid points. Let the altitude at a grid point (x i , y j ) be denoted by h ij . Since the values of the spline function coincide with the altitudes at the grid points, we have

Selection of the next contour point
We focus on the fourth phase which is the mainstay of the proposed algorithm. In the course of construction, the polygon of a contour enters a cell at a contour point on a side which we call an entry point. It goes out of the cell at one of the other contour points on the sides which we call an exit point for the entry point. Thus, if the number of the contour points on the sides is two, the exit point is automatically determined. Otherwise, we select the exit point using the following algorithm.

Coordinate transformation
To facilitate the algorithm, each time we enter a new (i, j)-cell, we transform the coordinates from (x, y) to (s, t) as specified in Table 1. Note that the range of the (i, j)-cell is As a result, the cell is rescaled as a unit square and the entry point is relocated on the bottom side. By substituting (s, t) into (x, y) in expression (3), the elevation function is formulated by All the coefficients in expression (26) are polynomials in t of degree at most 3. The range of the

The case D(t) < 0 for 0 ≤ t ≤ 1
For the exit point, we select a point with the smallest t among the remaining contour points on the sides ( Figure 10).

The case D(t) > 0 for 0 ≤ t ≤ 1
We select the exit point using the numerical order of the three real roots of the cubic equations. The exit point is determined as one with the smallest t among the contour points on the sides which have the same order as the entry point. To obtain the order of a root, no information is required on the other roots. We just use the following fact: For a monic cubic function we have Thus, we determine the order according to Table 2.
In the case of Figure 11, if the contour line enters the right-hand cell at the point A which order is 2, then it leaves the cell at the point D which order is also 2. On the other hand, if the contour line enters the right cell at the point B which order is 3, then it goes out of the cell at the point C which order is also 3. Note that the order of a root is determined over the range −∞ < s < +∞.

The case D(t d ) = 0 for some t d with 0 < t d < 1
The process of selecting the exit point varies according to the sign of D ′ (t d ) and q(t d ): (i) The case D ′ (t d ) < 0 and q(t d ) > 0 If the entry point is of order 1, among the other contour points, we select a point of order 1 with the smallest t < t d if any (Figure 12). If there is no such point, we select a point with the smallest t > t d (Figure 13).
If the entry point is of order 2(3), we select a point of order 2(3) with the smallest t < t d if any (Figure 14). If there is no such point, we select a point of order 3(2) with the largest t < t d (Figure 15).
(ii) The case D ′ (t d ) < 0 and q(t d ) < 0 If the entry point is of order 3, among the other contour points, we select a point of order 3 with the smallest t < t d if any (Figure 16). If there is no such point, we select a point with the smallest t > t d (Figure 17).
If the entry point is of order 1(2), we select a point of order 1(2) with the smallest t < t d if any ( Figure 18). If there is no such point, we select a point of order 2(1) with the largest t < t d (Figure 19).
(iii) The case D ′ (t d ) > 0 and q(t d ) > 0 Among the contour points other than the entry point, we select a point with the smallest t < t d if any (Figure 20). If there is no such point, we select a point of order 1 with the smallest t > t d (Figure 21).

Figure 20. The case D ′ (t d ) > 0 and q(t d ) > 0, enter as a unique real root (t < t d ), exit at (t < t d ).
(iv) The case D ′ (t d ) > 0 and q(t d ) < 0 Among the contour points other than the entry point, we select a point with the smallest t < t d if any (Figure 22). If there is no such point, we select a point of order 3 with the smallest t > t d (Figure 23).

The case a(t a ) = 0 for some t a with 0 < t a < 1
The process of selecting the exit point varies according to the sign of a ′ (t a )b(t a ) and D(t): (i) The case D(t) < 0 for 0 < t < 1 Among the contour points other than the entry point, we select a point with the smallest t < t a (Figure 24).
(ii) The case a � (t a )b(t a ) > 0 and D(t) > 0 for 0 < t < 1 If the entry point is of order 1(2), among the other contour points, we select a point of order 1(2) with the smallest t < t a if any ( Figure 25). If there is no such point, we select a point of order 2(3) with the smallest t > t a (Figure 26).
If the entry point is of order 3, we select a point of order 3 with the smallest t < t a (Figure 27).

Figure 23. The case D ′ (t d ) > 0 and q(t d ) < 0, enter as a unique real root (t < t d ), exit at order 3 (t > t d ).
Figure 24. The case D(t) < 0, enter as a unique real root (t < t a ), exit at (t < t a ).   (iii) The case a ′ (t a )b(t a ) < 0 and D(t) > 0 for 0 < t < 1 If the entry point is of order 2(3), among the other contour points, we select a point of order 2(3) with the smallest t < t a if any ( Figure 28). If there is no such point, we select a point of order 1(2) with the smallest t > t a (Figure 29).
If the entry point is of order 1, we select a point of order 1 with the smallest t < t a (Figure 30).

The case D(t d ) = 0 and a(t a ) = 0 for some t d , t a with 0 < t d , t a < 1
We perform the procedures of 3.2.4 and 3.2.5 according to the numerical order of t d , t a . In the case of more t d and t a , we iterate the procedures accordingly.  To illustrate the process, we consider a cell of the grid with solutions t a , t d for a(t a ) = 0 and D(t d ) = 0, respectively, assuming a ′ (t a )b(t a ) > 0 and D(t) > 0 (0 < t < t d ), D ′ (t d ) < 0 and q(t d ) < 0. In this case, we perform a procedure composed of 3.2.5-(ii) and 3.2.4-(ii). We start from an entry point of numerical order 1 shown in Figure 31. Firstly, we perform the process of 3.2.5-(ii) since t a < t d . No contour points of numerical order 1 are within the range of 0 < t < t a . Thus, the root curve containing the entry point is connected to the root curve of order 2 within the range of t a < t < t d . Secondly, we perform the process of 3.2.4-(ii) considering the intersection between the root curve and line t = t a as an entry point of order 2. No contour points of numerical order 2 are within the range of t a < t < t d . Thus, we select the contour point of numerical order 1 within the range of t a < t < t d as the exit point.

Selection of the next cell
In generating a polygon for a connected component of contour, we start from the entry point in the initial cell, and link the entry point of the current cell to the correct exit point and add it to the polygon. This process is repeated until returning to the initial point. In this process, we need to select the new cell for the entry point which is also the exit point in the current cell. In other words, we transfer

2.5-(ii) and 3.2.4-(ii).
to the next cell sharing the contour point with the current cell. If the contour point is not a vertex, these cells share a side. In this case, the next cell is clearly determined.
However, if the contour point is a vertex, the next cell is selected by the gradient of the elevation function which is perpendicular to the contour. Namely, if the gradient is in the direction of t-axis, we select the cell adjacent to the current cell in the direction of s-axis. If the gradient is in the direction of s-axis, we select the cell adjacent to the current cell in the direction of t-axis. If the gradient is neither in the direction of s-axis nor t-axis, the next cell is selected diagonal to the current cell. Each time the next cell is determined, we conduct the coordinate transformation in Section 3.2.1.

Case study of Yakushima island
In this section, we use the JMC50 m mesh data which is a proprietary product of Japan Map Center. These data are a remake of the altitude data on 50 m mesh of Geospatial Information Authority of Japan currently out of print. In this data, the altitudes of oceanic area are expressed by −999.9 m. We change it to −1 m for the purpose of interpolation following the precedent of Goto and Shimakawa (2015).
As a case study, we apply the algorithm described above to draw the contour map of Yakushima island which is located at the south part of Kagoshima prefecture. We use an altitude data on 498 × 598 grid. The grid interval is 50 m × 50 m. The altitude of each grid point is measured by meters. The highest point of Yakushima is at an altitude of 1,935 m.
We set m = 498, n = 598, and δ = 50 in expression (3) and obtain the spline function which is shown in Figure 32. As a result of the entire algorithm, the contour map is completed as shown in Figure 33 where the contours are drawn every 100 m from 0 to 1,900 m. To see the polygonal structure of the contours, we focus on area A in Figure 33, and enlarge it in Figure 34, where we add the contours for every 10 m from 0 to 100 m. In this scale, the piecewise linear structure of the polygons is recognizable. Although the discussion about the superiority of bicubic splines over bilinear splines is out of scope of our research, we compare these two splines by applying them to an area in Yakushima island in Figure 35. The contour of the elevation function by bicubic spline is colored black, while the one by bilinear spline is colored gray. The central cell of Figure 35 is magnified in Figure 36. The polygons constructed by our method are delineated by black solid lines, while the one by the method of Lopes and Brodlie (1998) is drawn by gray dashed lines in Figure 37 and Figure 38. Both methods faithfully trace the contour of spline functions, respectively. As a consequence, the advantage of bicubic spline interpolation is reflected on our polygon.

Figure 33. Contour map of Yakushima island.
We measure the time of computation for constructing each polygon as shown in Figure 39 where the level axis represents the number of points in a polygon and the vertical axis the computational time. The maximum computational time is measured 242.39 s for a polygon of 500 m consisting of 5,000 points. To compare the computational time between a cell with two contour points on its sides and a cell with more than two, let n 1 be the number of the former and n 2 be the number of the latter. The regression of computational time t on n 1 and n 2 is t = 0.046 n 1 + 0.108 n 2 − 10.68 with R 2 = 0.993. Although our algorithm of selecting the exit points seems complicated, the computational time is not so affected by the number of the contour points on the sides of a cell: The ratio of the regression coefficients for n 1 and n 2 above is about 1:2, while the ratios of n 1 and n 2 is around 95:5.