Contour detection with bicubic spline surfaces

ABSTRACT

Contour detection has a rich history in multiple fields such as geography, engineering, and earth science. The predominant approach is based on piecewise planar tessellation and now being challenged concerning the extraction of contour objects for non-linear elevation functions, particularly with respect to bicubic spline functions. A storage-efficient method was developed in previous research, but the detection of the complete set of contour objects is yet to be realized. Although intractable, theoretical underpinnings pertinent to curvature resulted in an approach to realize the complete detection of objects. Given a digital elevation model dataset, in this study, a bicubic spline surface function was first determined. Thereafter, candidate initial points on the edges across the region of interest were identified, and the recursive disaggregation of rectangles was repeated if the non-existence of a solution could not be assured. A developed tracking method was then applied. During advancement, other initial points on the same contour curve were identified and eliminated to circumvent duplicate detection. The completeness of the outlets provides analytical tools for elevation and other geographical assessments. Demonstrative experiments included the development of a three-dimensional contour-based network and slope assessments. The latter application transforms the slope analysis type from raster-based to vector-based.

Highlights

• Detection of a complete set of contour objects amenable to bicubic spline surfaces.

• Small closure inside a single patch is detectable if size exceeds the standard.

• Curvature & tolerances central to step length adjustment and tangent angle determination.

• Redundant initial points are identified and eliminated during the tracking process.

• Various potential applications in addition to geographical elevations.

KEYWORDS:

• Contour detection
• bicubic spline surface
• curvature
• tolerance
• patch disaggregation

1. Introduction

Generating a cartograph from the observed elevations is one among a wide range of applications of contour extraction. Along with a surface model in a Cartesian planar space, it is denoted by $z=f\left(x,y\right)$, and the elevation is calculated across the region of interest. Conversely, given an elevation of $z=\mathrm{h}$, the loci $\left(x,y\right)$ that yield this height, i.e. $f\left(x,y\right)=\mathrm{h},$ comprise curves referred to as contours (or isolines in several contexts). Contour identification and boundary detection remain fundamental components in various scientific and engineering fields such as topography (Oka, Garg, and Varghese 2012; Xu, Chen, and Yao 2016), meteorology (Matuschek and Matzarakis 2011; Malamos and Koutsoyiannis 2016), hydrology (O'Loughlin, Short, and Dawes 1989; Liu et al. 2012), disaster prevention and reduction (Parsian et al. 2021; Tadini et al. 2022), urban planning/security management (Qiao et al. 2019; Garfias Royo, Parikh, and Belur 2020), resource exploration (Haijun and Yongping 2010; Rui et al. 2011), biometrics (Guan et al. 2006; Li et al. 2016), image processing (Gong et al. 2018; Henry and Aubert 2019), industrial design (Andújar et al. 2005; Victoria, Martí, and Querin 2009), and electromagnetics (Hanyan et al. 2016; Yang and Dong 2020). Although the process for obtaining the elevation $z$ from location $\left(x,y\right)$ is straightforward, the opposite process of identifying loci $\left(x,y\right)$ from the elevation $\mathrm{h}\left(=z\right)$ is cumbersome. The resultant comprises a set of curved segments accompanied by a compromise in practice.

The predominant approach is a piecewise linear approximation: given observed or simulated values on a square or rectangular lattice, the entire surface comprises a set of triangles, in addition to continuity on all edges. Due to piecewise planarity, each contour entity is comprised of successive straight segments, i.e. a polygon. There are two major approaches to contour manipulation. The first type is based on tracking (Lihua, Ai-jun, and Lu-ming 2008; Rui, Song, and Ju 2010; Jang and Bang 2017): identifying an initial point $\left(x,y\right)$ that follows $f\left(x,y\right)=\mathrm{h}$, and the following point is determined to link these two points. This is repeated until the control point returns to the initial location or encounters a boundary, and a contour entity is composed. The second is referred to as the marching squares algorithm (Maple 2003; Neto, Santos, and Vidal 2016). By identifying intercepts on the borders of all square or rectangular patches, each patch is associated with a delineation pattern with reference to a lookup table. This is significant in that the associations are independent of each other, and therefore, the framework allows for parallel processing (Martínez-Frutos and Herrero-Pérez 2017; Tan et al. 2017; Zhou and Li 2021).

In addition to piecewise-linear tessellations, non-linear surface models are promising for achieving smoothness on both the original surface and contour curves. Bilinear spline interpolation (Schneider 2005; Shangguan et al. 2013) involves a simple non-linear model in which each contour segment is a portion of a hyperbola in a single patch. Although the segment in a patch is geometrically a curve, the first derivatives $f_x$ and $f_y$ are not continuous on the borders of the mesh. This entails the discontinuity of the tangent vector across borders. In particular, the interpolation scheme guarantees $G_0$ continuity (Song et al. 2016) across the entire region of interest, whereas it does not guarantee $G_1$ continuity on borders. Furthermore, the number of contour portions within a patch is one, if any. This translates into the limited reproductivity of the models.

Bicubic spline interpolation is an alternative method for enhancing smoothness. With 16 coefficients for a single patch, a higher degree of smoothness is achieved, in addition to richer reproductivity. The continuities of the three derivatives $f_x$, $f_y$, and $f_{xy}$ are guaranteed on borders, for which the tangent vector and its associated quantities are continuous across the entire region. Along these lines, interpolation is particularly common for surface delineation (Popowicz and Smolka 2014; Vannier et al. 2019). Aiming at further smoothness and flexibility, higher-order spline interpolations such as biquartic splines were addressed in Mino, Szabó, and Török (2016), Karim et al. (2021); however, the scope of this study was limited to the cubic order considering the main objective.

In contrast to surface delineation, contour detection for higher-order interpolation functions is challenging. The core objectives are threefold: (1) identify points that follow $f\left(x,y\right)=\mathrm{h}$, (2) identify pairs of two adjacent points, and (3) extract all contour portions. There are challenges to contour tracing for bicubic spline surfaces, as addressed in previous studies. In particular, the following point was determined using a first-order derivative of function $f$ (Goto and Shimakawa 2013) and three curvature-based quantities (Goto and Shimakawa 2016). Both involved the adjustment of the step length in accordance with the steepness of the curve around the current control point, which contributed to a decrease in the number of control points. In particular, the latter achieved significant storage efficiency. The frameworks constructed in Goto and Shimakawa (2017) were able to reproduce smoother curves, thus achieving $G_1$ and $G_2$ continuities.

Although the abovementioned studies demonstrated efficient computations, high storage efficiencies, and the reproduction of smooth curves, the identification of sufficient initial points for detecting all contour objects was not considered. This may cause several objects to be undetected in opposition to the geometrical presence. In particular, the contour outlets generated from the above frameworks may not be amenable to the original surface. Although it is intractable to detect a complete set of objects, the concept of tolerance is critical to solving this problem. The presence of small closures, even within a single patch, such as bulges and depressions, is detectable if the sizes are greater than a given tolerance size.

Furthermore, this study aimed to achieve an efficient elimination of initial points that belong to single-contour objects. This is in association with objectives (2) and (3). A single contour curve passes over multiple points, and multiple points may be detected as potential initial points. Once the locus of a single contour is determined, other initial points on the same curve are eliminated to avoid duplicate detection. In a previous study (Goto and Shimakawa 2016), duplication was checked by two feature quantities: the centroid and area. This process is secure; however, it is time-consuming because it can only be detected after the entire shape is determined. To efficiently eliminate redundant points, a method for gleaning neighboring points during the tracking process was developed in this study. A three-dimensional arc-shaped region regulated by piecewise arc approximation is central to achieving this objective.

Using the proposed framework, a complete contour map that is amenable to the original terrain surface is obtained. The detected control points comprising contour portions are storage-efficient, as they are obtained considering the steepness of the curves. Furthermore, in addition to an external module for detecting the steepest descending/ascending curves, in Goto (2021), a three-dimensional consistent contour-based network (CBN) (Maunder 1999; Rulli 2010) was generated. The significance of this framework is discussed in a further section.

2. Surface model manipulation

Contour detection across a region of interest translates into identifying loci $\left(x,y\right)$ that satisfy the following: $f\left(x,y\right)=\mathrm{h},$where $f\left(x,y\right)$ and $\mathrm{h}$ are the elevation function and target height. The following symbols associated with $f\left(\equiv f\left(x,y\right)\right)$ are defined as follows:

• $\mathrm{\nabla }f\equiv [f_x{f}_y{]}^{\mathrm{T}}$: gradient of $f$,

• $\mathit{n}\equiv \left[f_x{f}_y{]}^{\mathrm{T}}/\right|\mathrm{\nabla }f|$: unit normal vector,

• $\mathit{t}\equiv \left[f_y-f_x{]}^{\mathrm{T}}/\right|\mathrm{\nabla }f|$: unit tangent vector rotated $\mathit{n}$ in the counterclockwise direction,

• $H\equiv \mathrm{\nabla }{\mathrm{\nabla }}^{\mathrm{T}}f=\left(\begin{array}{cc}{f}_{xx}& f_{xy}\\ f_{xy}& f_{yy}\end{array}\right)$: Hessian of $f$.

2.1. Elevation functions

Among the two major variants for composing the elevation function as follows: (1) $z=f\left(x,y\right).$(1) there are tessellations comprised of triangular or rectangular patches. With respect to piecewise linear tessellations, whether triangular or rectangular, the elevation function for each patch has the following form: $f=a_0+a_{1}x+a_{2}y,$where the resultant local locus is a straight segment. Although useful in composition, the entire outlet comprises a set of broken line segments that fail to attain smoothness. In addition to a bilinear spline surface model expressed as follows: $f=a_0+a_{1}x+a_{2}y+a_{3}xy,$the locus in a single patch is a hyperbola. Although this may attain smoothness, the level of geometrical continuity remains at $G_0$, i.e. the tangent vector varies discontinuously across the borders of the patches.

The primary approach accompanying higher-level continuity relies on a bicubic spline function: (2) $f=\sum _{j=0}^3\sum _{i=0}^3{a}_{ij}{x}^i{y}^j.$(2) Given four elevations at the four corners of each rectangular patch, the continuities of $f$, $f_x$, $f_y$, and $f_{xy}$ are constrained on four borders, for which the 16 coefficients in (2) are determined uniquely. The entire surface achieves $G_1$ continuity across the region of interest, for which the visualized outlets of $z$ are smooth. Consequently, loci detection is intractable, as the bivariate polynomial with 16 coefficients renders it challenging to detect complete sets of loci.

When an initial point on or close to the exact locus is given, the frameworks developed in Goto and Shimakawa (2013, 2016) can be used to determine the step length by considering the steepness or roundness of the locus, thus allowing for the following control point to be determined. In the latter study, three curvature-based quantities were introduced for the identification of control points in a storage-efficient manner. However, these frameworks were designed to be applicable to general elevation functions in the form of (1) and operate while $f_x$ and $f_y$ are continuous. In exchange for usability, the detection of initial points relies on methods other than the bicubic spline function (2), for which the reproduced curves may not comprise a complete set of contours. In particular, certain entities may remain undetected, which is mathematically unavoidable.

2.2. Enhancing reproductivity

In cartographic contexts, this limitation may hinder the use of higher-order spline functions for the generation of contour- or CBN-based topographical maps. Although seemingly intractable for general elevation functions, the proposed solution allows for the detection of complete sets of contours and several restrictions. The first limits the elevation function to the bicubic-based type, i.e. Eq. (2(2) $f=\sum _{j=0}^3\sum _{i=0}^3{a}_{ij}{x}^i{y}^j.$(2) ), thus providing analytical higher-order derivatives to assess the steepness of the curves. The second restriction is to introduce a minimum standard with respect to the detectable scale, which is referred to as tolerance. Starting with the patch size of the original data, the potential existence of a closure within the patch is inspected. If the existence of a closure cannot be assured, it is subdivided for further inspection. This is repeated until an intersection of a contour is detected or the resultant disaggregation falls below the tolerance level. Depressions or bulges with sizes smaller than the tolerance level may not be detected and are therefore accepted. A smaller tolerance enhances the reproductivity and the computation time increases.

2.3. Curvature assessment

Contour curves break only on edges in piecewise linear tessellations, whereas they curve everywhere with non-linear elevation functions. Thus, step-length adaptation is critical to contour tracking. The method developed in Goto and Shimakawa (2016) mainly aims at storage efficiency, in which three curvature-associated quantities are computed with partial derivatives $f_x$ and $f_y$ at two successive points. These quantities enhance storage efficiency; however, they do not address the locus between the control points. Thus, the tracing method does not provide an approach to eliminate redundant initial points on the same trajectory. Accordingly, in this study, a local curvature vector was derived by utilizing analytical first- and second-order derivatives of $f$. This was achieved by limiting the elevation function to the bicubic model (2).

Curvature is defined as the change in the angle of the tangent vector per infinitesimal displacement. Many variants exist with respect to the trajectory and reference plane, and over dozen types of complete curvatures were summarized in Shary, Sharaya, and Mitusov (2002). Among others, the most popular quantities are profile and planform curvatures (Evans 1972; Krcho and Haverlík 1973), which are defined as follows: ${\kappa }_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{f}}=-{\displaystyle \frac{1}{{\left(1+{\left|\mathrm{\nabla }f\right|}^2\right)}^{3/2}}{\mathit{n}}^{\mathrm{T}}H\mathit{n},{\kappa }_{\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{n}}=-{\displaystyle \frac{1}{\left|\mathrm{\nabla }f\right|}{\mathit{t}}^{\mathrm{T}}H\mathit{t}.}}$The latter is referred to as plan curvature and is closely related to contour tracking in this study. Although curvatures are generally defined as scalar quantities, a vector-based curvature for direct application to trajectory tracking was derived in this study.

Let $\mathit{x}$ and $\delta l$ be the current control point and the step length, respectively. A curvature vector along the unit tangent vector $\mathit{t}\left(\mathit{x}\right)$ is then defined as follows: $\mathit{\kappa }=\underset{\delta l\to 0}{lim}{\displaystyle \frac{\mathit{t}\left(\mathit{x}+\delta l\cdot \mathit{t}\right)-\mathit{t}\left(\mathit{x}\right)}{\delta l}.}$The numerator of the right-hand side is approximated to the first-order of $\delta l$ by the following: $\mathit{t}\left(\mathit{x}+\delta l\cdot \mathit{t}\right)-\mathit{t}\left(\mathit{x}\right)\simeq \delta l\cdot J\left(\mathit{x}\right)t\left(\mathit{x}\right),$where $J\left(\mathit{x}\right)$ is a Jacobian matrix with respect to $\mathit{t}\left(\mathit{x}\right)$. The specific representation of $J\left(\mathit{x}\right)$ is as follows: $J\left(\mathit{x}\right)={\left[\left(\begin{array}{c}\mathrm{\partial }/\mathrm{\partial }x\\ \mathrm{\partial }/\mathrm{\partial }y\end{array}\right)\left({\displaystyle \frac{f_y}{\left|\mathrm{\nabla }f\right|}-{\displaystyle \frac{f_x}{\left|\mathrm{\nabla }f\right|}}}\right)\right]}^{\mathrm{T}}={\displaystyle \frac{1}{{\left|\mathrm{\nabla }f\right|}^3}\left(\begin{array}{cc}{f}_x\left(f_x{f}_{xy}-f_y{f}_{xx}\right)& f_x\left(f_x{f}_{yy}-f_y{f}_{xy}\right)\\ f_y\left(f_x{f}_{xy}-f_y{f}_{xx}\right)& f_y\left(f_x{f}_{yy}-f_y{f}_{xy}\right)\end{array}\right).}$Based on the above, the curvature vector along the tangent vector can be derived as followed: $\begin{array}{rl}\mathit{\kappa }& =J\left(\mathit{x}\right)\mathit{t}\left(\mathit{x}\right)\\ & ={\displaystyle \frac{1}{{\left|\mathrm{\nabla }f\right|}^4}\left(\begin{array}{cc}{f}_x\left(f_x{f}_{xy}-f_y{f}_{xx}\right)& f_x\left(f_x{f}_{yy}-f_y{f}_{xy}\right)\\ f_y\left(f_x{f}_{xy}-f_y{f}_{xx}\right)& f_y\left(f_x{f}_{yy}-f_y{f}_{xy}\right)\end{array}\right)\left(\begin{array}{c}{f}_y\\ -f_x\end{array}\right)=-{\displaystyle \frac{1}{{\left|\mathrm{\nabla }f\right|}^4}\left(\begin{array}{c}{f}_x\left(f_y^2{f}_{xx}-2f_x{f}_y{f}_{xy}+f_x^2{f}_{yy}\right)\\ f_y\left(f_y^2{f}_{xx}-2f_x{f}_y{f}_{xy}+f_x^2{f}_{yy}\right)\end{array}\right)}}\\ & ={\displaystyle \frac{-1}{\left|\mathrm{\nabla }f\right|}{\displaystyle \frac{\left(f_y-f_x\right)}{\left|\mathrm{\nabla }f\right|}\left(\begin{array}{cc}{f}_{xx}& f_{xy}\\ f_{xy}& f_{yy}\end{array}\right)\left[{\displaystyle \frac{1}{\left|\mathrm{\nabla }f\right|}\left(\begin{array}{c}{f}_y\\ -f_x\end{array}\right)}\right]\left[{\displaystyle \frac{1}{\left|\mathrm{\nabla }f\right|}\left(\begin{array}{c}{f}_x\\ f_y\end{array}\right)}\right]=-{\displaystyle \frac{1}{\left|\mathrm{\nabla }f\right|}{\mathit{t}}^{\mathrm{T}}H\mathit{t}\cdot \mathit{n}={\kappa }_{\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{n}}n,}}}\end{array}$where the final equality highlights a direct relationship between the vectorized curvature and scalar planform curvature.

3. Contour tracking

3.1. Advancement of control point

Given an initial point, advancement along the contour continues until the current control point returns to the initial point or encounters a boundary. Each advancement comprises determinations of the step length and tangent angle denoted by $\delta l$ and $\theta$, respectively. In contrast to piecewise linear surface models, the following control point is determined irrespective of the mesh borders. As depicted in Figure 1, considering two orthonormal vectors, namely, normal and tangent vectors, the positional relationship is represented by $\delta l$ and $\theta$. This representation is similar to a polar coordinate expression.

Figure 1. Positional relationships of the current and following control points, in addition to variables for advancement.

The step length $\delta l$ is first adapted to confine the deviation from the exact (yet unknown) curve into a standard referred to as the tolerance, as denoted by Tolx. Considering the advancement from the current control point to the following point, let $\delta \mathit{x}$ and $\theta$ denote the shift from the current control point and the tangent angle from $\mathit{t}$ (unit tangent vector) to $\delta \mathit{x}$, respectively. The shift vector is then expressed as follows: (3) $\delta \mathit{x}=\left(\cos \theta \cdot \mathit{t}+\sin \theta \cdot \mathit{n}\right)\delta l.$(3) Denoting the current point by $\mathit{x}$, the elevation at point $\mathit{x}+\delta \mathit{x}$ is estimated using the Taylor expansion with respect to $\delta l$ as follows: $f\left(\mathit{x}+\delta \mathit{x}\right)=f\left(\mathit{x}\right)+\sin \theta \left|\mathrm{\nabla }f\left(\mathit{x}\right)\right|\delta l+\mathrm{O}\left(\delta l^2\right),$where $\theta$ denotes the initial tangent angle. When $\delta l$ is determined and fixed, $\theta$ is then adjusted to contain $\left|f\left(\mathit{x}+\delta \mathit{x}\right)\right|$ within the tolerance. Angle $\theta$ is updated based on the Taylor expansion of $f\left(\mathit{x}+\delta \mathit{x}\right)$ as a function of $\theta$ as follows: $f\left(\theta +\delta \theta \right)=f\left(\theta \right)+{\displaystyle \frac{df\left(\theta \right)}{d\theta }\delta \theta +\mathrm{O}\left(\delta {\theta }^2\right)\simeq f\left(\mathit{x}+\delta \mathit{x}\right)+\left|\mathrm{\nabla }f\left(\mathit{x}+\delta \mathit{x}\right)\right|\delta l\cdot \cos \theta \delta \theta .}$To update $f\left(\theta +\delta \theta \right)$ closer to 0, $\delta \theta$ should be set as follows: (4) $\delta \theta =-{\displaystyle \frac{f\left(\mathit{x}+\delta \mathit{x}\right)}{\cos \theta \left|\mathrm{\nabla }f\left(\mathit{x}+\delta \mathit{x}\right)\right|\delta l}.}$(4) This process is repeated until $\left|\delta \theta \right|$ is contained in the standard. The difference of $\delta \mathit{x}$ impacted by the above is inspected until it is smaller than Tolx, which results in the following: $\left|\mathrm{\Delta }\left(\delta \mathit{x}\right)\right|=\delta l\left|\delta \theta \right|={\displaystyle \frac{\left|f\left(\mathit{x}+\delta \mathit{x}\right)\right|}{\cos \theta \cdot \left|\mathrm{\nabla }f\left(\mathit{x}+\delta \mathit{x}\right)\right|}<\mathrm{T}\mathrm{o}\mathrm{l}\mathrm{x}.}$The elevation tolerance is then formulated as follows: (5) $\left|f\left(\mathit{x}+\delta \mathit{x}\right)\right|<\mathrm{T}\mathrm{o}\mathrm{l}\mathrm{x}\cdot \cos \theta \cdot \left|\mathrm{\nabla }f\left(\mathit{x}+\delta \mathit{x}\right)\right|.$(5) Although the right-hand side enforces $\cos \theta >0$, i.e. $\left|\theta \right|<\pi /2$, a safer constraint is applied as follows: (6) $\left|\theta \right|\le \pi /4.$(6) If the sign of $f\left(\mathit{x}\right)$ changes during the iteration, there is a solution between two successive updates. Let $\delta {\mathit{x}}_m$ be the shift vector at the $m$th step. There is a solution if the following is true: (7) $f\left(\mathit{x}+\delta {\mathit{x}}_{m-1}\right)\cdot f\left(\mathit{x}+\delta {\mathit{x}}_m\right)\le 0.$(7) In addition, there are constraints derived by (5) if the following is true: (8) $\begin{array}{rl}& \left|f\left(\mathit{x}+\delta {\mathit{x}}_{m-1}\right)\right|<\mathrm{T}\mathrm{o}\mathrm{l}\mathrm{x}\cdot \mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{m-1}\cdot \left|\mathrm{\nabla }f\left(\mathit{x}+\delta {\mathit{x}}_{m-1}\right)\right|,\\ & \left|f\left(\mathit{x}+\delta {\mathit{x}}_m\right)\right|<\mathrm{T}\mathrm{o}\mathrm{l}\mathrm{x}\cdot \mathrm{c}\mathrm{o}\mathrm{s}{\theta }_m\cdot \left|\mathrm{\nabla }f\left(\mathit{x}+\delta {\mathit{x}}_m\right)\right|.\end{array}$(8) The update of $\theta$ is then discontinued to proceed to the following control point.

3.2. Restriction and adaptation of the step length

Key restrictions for the step length $\delta l$ are associated with the deviation of the exact curve. In addition to the piecewise arc approximation, the tangent angle is related to the step length $\delta l$ and plan curvature $\kappa$ as follows: (9) $\mathrm{s}\mathrm{i}\mathrm{n}\theta ={\displaystyle \frac{\kappa }{2}\delta l.}$(9) Combined with Constraint (6), $\delta l$ should be restricted as follows: (10) $\delta l\le \sqrt{2}/\left|\kappa \right|.$(10) The tangent angle determined by Iteration (4) applied to the surface model includes the influence of the approximation. Given that the curvature $\kappa$ in (9) is a local quantity, the value of $\theta$ obtained by (9) differs from that computed by Iteration (4) in addition to a finite value $\delta l$. Thus, the step length $\delta l$ should be constrained by considering the approximation error. Let ${\kappa }^{\mathrm{\prime }}$ be the plan curvature on $\mathit{x}+\delta \mathit{x}$, and ${\theta }^{\mathrm{\prime }}$ be the tangent angle estimated by ${\kappa }^{\mathrm{\prime }}$. Based on Figure 2, the distance between the two following points regulated by $\kappa$ and ${\kappa }^{\mathrm{\prime }}$ is as follows: (11) $2\delta l\left|\mathrm{s}\mathrm{i}\mathrm{n}{\displaystyle \frac{\theta -{\theta }^{\mathrm{\prime }}}{2}}\right|.$(11) It should be noted that the Taylor expansions of $\sin x$ and ${\sin }^{-1}x$ around $x=0$ are respectively expressed as follows: $\sin x=x+\mathrm{O}\left(x^3\right)\mathrm{a}\mathrm{n}\mathrm{d}{\sin }^{-1}x=x+\mathrm{O}\left(x^3\right).$

Figure 2. Difference between tangent angles estimated by two curvatures.

Moreover, Equation (11) is evaluated in the following form: $2\delta l\left|\mathrm{s}\mathrm{i}\mathrm{n}{\displaystyle \frac{\theta -{\theta }^{\mathrm{\prime }}}{2}}\right|=2\delta l\left|\mathrm{s}\mathrm{i}\mathrm{n}{\displaystyle \frac{{\sin }^{-1}\left({\displaystyle \frac{\kappa }{2}\delta l}\right)-{\sin }^{-1}\left({\displaystyle \frac{{\kappa }^{\mathrm{\prime }}}{2}\delta l}\right)}{2}}\right|={\displaystyle \frac{\delta l^2}{2}\left|\kappa -{\kappa }^{\mathrm{\prime }}\right|+\mathrm{O}\left(\delta l^4\right).}$To incorporate this distance into a standard, i.e. Tolx, the step length $\delta l$ is constrained with respect to the following: (12) $\delta l\le \sqrt{{\displaystyle \frac{2\cdot \mathrm{T}\mathrm{o}\mathrm{l}\mathrm{x}}{\left|\kappa -{\kappa }^{\mathrm{\prime }}\right|}}}.$(12) During advancement, aversion to both understepping and overstepping is essential, for which $\delta l$ is constantly adjusted. To circumvent understepping, tracing acceleration is attempted as follows: (13) $\delta l\leftarrow {\mathrm{c}}_1\cdot \delta l::({\mathrm{c}}_1>1).$(13) In contrast, overstepping is averted with the maximum step length, which is denoted by $\delta l_{\mathrm{M}}$. The step length is then restricted as follows: (14) $\delta l\le \delta l_{\mathrm{M}}.$(14) Combining (12), (13), and (14), the step length is elongated and constrained based on the following: (15) $\delta l\leftarrow min\left({\mathrm{c}}_1\delta l,\delta l_{\mathrm{M}},\sqrt{{\displaystyle \frac{2\cdot \mathrm{T}\mathrm{o}\mathrm{l}\mathrm{x}}{\left|\kappa -{\kappa }^{\mathrm{\prime }}\right|}}}\right).$(15) Once $\delta l$ is adjusted to follow (12), $\delta l$ is fixed and seeks $\theta$ until (7) and (8) hold. If such an $\theta$ cannot be determined within a certain number of iterations, $\delta l$ is then contracted to re-attempt determination by the following: (16) $\delta l\leftarrow {\mathrm{c}}_2\cdot \delta l\left(0<{\mathrm{c}}_2<1\right).$(16) Empirically, the constants ${\mathrm{c}}_1$ and ${\mathrm{c}}_2$ are set as 2 and 1/2, respectively.

Inspection of (4) implies potential divergence where $\left|\mathrm{\nabla }f\left(\mathit{x}+\delta \mathit{x}\right)\right|\ll 1$, that is, the following point is inadvertently positioned around a stationary point. To circumvent the stall around such critical points, contraction based on (16) is applied. Although empirical, this is performed unless the following is true: (17) $\left|\mathrm{\nabla }f\left(\mathit{x}+\delta \mathit{x}\right)\right|\delta l\ge \left|f\left(\mathit{x}+\delta \mathit{x}\right)\right|.$(17)

3.3. Identifying initial control points

To initiate tracking using the framework developed in the previous subsections, the determination of an initial point is required. This translates to identifying a solution, such that the following is true: $f\left(x,y\right)=\mathrm{h}.$A straightforward method is to identify all solutions on the edges across the entire region, which is possible if the elevation function is a two-dimensional polynomial. Solutions on the edges perpendicular to the y-axis are obtained for a fixed constant y, and those on the x-axis are obtained for a fixed constant x. For the bicubic spline function (2), the first task is reduced to identifying the real roots of a cubic equation for x. Although analytical methods are available, a numerical method was adopted to achieve faster computation. Computing the eigenvalues of the companion matrices (Edelman and Murakami 1995) is a common technique. All real roots on both the horizontal and vertical edges of the entire mesh are obtained as candidate initial points.

One notable difference from piecewise linear tessellations is that only a single initial point is sufficient to identify a single contour object. This may render multiple remaining points redundant, and those on the same contour should be eliminated to circumvent duplicate detection. A technique based on piecewise arc approximation was considered in this study.

Another difference from piecewise linearity is the challenge associated with identifying complete sets of objects in the entire region. Although the majority of contours intersect with edges, there may be small-sized closures that do not intersect with any of the edges. These cannot be detected only with initial points on the original edges for which the disaggregation of a patch is necessary for the search for intersections.

3.4. Gleaning redundant initial points

When a single initial point is given, no other point on the same contour object is required, which is in contrast to piecewise linear tessellations. During single-step advancement, points proximal to the approximated arc are identified to eliminate from the list of candidate initial points. The concept of the overhang vector (Goto 2021), among others, is utilized to achieve this.

Let ${\mathit{x}}_i$ and ${\mathit{x}}_{i+1}$ be two control endpoints. Consider an approximated arc and candidate initial point ${\mathit{p}}_j$, as shown in Figure 3. Let $\mathit{q}\equiv {\mathit{x}}_{i+1}-{\mathit{x}}_i$ and $\mathit{r}\equiv {\mathit{p}}_j-{\mathit{x}}_i$. The orthogonal projection of $\mathit{r}$ onto $\mathit{q}$ is $s\cdot \mathit{q}$, where the following is true: $s={\mathit{r}}^{\mathrm{T}}\mathit{q}/{\mathit{q}}^{\mathrm{T}}\mathit{q}.$If $0\le s\le 1$, the projection of ${\mathit{p}}_j$ to segment $\mathit{q}$ is on segment $\left[{\mathit{x}}_i,{\mathit{x}}_{i+1}\right]$. An overhang vector for $\mathit{r}$, as denoted by ${\mathit{h}}_r$, is introduced to quantify the distance from the segment as follows: (18) ${\mathit{v}}_r=\mathit{r}-s\mathit{q}=\mathit{Q}r,\mathit{Q}=\mathit{I}-{\displaystyle \frac{\mathit{q}{\mathit{q}}^{\mathrm{T}}}{{\mathit{q}}^{\mathrm{T}}\mathit{q}}.}$(18) The overhang of the approximated arc is obtained geometrically. Denoting this by ${\mathit{v}}_a$, it is quantified by the following: (19) ${\mathit{v}}_a={\displaystyle \frac{\left|\mathit{q}\right|}{2\mathrm{s}\mathrm{i}\mathrm{n}\theta }\left(\cos \phi -\cos \theta \right){\displaystyle \frac{\mathit{Q}r}{\left|\mathit{Q}r\right|}\left(-\theta \le \phi \le \theta \right),}}$(19) where $\cos \theta ={\displaystyle \frac{{\mathit{t}}^{\mathrm{T}}\mathit{q}}{\left|\mathit{q}\right|},\sin \phi =\left(2s-1\right)\sin \theta .}$In the above, $\mathit{t}$ is the unit tangent vector, and $\left|\mathit{t}\right|=1$.

Figure 3. Gleaning points proximal to an arc portion.

Let Tolh ($\ge$Tolx) be the tolerance for the overhangs. Subsequently, the set of proximal points ${\mathit{p}}_j$ that satisfies. (20) $\left|{\mathit{v}}_r-{\mathit{v}}_a\right|\le \mathrm{T}\mathrm{o}\mathrm{l}\mathrm{h},$(20) is obtained to eliminate them from the list of initial points. Besides the inspection by the overhang, points proximal to endpoints are obtained in the case of $s\simeq 0$ and $s\notin \left[0,1\right]$. They are identified by the following: (21) $\left|\mathit{r}\right|\le \mathrm{T}\mathrm{o}\mathrm{l}\mathrm{x}.$(21) The two tolerances, Tolx and Tolh, reflect the extent of the minimum computational resolution limit and the maximum error for piecewise arc approximation, respectively.

3.5. Termination or discontinuation

If the current control point approaches either the initial point or any of the borders of the analysis region, the step length may be adjusted to circumvent passing over the initial point or exceeding the border.

Let ${\mathit{x}}_0$ be the initial point. When the current step length is constrained by 3.1 follows: $\delta l>\left|{\mathit{x}}_i-{\mathit{x}}_0\right|,$whether or not the next point will potentially return to ${\mathit{x}}_0$, the step length is shortened according to $\delta l\leftarrow \left|{\mathit{x}}_i-{\mathit{x}}_0\right|,$during preparation for potential closure. In conjunction with this adaptation, if the following point is identical to ${\mathit{x}}_0$, then a closed contour is detected, and the advancement is terminated.

Let $\mathrm{\Omega }$ and $d\left(\mathit{x},\mathrm{\Omega }\right)$ be the region of analysis and distance from $\mathit{x}$ to $\mathrm{\Omega }$, respectively. If the following is true: $\delta l>d\left(\mathit{x},\mathrm{\Omega }\right),$$\delta l$ is contracted by the following: $\delta l\leftarrow d\left(\mathit{x},\mathrm{\Omega }\right),$to avoid inadvertently exceeding the border. Further, if the following is true: (22) $d\left(\mathit{x},\mathrm{\Omega }\right)<\mathrm{T}\mathrm{o}\mathrm{l}\mathrm{x},$(22) the advancement is discontinued, assuming the control point reaches a border. The tracking is then reversed, and advancement is continued until the control point reaches another border. An open contour is detected.

3.6. Identifying small contours inside a patch

In the proposed algorithm, candidate initial points are searched for only at the edges of the original mesh. Small-sized closures remain undetected if they do not intersect with any of the edges. To identify small loci, $f\left(x,y\right)=0$ for bicubic spline surface functions, and patch disaggregation is conducted for closer inspection. Given that the identification of a potential solution is intractable, the objective is to confirm non-existence, thus averting excessive disaggregation.

For simplicity, in this study, we investigated a square patch with dimensions of $1\times 1$. This assumption does not lose generality. In conjunction with scaling, the same inspection is applicable to rectangular and/or disaggregated patches. Thereafter, Eq. (2(2) $f=\sum _{j=0}^3\sum _{i=0}^3{a}_{ij}{x}^i{y}^j.$(2) ) can be re-written in matrix form as follows: (23) $f\left(x,y\right)={\mathit{x}}^{\mathrm{T}}A\mathit{y},$(23) where ${\mathit{x}}^{\mathrm{T}}=\left[1x{x}^2{x}^3\right]$ and $\mathit{y}=[1y{y}^2{y}^3{]}^{\mathrm{T}}$. If there is no intersection at edges $x=\left\{0,1\right\}$ and $y=\left\{0,1\right\}$, the following is true on $\mathrm{\forall }\left(x,y\right)\in \left(\left\{0,1\right\},0\le y\le 1\right),\left(0\le x\le 1,\left\{0,1\right\}\right)$: $f\left(x,y\right)>0\mathrm{o}\mathrm{r}f\left(x,y\right)<0.$Hereafter, the inspection is limited to $f\left(x,y\right)>0$ because the same framework is applicable to the opposite case.

Given a fixed $y_0\left(0\le y_0\le 1\right)$, $f\left(x,y\right)$ is a cubic function with respect to $x$ and has a maximum of two extremums. The locations of the potential extremums are the real roots of the quadratic equation: (24) $\mathit{x}{{}^{\mathrm{\prime }}}^{\mathrm{T}}A{\mathit{y}}_0=0,$(24) where $\mathit{x}{{}^{\mathrm{\prime }}}^{\mathrm{T}}=\left[012x3x^2\right]$ and ${\mathit{y}}_0=[1y_0{y}_0^2{y}_0^3{]}^{\mathrm{T}}$. If the real roots are in the range of $\left[0,1\right]$, the positiveness of the values on the points is checked. This is performed for the bottom and top borders, i.e. $y_0=0,1$. The formal representation is as follows: (25) $f\left(x,y\right)>0,x\in \left\{0,1\right\},$(25) for $y=0,1$ and real root(s) between $0\le y\le 1$ in (24). In an analogous manner, the same type of investigation is conducted for the left and right borders, where $x_0=0,1$. If any of the conditions do not hold, intersections exist such that $f\left(x,y\right)=0$, and the patch is not disaggregated.

Further conditions for ensuring non-existence rely on a tangent surface on the borders. The tangent surface on a fixed $x=x_0$ is expressed as follows: $g\left(x,y:;x=x_0\right)\equiv f_x\left(x_0,y\right)\left(x-x_0\right)+f\left(x_0,y\right)=[{x^{\mathrm{\prime }}|}_{x=x_0}\left(x-x_0\right)+\mathit{x}{]}^{\mathrm{T}}A\mathit{y},$where $x^{\mathrm{\prime }}{|}_{x=x_0}^{\mathrm{T}}=\left[012x_{0}3x_0^2\right]$. If the elevation of this surface based on the left border at $x_0=0$ is constantly above 0 on the opposite side, i.e. $x=1$, then the following is true: (26) $g\left(1,y;x=0\right)=\left[\left[0100\right]+\left[1000\right]\right]A\mathit{y}=\left[1100\right]A\mathit{y}>0.$(26) Conversely, considering the tangent surface on the opposite border at $x_0=1$, the elevation on the left border at $x=0$ is expressed as follows: (27) $g\left(0,y;x=1\right)=\left[\left[0123\right]\left(-1\right)+\left[1111\right]\right]A\mathit{y}=\left[10-1-2\right]A\mathit{y}>0.$(27) Similar to the investigation of $f\left(x,y\right)>0$ above, Conditions (26) and (27) for all $y\left(0\le y\le 1\right)$ are checked. The left-hand sides are cubic polynomials with respect to $y$, for which the positiveness can be checked in a manner analogous to that for $f\left(x,y\right)$. The corresponding derivatives with respect to $y$ are those multiplied by ${\mathit{y}}^{\mathrm{\prime }}=[012y3y^2{]}^{\mathrm{T}}$ to replace $\mathit{y}$. For example, the stationary points on the left-hand side of Eq. (26(26) $g\left(1,y;x=0\right)=\left[\left[0100\right]+\left[1000\right]\right]A\mathit{y}=\left[1100\right]A\mathit{y}>0.$(26) ) are located at the roots of the following: (28) $\left[1100\right]A[012y3y^2{]}^{\mathrm{T}}=0.$(28) The elevations are checked only if there are real roots in the range of $0\le y\le 1$. In addition to the two conditions of $g\left(1,y:;0\right)>0$ at $y=0$ and $y=1$, if the following is true for $y=0,1$: (29) $g\left(1,y:;x=0\right)>0,$(29) and the real root(s) of (28) between $0\le y\le 1$, the necessary conditions for non-existence are then satisfied.

Only Conditions (25) and (29) do not guarantee non-existence. Figure 4 illustrates the three scenarios of tangent surfaces on the left border, where $x=0$. Moreover, $f\left(0,y\right)>0$ and $g\left(1,y:;x=0\right)>0$ are true for all cases, and Cases (b) and (c) have intersections at $f\left(x,y\right)=0$ in the range of $0\le x\le 1$. Such cases occur when the tangent point is on a concave surface. To avert misidentification, another tangent surface on the opposite border, i.e. (27), should be verified: $g\left(0,y:;x=1\right)>0$ does not hold, whereas $f\left(1,y\right)>0$ holds. All the inspections above are conducted for the left and right borders, that is, $x=\left\{0,1\right\}$, in addition to the bottom and top borders, that is, $y=\left\{0,1\right\}$.

Figure 4. Existence of cubic function roots.

If all conditions are not satisfied, the non-existence of a solution inside the patch is assured, and further disaggregation is discontinued. Otherwise, the current patch is disaggregated to search for potential closures. The simplest type of disaggregation was adopted in this study, where the square patch was disaggregated into four square patches divided by two borders, that is, $x=y=1/2$. The same inspection is performed for each disaggregated patch. If an intersection $f\left(x,y\right)=0$ is found on the borders, further disaggregation is skipped. This is valid because the bicubic spline function in (23) cannot contain two closures inside the original patch. Thus, disaggregation is conducted recursively until the edge length contracts below the standard.

Figure 5 presents an example of the disaggregation of the original patch with dimensions of $1\times 1$. As several positive conditions do not follow the outermost borders, the original patch is subdivided into four pieces. All conditions are then satisfied for the upper two patches, which ensures the non-existence of a solution, and no further disaggregation is conducted. In contrast, the lower two patches are disaggregated into four $\left(1/4\right)\times \left(1/4\right)$-sized patches. With respect to the eight $\left(1/8\right)\times \left(1/8\right)$-sized patches in the lower right, intersections are found on the edges, for which the recursive subdivision is discontinued, and the coordinates of the intersections are added as candidate initial points. If the non-existence conditions continue to be satisfied, the recursive process is repeated until the edge length is below the standard Tolp. A smaller tolerance may aid in the identification of small-sized contour objects, with an impact on the computation time.

Figure 5. A small closure detected by the recursive disaggregation of patches.

3.7. Entire flow

The individual components considered above were organized for enhanced visibility and accessibility. The flow of the proposed framework is outlined in Figure 6. First, a surface model is determined using a digital elevation model dataset. The target elevation $h$ is subtracted to focus on solving $f\left(x,y\right)=0$. Intersections on the borders of the patches are then identified as candidate initial points, which is exhaustively conducted. If non-existence in a target patch cannot be assured, it is disaggregated for further inspection. After identification, an initial point ${\mathit{x}}_0$ is selected from the list of points $\mathrm{\Psi }$, which required no elaboration. The selection of the first point is sufficient. The tracking process is then initiated. It should be noted that step-length adaptation requires elaboration considering both the curvature error and tolerances. When the step length $\delta l$ is determined, the identified points in $\mathrm{\Psi }$ that lie on the same contour are eliminated and canceled using as initial points. The control point $\mathit{x}$ is advanced, which is repeated until it returns to the initial point ${\mathit{x}}_0$ or encounters a boundary of the target region $\mathrm{\partial }\mathrm{\Omega }$. In each advancement, step length elongation is applied where and if possible.

Figure 6. Proposed algorithm for identifying complete sets of contour objects.

The procedures described above are used to detect and identify a complete set of contour objects for a single target elevation. Most practical outlets are accompanied by repetitive applications of this module for varying elevations, which are typically at a constant interval.

4. Demonstrative experiments

4.1. Study area

An isolated island located in the south of central Tokyo, Japan, was selected as the study area. The island is referred to as Izu-Oshima, and its location is shown in Figure 7. It has an area of approximately 91.06 $\mathrm{k}{\mathrm{m}}^2$, and the average annual precipitation and temperature during 1991–2020 were 2858 mm and 16.4${}^{\circ }\mathrm{C}$, respectively. At the center of the island is an active volcano, namely, Mt. Miharayama, with an elevation at the summit of 764 m. A major eruption occurred in 1986, in addition to volcanic tremors and minor eruptions. There is a large caldera vent around the center along this line.

Figure 7. Location of Izu-Oshima island, Tokyo, Japan.

4.2. Experiment conditions

In this study, a digital elevation dataset was used, namely, ‘Digital Map 50 m Grid (Elevation)’, which was published by the Geospatial Information Authority of Japan. The mesh size was 2.25 s $\times$ 1.5 s in the longitude and latitude directions, respectively, which corresponds to approximately 57.24 m $\times$46.22 m. Although in recent years, higher resolution datasets, e.g. 0.4 s $\times$ 0.4 s, have been made available, coarser resolutions were used to exploit the significance of the proposed framework. The elevation of the ocean was set as −1 m to estimate the coastline by searching for contours with elevations of h = 0 m. Contour detection was conducted for elevations varying from 0–750 m at intervals of 10 m or 25 m.

A personal computer (PC) equipped with an Intel® Core™ i7-1185G7, 3.00 GHz processor, and 32 GB of random access memory (RAM) running Windows 11 Professional was used. The associated program codes were implemented in MATLAB R2023a.

4.3. Contour detection

The cartographic contour extraction was performed in a straightforward manner. Figure 8 presents an overview of the outlet detected using the proposed framework, in which the contours were colored in varying gradations from blue to red. The elevation interval was 10 m. Several small contours were highlighted in red, thus indicating the presence of a depression inside the region. All extracted objects were oriented, for which convexity or concavity was assessable.

Figure 8. Contours identified using the proposed framework (Framework A).

Figure 9 presents the outlet computed using the preceding algorithm developed in Goto and Shimakawa (2016). The tracking algorithm considers the roundness of the curve, which is evaluated using three curvature-associated quantities. Although the smoothness of curves and storage efficiency are core advantages, the method to obtain initial points requires improvement, and the candidate points have been sought on a piecewise linear surface model. A direct comparison of Figure 8 and Figure 9 reveals a difference: a greater number of depression contours in red can be observed in Figure 8. There are steep cliffs and jaggy coastlines around the northern and southwestern areas, as designated by magenta boxes. Figure 8 reveals small contours around the reefs, which are not present in Figure 9.

Figure 9. Contours identified using the method established by Goto and Shimakawa (Goto and Shimakawa 2016) (Framework B).

To observe the significance of this application, four types of contour identification were tested: (A) the established framework, (B) a framework focusing on roundness developed by Goto and Shimakawa (2016), (C) a piecewise triangular tessellation, and (D) a bilinear spline model. Table 1 reports the profiles of the objects reproduced by the four frameworks with elevation intervals of 25 m, in addition to the experimental parameters. For the first two, namely (A) and (B), the number of independent polygons, duplicate detections, total number of control points, and average step length are reported. For the latter two frameworks, (C) and (D), all initial points were used to link each other, for which there were no redundant initial points.

Table 1. Contour profiles identified using the three frameworks.

A duplicate detection indicates that the identified contour object was already detected and listed by a preceding process initiated from a different initial point. This was deemed true if the area, perimeter, and centroid of an identified object were all close within a threshold to an entry in the list of independent objects. Along with this standard, Framework (B) encountered frequent duplications in various elevations, while the proposed framework (A) did not exhibit such inefficiency. This enhancement was achieved by piecewise arc approximation between successive control points, along with the elimination technique of initial points established in Section 3.4. Although duplicate detection might occur with different parameters – for instance, tighter tolerances Tolx and Tolh – the current configuration was able to avert it.

It should be noted that the number of objects identified by (A) was greater than that identified by (B) at most elevations. For example, with respect to an elevation of h = 0, Framework (A) detected 67 independent contours, whereas only two were detected using Framework (B). The main difference was the small contours distributed around the northern and western shores. At an elevation of h = 750 m, Framework (A) detected one small contour, whereas (B) did not. This is because of the non-existence of an intersection in the piecewise linear surface for Framework (C), whereas Framework (B) was dependent on identifying the initial points. Although the maximum step length $\delta l_{\mathrm{M}}$ was set as the same value for Frameworks (A) and (B), the average step length in Framework (A) was half of that in Framework (B). Consequently, this enhanced the identification of the other initial points along the trajectory. For Frameworks (C) and (D), the step length was restricted by the sizes of the triangular and rectangular patches. The length of the former was slightly smaller than half of the longer edge length, while the length of the latter was slightly smaller than the shorter edge length.

To inspect the computational efficiency in Framework (A), the computation times were measured for varying tolerance for patch disaggregation, namely Tolp. The times were measured by the ‘cputime’ function, with which the resulting time may fluctuate by each trial depending on the CPU load. Table 2 reports the computation times in seconds consumed by three parts: (1) finding initial points on the original mesh, (2) those on disaggregated patches, and (3) contour tracking initiated from a detected initial point. The depth of disaggregation indicates the maximum number of recursions to find the root(s) or confirm non-existence. The elapsed times for finding initial points on the original mesh were shorter than 1 min in all cases, while those for contour tracking, including gleaning of redundant initial points, were approximately 20 s. By contrast, the consumed times to find initial points on disaggregated patches were shorter than 2 s, even with the smallest Tolp = 0.001 m. The number of detected contours was 264, with the largest Tolp = 10 m, while those with smaller ones were all 266: two small contours were detected with Tolp $\le$1 m. Even with a small Tolp, the computation time in total would not increase significantly, and the target precision should rather be determined by the resolution of the dataset. Considering the mesh size, Tolp = 1 m may suffice in this experiment.

Table 2. Computation times with varying tolerances for patch disaggregation.

Tolerance: Tolp (m)	10	1	0.1	0.01	0.001
Computation time (sec)
(1) Original mesh	52.95	52.38	52.30	51.22	51.36
(2) Disaggregated patches	1.33	1.20	1.36	1.33	1.83
(3) Contour tracking	22.91	20.33	19.97	19.95	18.84
Depth of disaggregation	2	3	7	9	14
Number of detected contours	264	266	266	266	266

A further comparison of the four frameworks was conducted. Figures 10 and 11 present the magnified outlets by the four frameworks focused on the northern and western reefs demarcated in Figure 8 with elevation intervals of 25 m. The red contours reflect the depression inside each object. In Figure 10, the runway of Oshima Airport is visible in the lower left, where no contour was identified. From a detailed inspection of Frameworks (A) and (B), it can be seen that both exhibited smooth curves, and only the former detected several small objects in the northern reef of the coastline. To compare the proposed framework with traditional frameworks, oriented contour tracking algorithms were applied to piecewise triangular and bilinear spline surface models in Frameworks (C) and (D), respectively. In contrast to the bicubic spline model, red objects indicating depressions emerged in the lower part, which would relate to the approximation error inevitable in lower-order surface models. Further, a close inspection of Framework (C) reveals a notable modeling error on the headland in the west northernmost. While only one narrow hill was detected in Frameworks (A), (B), and (D), it is inadvertently separated into two hills in Framework (C) due to the ambiguity specific to piecewise triangular modeling.

Figure 10. Identified contours around the northern reef using the four methods.

Figure 11. Identified contours around the southwestern reef using the four methods.

A comparison of Frameworks (A) and (B), as shown in Figure 11, reveals a trend similar to that in Figure 10. Framework (A) could identify small objects along the western coastline compared with Framework (B), and Framework (C) reproduced jagged coastlines compared with the other ones. In Framework (D), a larger number of red objects were detected along a gutter toward the southwest, yet it would be deemed an approximation error because the elevation would decrease in a monotonic manner along a runoff or stream. In fact, most contours around the gutter are receded and not closed in Frameworks (A) and (B).

4.4. Three-dimensional consistent contour-based network

Another application of the proposed framework was tested. A three-dimensional manipulation framework for a bicubic spline surface model established in Goto (2021) was applied. It should be noted that streamlined tracing along the steepest upward and downward directions aids in reproducing a three-dimensionally consistent contour-based network (3D-CBN). The framework proposed in this study was used as an external module. Figure 12 presents an overview of the outlook of the 3D-CBN, whereas Figure 13 presents a magnification of the cliffs in the north and the scoria cone in the south, as demarcated in Figure 12. Frameworks (A) and (B) were assessed for comparison, and the outlets are delineated with elevation intervals of 10 and 25 m in the northern and southern outlets, respectively. The green objects indicate contours identified by the proposed framework, whereas cyan and magenta objects indicate the steepest upward and downward streamlines, respectively. In the upper CBNs, two major steep hills were observed in the left and right parts close to the coastline, whereas the upward streamlines converged to local maxima without oscillation. Streamline tracing is continued until it either encounters a stationary point or reaches the lower or upper limit of the elevation interval. Upon reaching such a limit, the tracing discontinues to evaluate if multiple streamlines are proximal or confluent. This may ensue that, if the contour detection is not complete, tracing may terminate inadvertently. In the upper outlet by Framework (B), such a flaw was observed in the right part demarcated by the yellow frame. By contrast, the reproduced CBN by Framework (A) is complete and consistent topologically.

Figure 12. Overview of the three-dimensionally consistent contour-based network.

Figure 13. Detailed contour-based network around the cliffs and scoria cone.

In the lower outlets of the scoria cone, the green contours are roundish, from which the steepest upward streamlines moved toward the top of the hill. There seemed little difference between Frameworks (A) and (B) in this area. It should be noted that the mesh size was approximately the quarter length of the scalebar, and the peak location of the cone could be identified in a precise manner in conjunction with the underlaid aerial photo. A driving school was observed at the top, and no contour was detected around it. This may indicate that the premises are built on a flat area.

4.5. Slope analysis

The proposed framework was utilized for further application to slope analysis. Given a continuous elevation function of $z=f\left(x,y\right)$, the slope was quantified by the magnitude of the gradient vector, i.e. $\left|\mathrm{\nabla }f\right|=\sqrt{f_x^2+f_y^2}$. The interpolated values of $\left|\mathrm{\nabla }f\right|$ on a regular rectangular lattice were provided as elevation values.

Figure 14 delineates the outlet of the slope analysis with a color gradation from green to red, in which the isovalues varied in intervals of 0.1. The values of $\left|\mathrm{\nabla }f\right|$ = 0.1, 0.5, and 1.0 correspond to approximately 5.7, 26.6, and 45.0 degrees, respectively. Red objects emanated around the northern, western, and eastern coastlines. The two regions in the northern cliffs and the scoria cone were magnified for closer inspection, and Frameworks (A) and (B) were again assessed for comparison. The cyan curves in Figure 15 correspond to geographical contours with an elevation interval of 25 m. The solid and dashed curves indicate that the values inside the curve are greater, i.e. steeper (respectively smaller and less steep) than those on the contour.

Figure 14. Overview of the island based on the slope analysis.

Figure 15. Detailed slope analysis around the cliffs and scoria cone.

Along the northern coastline, the red curves existed, which exhibit the steepness of the cliffs. While relative steepness might be assessed by observing the density of contours, the outlet delineated in Figure 15 can quantitatively detail the extent. In view of the upper two outlets, successive closures on the left side of the runway signify that the shape of the hill is acute. In contrast, dashed contours were observed to the west of the golf course, for which the slope turned gentle toward the ridge. While Frameworks (A) and (B) might have reproduced similar outlets, there were differences in the ability to detect small closures: the green and light green contours detected by Framework (A) designated by red arrows did not exist in the outlet by Framework (B).

In the lower two results of the scoria cone, the cyan curves indicate again elevations at intervals of 25 m. The height of the cone is thus some 150 m from the base. Approximately 100 m north of the top of the hill, there was a series of solid brown contours, which implies the existence of a locally steep slope. The yellow ‘8’-shaped dashed contour located south of the hilltop indicates the slope is gentler than the surrounding locations. Along these lines, a series of solid objects implies that the surface is acute, whereas those of the dashed lines imply that the slope turns gentle toward the inside. The ability to detect small closures by Framework (A) was significant in the lower outlets. The green and light green closures designated by red arrows were not detected by Framework (B). It is worth noting that Framework (B) relies on piecewise triangular tessellation to find candidate initial points, in which greater approximation errors might be contained in the surface model. The slope measure $\left|\mathrm{\nabla }f\right|$ includes differentiations of the elevation function $f\left(x,y\right)$, for which it would vary more sharply than the elevation itself. While isolines for this type of quantity would contain many small closures, the established Framework (A) provided a complete detection scheme suited for bicubic spline models.

The outlet comprises a set of signed polygons and not a raster image. Current popular analyses using geographic information system (GIS) functionality are inextricably related to raster analysis, for which the outlet turned jagged when magnified. In contrast, by utilizing the proposed framework, the method of slope analysis was transferred to vector- or polygon-based; the overlay of other vector-based objects, including geographical contours, would evolve analytical techniques to exploit better insight: an additional option for delineation is gained.

5. Conclusions

The core objective of contour detection is to identify loci $\left(x,y\right)$ such that $f\left(x,y\right)=0$. Piecewise linear tessellation provides a straightforward solution to the resulting commodity in GIS-based analyses. All intercepts are determined and used as control points that connect two neighboring points. Its structural simplicity simplifies the calculation and delineation objectives, thus yielding byproducts of jagged outlets when magnified. In contrast, bicubic spline tessellation presents a challenge with respect to both theoretical and computational aspects. In addition to the non-linearity of the elevation function, the number of independent trajectories on and in a unit patch is yet to be determined.

Given a fixed constant x or y, interceptions can be computed by solving a cubic polynomial. A preliminary tracking method on bicubic spline surfaces was established in Goto and Shimakawa (2016), in addition to a detailed formalization of curvature-associated quantities. It is yet to be determined which interceptions should be used as initial points to identify independent objects, thus averting duplicate detection if possible.

Significant achievements made in this study include the complete detection of candidate initial points and the elimination of those lying on the same contour object. Mathematically, the former problem is intractable when confirming the existence of a solution on or inside a region. In this study, a computational inspection method for verifying the non-existence of a solution was developed. In particular, if non-existence cannot be ensured in a current region, it is subdivided into four portions, which are inspected recursively. Tolerance is essential to circumvent excessive disaggregation significantly below the precision of the original dataset.

The elimination of the initial points on the same contour curve is dependent on the buffered elimination along the approximated arc portion. To attain a higher elimination rate, a slight adjustment of the step length is necessary, and excessive contraction results in computational and storage inefficiencies, whereas improper elongation fails to eliminate the initial points. Tolerance is critical in this procedure. In particular, for the comparison of local curvatures at two distinct points, the step length is adapted to contain the error by piecewise arc approximation into the standard, i.e. tolerance. With the exhaustive detection of the initial points, this method is invaluable for the improved prevention of the redundant detection of identical objects. The experimental results in Table 1 provide a practical benefit: no duplicate detection of the same object was observed.

The evolution of contour detection for bicubic spline surfaces, particularly concerning the completeness of identified objects, provides new possibilities for numerous applications to GIS-based analyses in various fields, namely, engineering, geography, and earth science. The two demonstrative examples address applications for generating 3D-CBN and slope analysis. The former experiment exhibited significant affinity with streamline tracking: the smoothness of the identified curves was significant, and the location identification of ridges and local peaks was well-demonstrated.

The other application to slope analysis provided richer possibilities and insightful outlets, particularly via detailed delineations. In contrast to the current slope analyses, which are inextricably related to raster-based assessment, the fundamental method of analysis has transitioned to vector-based assessment. In addition to quantifying the extent of steepness and identifying the tracking direction, the detected object is accompanied by the orientation of the closure. The sign of a closed object provides a measure of whether the slope inside the object is greater or smaller. Inside the contours with negative and smallest absolute values, the terrain around the region may be considered almost flat. Although the steepness may be assessed upon observing the density of geographical contour curves, this is unsuitable for quantitative assessments. To date, slope analysis has been employed for this objective; however, a typical or averaged value is associated with a single patch, for which detailed analysis was ruled out, where only previous or coarse datasets are available. The final experiment demonstrated a promising application that was addressed and implemented in raster-based analysis contexts.

The proposed framework was developed on a mathematical basis, for which the completeness of object detection, including small-sized closures, was achieved. Moreover, storage and computational efficiencies were minimized when compared with the framework proposed in Goto and Shimakawa (2016). This is a byproduct of completeness, and further advancements with regard to inspecting the non-existence or existence of a solution on or inside a rectangular region will be conducted in future research.

Disclosure statement

No potential conflict of interest was reported by the author(s).

This research was funded by the JSPS KAKENHI Grant Number 21K01021.

Data availability statement

The datasets that support the findings of this study are openly available in ‘Mendeley Data' at http://dx.doi.org/10.17632/f4hfsfrsy4.1. All of them are preserved in the ESRI’s Shapefile format.

Additional information

Funding

This work was supported by Japan Society for the Promotion of Science [grant number 21K01021].