A combined contour lines iteration algorithm and Delaunay triangulation for terrain modeling enhancement

ABSTRACT Digital Elevation Models (DEMs) play a crucial role in civil and environmental applications, such as hydrologic and geologic analyses, hazard monitoring, natural resources exploration, etc. Generally, DEMs can be generated from various data sources, such as ground surveys, photogrammetric stereo methods, satellite images, laser scanning, and digitized contour lines. Compared with other data sources, contour lines are still the cheapest and more common data source becausethey cover all areas, at different scales, in most countries. Although there are different algorithms and technologies for interpolation in between contour lines, DEMs extracted solely from contours still suffer from poor terrain quality representation, which in turn negatively affects the quality of analytical applications results. In this paper, an approach for improving the digital terrain modeling based on contour line densification and Delaunay triangulation is presented to acquire a more suitable DEM for hydrographic modeling and its applications. The proposed methodology was tested using a variety of terrain patterns in terms of intensity: hilly, undulated, and plain (1:25,000 topographic map, 5 m contour interval). The precision of the extracted GRID model increases as the number of added contours increases. Adding four contour lines, the Root Mean Square Error (RMSE) of examining points were 0.26 m, 0.29 m, and 0.05 m for hilly, undulated, and plain samples, respectively, and the Mean Absolute Error (MAE) were 0.50 m, 0.48 m, and 0.17 m. The convergence probabilities between extracted and original flow lines for the same regions were 96.91%, 94.93%, and 84.03%. Applying the methodology, experimental results indicate that the developed approach provides a significant advantage in terrain modeling enhancement, generates DEMs smoothly and effectively from contours, mitigates problems and reduces uncertainties.


Introduction
Digital Elevation Model (DEM) is an array representation of squared pixels; each pixel has associated with its elevation data.It is considered as the most prevalent and widespread model that can be used in many civil and environmental applications, such as geological studies (Liu et al. 2021;Ricchetti 2001;Yang, Meng, and Zhang 2011), hydrological modeling (Hopkinson, Hayashi, and Peddle 2008;Le Coz et al. 2009;Li et al. 2017), disaster analysis (Demirkesen et al. 2007;Tsai et al. 2010;Xu and Dong 2021), and agriculture applications (Tijskens, Ramon, and De Baerdemaeker 2003;Chidi et al. 2021).DEM generation with various data sources has attracted a lot of attention in the research community (Shen, Meng, and Zhang 2016;Zhou and Zhu 2013;Kamel, Miky, and El Shouny 2020;Bhushan et al. 2021;Habib 2021).Besides the digitized contour lines method, DEMs can be generated from various data sources; ground surveys, photogrammetric stereo methods, satellite images, laser scanning, radar interferometry, and Global Positioning System (GPS) (Wang, Holland, and Gudmundsson 2018;Soni et al. 2021).Each data type has its methodology however, each methodology has its advantages and challenges.Indeed, the scale and level of detail required in an individual project will determine the spatial resolution and accuracy of the DEM required and so the suitable methodology to be followed.For very large-scale projects covering thousands of kilometers, coarse resolution DEMs with a spatial resolution of tens of meters are most appropriate.while at the very detailed level, a fine resolution that which typically has a spatial resolution of about one meter with high accuracies is needed.
Surveying works needed for the construction of high-resolution DEM spatially for small areas may be very expensive because of locational charges or even may not be available due to flying restrictions when using drowns for example.In such circumstances, DEMs from satellite platforms that have lower vertical accuracies but cover many large areas and are almost free of charges may be the best option.Although these platforms provide global coverage and their data is already used in various applications (Yang et al. 2015), they are suffering from great shortcomings associated with DEMs production.For example, the radar interferometry technique is challenged by high relief terrain and the distortion associated with the imaging geometry of side-looking radars (Kubanek, Westerhaus, and Heck 2015) while, the stereo-optical technique is challenged by cloud cover.On the other hand, accurate interpolation from contours is still an important way to generate DEMs in many areas; where high-resolution DEMs are not available (Taud, Parrot, and Alvarez 1999;Werbrouck et al. 2011), in cases where historical contour maps are used to evaluate topographic change (Carley et al. 2012), and for the generation of regional DEMs.Compared with other data sources, contour lines are still the cheapest and more common data source because, they cover all areas at different scales in most countries.
There are many interpolation algorithms adopted in Digital Terrain Models (DTMs) production from contour lines (Hutchinson 1988;Taud, Parrot, and Alvarez 1999;Ardiansyah and Yokoyama 2002;Gousie and Franklin 2003;Zheng et al. 2016;Li et al. 2017;Zhao 2020).However, most of these methods are heavily affected by the discontinuity or absence of altitude data between contour lines.This may cause systematic modeling errors in the generation of DEM.The appearance of a number of lines that represent imaginary topographic elements that are not identical to the terrain reality, and the appearance of pits and mock terrain are the most common errors that arise in grid terrain modeling.Furthermore, the inability to derive the local peak points; top of a hill (saddle) and bottom of a valley (perigee) which appear as completely flat areas.These errors have an impact on modeling runoff paths as they constitute interruptions in the water streams.Also, the lack of elevation information due to the long horizontal spacing between contour lines creates unreal features in the DEM, which in turn affects the vertical resolution of the DEM produced.Hence, the importance of removing these defects in order for the digital model to contribute to a reliable description of the earth's surface which will allow derivation of realistic streamlines along the topographic surface.It is possible through a perfect modeling procedure based on high-density contour lines to obtain more useful information for spatial analysis purposes than that provided by digital models with a high spatial resolution (Bennet 1991).
In this work, an approach for improving the digital terrain modeling to acquire a more suitable DEM for hydrographic modeling and its applications is presented.This proposed approach aims to mitigate the terrain modeling errors based on the concept of densification of elevation data represented by contour lines and produce smooth graduation of elevation values over smaller distances.That provides a structure of a digital terrain model capable of representing the morphological features without disappearing or interruption.
Briefly, this paper provided the following main contributions: • Improvement of the digital terrain modeling by proposing an approach that maintains a high statistical and visual accuracy to acquire a more suitable DEM for hydrographic modeling.The article is structured as follows.Section 2 gives a brief summary of the relevant work.Section 3 discusses the research methodology of the suggested solution.Results are presented in Section 4. Sections 5 and 6 address the discussion of results and conclusions, respectively.

Overview of contour line interpolation methods
To interpolate a point elevation from adjacent contour lines, a large variety of interpolation algorithms have been proposed or developed.The accuracy of the interpolated elevation is strongly impacted by the degree of terrain complexity and estimation method.The estimation method can be classified as fittedfunction methods (Lagrange interpolation, collocation, minimum curvature splines, Kriging, relaxation, approximation), distance-based methods (Inverse Distance Weighting (IDW), Inverse Distance Weighting Gradients (IDWG)), triangle-based methods (linear facets, nonlinear triangular patches), rectangle-based methods (bilinear patches, Hermite patches, Bezier patches, B-spline patches), and neighborhood-based methods (linear interpolation, nonlinear interpolation) (Watson 1992).According to Taud, Parrot, and Alvarez (1999), the above methods can be classified into two groups which are: analytic function-based methods and direct summation-based methods.Some researchers have demonstrated the strong potential of using direct summation-based methods (Guo et al. 2010;Li et al. 2017;Lu and Wong 2008;Soycan and Soycan 2009).
From another point of view, grid interpolation methods from contour lines can be classified into two types: point methods and line methods (Ardiansyah and Yokoyama 2002).In point methods, discrete vertices along contour lines are to be used to estimate the elevation of the unknown grids.The estimated value is mainly based on the characteristics of the point.Triangulated Irregular Network (TIN) and minimum curvature are the most popular interpolation point methods that can be used to generate DEM from contours.On the other hand, the line method uses all points along the contours, it depends on the line characteristics to estimate the grid elevation.This means that all flow lines, skeleton lines and also intermediate contour lines can help in DEM generation (Gousie and Franklin 2003).Watson (1992) has used the weighted moving average technique to interpolate the elevation value on a grid through the geometric relationships between adjacent contour lines.In this method, eight intersection points between both cardinal and intermediate directions rays are drawn at a certain point and their closest adjacent contour lines are identified, then the inverse distance weight method was used to compute the elevation value from these points.Taud, Parrot, and Alvarez (1999) presented a raster-based interpolation algorithm using an iterative procedure to produce an extension of contour lines by applying alternate fourand eight-connected erosions of the background until the generated surfaces become contiguous.However, this method gives an error when processing the inner data of a closed contour.Ardiansyah and Yokoyama (2002) have used a monotone interpolation provided by a simple cubic Hermit function to evaluate the elevation values of a slope grid points along the steepest slope segment chain.They have improved methodology for generation DEM from contours.The used method is characterized by 16 search directions to find the steepest slope from a grid point to the nearest and the second nearest contour line.Xie et al. (2003) developed a special kernel filter to insert the height value of steep slope areas.That kernel can interpolate a grid cell using the following procedure: in the beginning, the contours are rasterized as a grid.Then, three types are defined concerning the relationships between contour lines and cells.These types are "No value" for no contour line passing this cell, "Single value" for one contour line passing this cell, and "Multi value" when more than contour line can pass this cell.The elevation of cells defined as "Multi value" is calculated using the proposed filter, whereas those of "Single value" are directly assigned based on the appropriate contour elevation.The value of a given "No value" cell is interpolated from its neighboring cells.Wang, Gao, and Wu (2005) developed the way of estimating the elevation value of cells defined as "Single value" by considering the spatial relationship between the associated contours and their neighborhood.A triangulated irregular network TIN constructed using two adjacent contours is used to detect the triangle within which the "Single value" cell center falls, then the three vertices of this triangle are used to interpolate the elevation value.An algorithm for the generation of DEMs from contour lines considering geomorphic features was proposed by Rui et al. (2016).That algorithm is based on identifying the point on a contour line that is nearest to the interpolated point.Then, from the interpolated point 16 points in eight directions on neighboring contour were obtained.Virtual control points, according to various geomorphic features, were added in each direction to prevent the interpolation result from exceeding the contour interval.The local surface was constructed using the spline function.Zhao (2020) presented an indirect interpolation model that considers local terrain similarity when interpolating unknown points.
In ArcGis software, Topo to Raster tool is used to interpolate elevation values for a raster while imposing specific constraints.It is based on the ANUDEM algorithm developed by Hutchinson (Hutchinson 1988).This algorithm is designed specifically to create a surface that is more closely to the natural drainage surface.In this algorithm, information in the contours is used to find the points of local maximum of curvature in each contour.A network of curvilinear streams and ridges intersecting these points is then derived using the initial elevation grid.Then, an iterative way is used to accurately locate the locations of stream and ridge lines and update the DEM.Zheng et al. (2016) presented an improved ANUDEM method that combines topographic correction using a high correlation of topological structure between contour lines from multi-scale DEM.
Shortcomings arise when using contours to generate DEMs can be summarized as: • The interpolation in complex areas may yield improper results when it is done along the direction of the steepest slope.• Since it is determined only from the grid point to the nearest contour line, the direction of interpolation may not represent the steepest slope.• The misidentification of slope grid points or the steepest slope segment cannot be calculated correctly.
• The interpolated data in the inner exhibits the same elevation value as the closed contour line.• Many imaginary flat surfaces with a slope equal to zero appeared.• Wrong formation of the triangles whose points were either built on the same contour line or on two contour lines withthe same height value.
• The inability to detect the local peak points which appear as completely flat areas.

Basic principle of fractional calculus
By standard definition, a triangulation of a polygon is a Delaunay triangulation if every triangle meets the empty-circle condition to maximize the internal minimum angles.This means that, given a triangle (P i ,P j , P k ) belonging to a Delaunay triangulation of a set of generating input points, no other point of the generating set is internal to the circle defined by P i ,P j ,P k .This guarantees that the nearest points are joined together.Also, in the Delaunay triangulation, the minimum angle of all the triangles is the maximum outside of all the triangulations.It means that the triangles are as "equilateral" as possible so that distal information is not linked together (Cao 2015).However, the Delaunay triangulation cannot guarantee that the preexisting segments like contours and line of maximum slope segments are preserved.The Delaunay triangulation can be applied to construct a TIN from a set of points coming from the vectorization of contour lines and isolated points of a topographic.When a TIN is generated from topographic contour lines, in some cases, a number of false flat triangles could be generated where nodes of equal height (belonging to the same contour are joined, forming horizontal artificial terraces along the contour line (Stevens, Manville, andHeron 2003, 2003).The flat regions are particularly liable to occur where the curvature of the contour lines is high.Also, the Delaunay triangulation can join points belonging to different noncontiguous contour lines, crossing the intermediate contour line.The main defect of Delaunay triangulation is that it does not force given segments (faults or other morphological features, e.g.ridges) to become edges of triangles (breaklines).
The concept of constrained Delaunay triangulation is introduced in order to prevent some false features from appearing.It can keep all preexisting segments in the model when relaxing the empty-circle criterion around restrained edges.This relaxation rule allows the circumcircles of triangles containing the preexisting edge to include other points.Obviously, a constrained Delaunay triangulation may not always satisfy the criteria for a Delaunay triangulation.Many available algorithms in the literature are designed to overcome the problem of flat areas (Eastman 1995 ;Favalli 2004).

Tools and mathematical notions
The contour lines in a topographic map are mathematically known as level curves which are curves of constant elevation above mean sea level.The level curves of a function f of two variables are the curves with equations f where k is a constant in the range of f .Thus, the level curves are just the traces of the graph of f in the horizontal plane z ¼ k projected down to the xy-plane, see Figure 1(a).
Since these contours are drawn for equally spaced values of C, the spacing of the contours themselves conveys information about the relative steepness at various places on the terrain surface; the land is steepest where the contour lines are closest together.It is somewhat flatter where they are farther apart.The slope along the level curve is not ascent or descent.When the route of steepest descent (Streams) or steepest ascent can be drawn by making it perpendicular to all the contour lines see Figure 1(b).
The slope of the tangent line The slopes of tangent lines on the surface at P a; b; c ð Þ to the traces C 1 and C 2 of S in the planes y = b and x = a, as shown in Figure 2(a), are the geometric expression of the partial derivatives of z ¼ f x; y ð Þ: They represent the rate of change of z in the x and y directions.But when the rate of change of z at x 0 ; y 0 ð Þ in the direction of an arbitrary unit vector u ¼ a; b ð Þ is needed, the surface S with the equation: Then the gradient of f is the vector function defined Ñf by: In like manner, the gradient vector Ñf x 0 ; y 0 ð Þ at a point P x 0 ; y 0 ð Þ on the level curve f x; y ð Þ ¼ k gives the direction of fastest increase of f .Also, by considering the tangent plane, it can be shown that Ñf x 0 ; y 0 ð Þ is perpendicular to the level curve that passes through P.
The surface function z ¼ f x; y ð Þ has two variables, then its partial derivatives f x and f y are also function of two variables, which are called the second partial derivatives of f : The gradient vector is represented by an arrow which indicates a quantity that has both magnitude (length of the arrow) and direction (the arrow points).The dot product (scalar product) of the gradient vector with another vector (such as the tangent vector) is not a vector, it is a real number.the dot product can be given as a geometric interpretation in terms of the angle θ between two vectors, which is defined to be the angle between the representation of gradient vector and tangent vector.
The two nonzero vectors are called perpendicular or orthogonal if the angle between them is θ ¼ π=2 (Figure 3), then the gradient vector and the tangent vector are orthogonal if and only if:

Principles of contour lines densification
Delaunay triangulation is ideal for regular data, whether these data are points or lines.For irregular data, such as contour lines whose curvature and divergence vary according to the terrain pattern they represent (hilly, undulated, and plain), some rules, constraints, and enhancement tools are needed.Furthermore, in some cases manual editing, via changing the direction of some triangle's sides, is essential to improve the triangulation output.The key challenges in this research are the dispensing of manual editing of Delaunay triangulation to fix all criteria and constraints, finding an apparent continuity of data that is not continuous, and ensuring that the lengths of the triangle's sides connect the height data are as short as possible.The proposed methodology aims at mitigating the errors that arise in topographic modeling based on the densification of height data represented by contour lines.That permits well graduation of height values along with small horizontal distances.Then the structure of the produced model becomes able to re-represent the morphological features without disappearing or interruption.
The introduced method is focusing on improving the digital TIN models through the creation of additional intermediate curves between the original contour lines.These new curves can be obtained by linear interpolation along the rays connecting the original curves' nodes and perpendicular to them.Therefore, an infinite number of rays between each two adjacent contour lines can be implemented according to the densification of nodes along the curves, or a specific number of rays using only the control nodes of each contour line can be executed.These rays may be continuous, i.e. they share the same nodes on the contour line in an ascending or descending way or they are independent of each other.Four different approaches are examined to find out the appropriate method to define rays between contour lines which mathematically represent rays that are perpendicular to the tangents at each node and reach the adjacent contour line i.e. the steepest slope.These approaches are: By applying linear interpolation processes along the rays and obtaining the additional contour lines, along with the original lines they will be the input data to produce an optimized terrain surface.So, it is necessary to judge the best approach that can be adopted in the final methodology to improve the terrain modeling.

The proposed algorithm
Based on the position and direction of the steepest slope, the proposed algorithm was designed to accurately extract the elevation data in between contour lines and use them for DEM generation.The algorithm begins by checking the nodes on each contour line to establish rays that are orthogonal to the tangent of the contour at each node.The shortest ray represents the steepest slope between two adjacent contour lines and the added spot points work as a very small area contour line represented by one point.In the second step, linear interpolation along the produced rays takes place.The number of segments on each ray is dependent on the contour interval which can be varied from four meters to one meter depending on the final results.After dividing the rays into small segments, densification of contour lines should be carried out.At this step, performing Delaunay triangulation is a very easy task because of the existence of smooth graduated data between the original contour lines.When the Delaunay triangulation (TIN model) has been completed, the false flat area is estimated.If the ratio of false flat area is acceptable according to the initial value or the contour interval equals one meter, the DEM can be generated from the enhanced TIN model, otherwise, the contour interval should be decreased by one and returned to the densification step.Generation of contour lines from the produced DEM model is a very important step to compare them with the original ones.If the interpolated contours' symmetric differences area does not exceed a threshold value equal (T), the contour interval should be decreased by one and returned to the densification step.The process is repeated until obtaining a DEM that meets all predefined criteria.The flowchart of the proposed algorithm is shown in Figure 4. Civil 3D software has been used for processing.
Input: contour lines extracted from 1:25,000 topographic map in a vector format, height values of spot points (m), initial contour interval value (m), the allowed ratio of false flat area as a percentage and symmetric differences area T (m 2 ).Output: DEM 1. Assigning elevation value for peak points.2. Delineating orthogonal rays between contours.3. Liner interpolation along the orthogonal rays using initial contour interval value.4. Construction of intermediate contour lines.5. Delaunay triangulation (TIN model) based on the intermediate contours concerning the associated contour interval.6. False flat areas ratios estimation for the produced TIN models.7. DEMs generation from the produced TIN models.8. Generation of the contour lines from the DEM. 9.The process is repeated until obtaining a DEM that meets all predefined criteria.

The best approach for orthogonal vectors execution
To select the best approach for orthogonal vector execution, an experiment using a set of contour lines was implemented.The used contour lines have a 5 m contour interval, they were extracted from 1:25,000 topographic map as a part of the study area of different terrain patterns varying between hilly, undulated, and plain.First, the topographical skeleton was generated for the selected region based on the original contours without any refinements.The results indicate that the sum of false flat area reaches 15.42%.The formation of such areas is mainly due to the wrong distribution of triangles sides: either the vertices of the triangle are located in the same contour line or two contours of the same height, or due to the absence of the vertical value, i.e. discontinuity as shown in Figure 5.
Thereafter, the four abovementioned approaches were examined to select the best one to perform the orthogonal vectors.Besides the general matching and congruence between the new intermediate contours and the original ones, the main criterion to choose the best method is which method can generate the minimum flat area when modeling the terrain surface?Preforming these vectors is subjected to the following rules: 1) first, these vectors are assumed to be not intersected at all, and 2) second, the distances between two adjacent contour lines along these directions are to be as minimum as possible.In the next step, the linear interpolation along the rays was performed to find out the locus of the new intermediate contours with a vertical interval of 1 m between the original contour lines.After that, the Delaunay triangulation method was applied to generate the TIN model using the whole contour lines without enforcing restrictions or applying any modeling enhancement technique.Finally, the false flat areas were estimated, and new contour lines were extracted from the generated TIN surface to compare them with the original ones.Figure 6 presents the results of the orthogonal vectors using the four methods (CAOV, CDOV, IAOV, IDOV).Table 1 shows a comparison between the statistics for each method in terms of the false flat areas as a percentage and the symmetric differences between intermediate contour lines and the new ones generated from TIN surface.It is clearly seen that all models were improved when the false flat areas were reduced from 15.42% in the original model generated from the original contours only (5 m contour interval) to less than 8% in all the enhanced models.However, when comparing the new contour lines extracted from TIN surface with the intermediate contours produced from densification, IAOV approach has proved that it is the best method among the four approaches to enhance the TIN modeling based on the densification of contour lines.
From the above experiment, it appears that the false flat areas are concentrated in the regions of small slopes due to the large horizontal distances separating the contour lines as well as the lack of necessary information expressing the slope changes, while they rarely appear in the locations where the original contour lines are very close.However, even if they appeared, they are concentrated where the contour line takes the shape of V letter.In these locations, in particular, the importance of vertices densification along the contour lines is evident to prevent the formation of triangles on the same contour lines during TIN generation.One of the most TIN modeling enhancement procedures is to find the streamlines (lines of changing greatest slops) which pass through these points because they provide a series of points of known height which help in reducing the length of triangle sides.It is also noticed that the locations of flat areas in most models are inside closed contours.To avoid this problem, known elevation points should be inserted or calculated internal point height by linear interpolation from the neighborhood information.By applying the IAOV method, for improving the terrain modeling, minimizing the false flat areas as much as possible, and achieving the lowest differences between newly generated contours and the intermediate ones, the suitable vertical spacing must be found by iteration between intermediate contour lines and the number of vertices or nodes along each contour line used for drawing the orthogonal vectors.In addition, the top of a hill or bottom of a valley which represents the center of each closed contour line, must be identified, and be assigned a height value that equal to the elevation of the contour that contains this spot point plus or minus half of the vertical interval in the case of top or bottom, respectively.Then, the intermediate contours within that closed curve (top or bottom) can be generated parallel to the origin closed contour line and with horizontal separating distances that vary according to the point elevation and the shortest distance from the point to the contour line.When applying this procedure to the previous experiment, the false flat areas decreased to 1.34% after densifying the contours to reach 0.5 m contour interval and to 0.9% when the contour interval was 0.25 as shown in Figure 7, where original contours appeared in brown color, blue and pink are related to densified contours to a vertical interval of 0.5 m and 0.25 m, respectively, and red color is for the false flat areas.This means that densification of contour lines iteratively using the IAOV method is an effective alternative to reinforcing the model with break lines (maximum slope change lines) which represent the basic morphology of the topographic surface.

Accuracy assessment of the proposed method
To investigate the accuracy of the proposed method, three samples (hilly, undulated, and plain) were selected from the study area as shown in Figure 8 and subjected to the application of the proposed methodology.The choice of the selected regions was based on the difference of the horizontal distance between contour lines for each type which permits studying the effect of terrain variation on the final resolution of the model as well as the effect of each site on the quality of terrain modeling after densification of contours.In literature, extensive critical reviews of the positional accuracy assessment methods applied to the case of DEMs are made, and it is noteworthy that the majority approach is based on the use of points, regardless of the standard or method that is applied (Alganci, Besol, and Sertel 2018;Szypula 2019;Polidori and Hage 2020;Mesa-Mingorance and Ariza-López 2020).In this study, DEMs generated by the proposed method will be evaluated from three points of view: vertical accuracy using check points along with a comparison with TOPO raster DEM that based on ANUDEM algorithm, visual analysis of the spatial distribution of contours generated from DEMs, and the ability of DEMs to identify drainage lines.

Vertical accuracy assessment
To assess the vertical accuracy of the proposed approach, 31 check points that are not included in calculation were used.The heights of these points were extracted from the TIN models and the differences with their corresponding true values were estimated.When comparing height differences, the following statistical measures were used: Mean Absolute Error (MAE) and Root Mean Square Error (RMSE) as the following equations. .where Z i is the height calculated from the model and H i is the true value for the check point.
The proposed methodology was implemented in its four stages to examine its ability to enhance the terrain modeling and improve the production of flow lines.In the first stage, the best approach for orthogonal vector execution IAOV was applied on the original contour lines (5 m interval) in the three regions: plain, undulated and hilly terrain.Then in the second stage, linear interpolation along orthogonal rays was estimated to obtain locus of the vertices used in contour lines densification with equal spacing: 2.5 m for one additional contour line, 1.67 m for two additional contour lines, 1.25 m for three additional contour lines, and 1 m for four additional contour lines.In the third stage, digital terrain models were generated according to the densification level.They are TIN1, TIN2, TIN3,  and TIN4 respectively, where TIN0 refers to the model generated based on the original contour lines.The differences in the compared elevations with their corresponding values extracted from TIN0 ranged from −0.89 m to 0.99 m in plain location, from −1.62 m to 1.51 m in undulated location and from −1.56 m to 1.37 m in hilly location.After applying contour densification, in TIN4, these differences reach −0.66 m to 0.46 m, −0.83 m to 0.58 m, and −0.20 m to 0.42 m for plain, undulated and hilly locations, respectively.The values of the MAE and RMSE were improved in all locations, Table 2 presents the precisions indicators calculated for each TIN model, while Figure 9 shows the effect of contour line densification on the resultant accuracy.It is noted that the precision of the model increases with the increase in the number of added contours (decreasing the vertical divergence).
The GRID elevations are calculated by weighted linear interpolation based on the closest triangles' nodes of the TIN network.On the other hand, three DEMs for the selected locations are generated in ArcGIS software using TOPO raster that based on the ANUDEM algorithm to be compared with DEM4 for the selected locations.The spatial resolution of all models was 5 m.Table 3 shows the geomorphometric parameters of DEMs and their statistics in the selected locations.The most commonly used geomorphometric local parameters are altitude, local relief, aspect, slope and curvature with their statistical measures.
The first parameter was altitude which is the vertical distance expressed in meters from the reference level-surface.It is found that there is no big difference between TOPO raster DEM and DEM4 when comparing altitude values in all locations.The next parameter was local relief which is altitude range between the highest and the lowest points expressed in meters.Regardless the improvement in plain location, the results were convergent for undulated and hilly locations when comparing TOPO raster DEM and DEM4 local relief results.The next parameter is slope, which is the maximum rate of change in value from that cell to its neighbors.When comparing slope results, it is found that DEM4 is acting better than TOPO raster DEM in plain location while this improvement was not significant in the other locations.To overcome the problem that the accuracy indicators are globally reported, and no indication of the error spatial distribution is given, a differential height map between the DEM4 and TOPO raster DEM is performed.Figure 13 depicts the spatial distribution of elevational changes between both models in the three selected regions, plain, undulated and hilly.The largest elevation differences occurred in hilly location.
Generally, one can conclude that DEM based on contour densification deals well with variable, it shows advantages than TOPO raster DEM especially in undulated and hilly regions.It can describe the topography smoothly and effectively.

Visually analyzing contour lines generated from DEMs
A recreation of contour lines based on the DEM models was done as well as the position and attitude of the greatest slope lines from the starting points of the watershed.A comparison between original contours and their counterparts generated from grid models in the three sample locations is depicted in Figure 14.
The outputs were evaluated visually and statistically.It is seen that the smaller the vertical spacing between the contour lines, the more converge and match the generated curve with the original contour line.It is worth noting that this convergence is affected by the topographical pattern considered, as well as the spatial resolution of the grid model.

The ability of DEMs to identify drainage lines
Since the distribution of the test points used in the precision calculation is not sufficient to meet the qualitative criteria for the quality of the streamline extraction, a special quantitative study was adopted to improve the quality standard in calculating the probability of matching the streams derived from grid models with their original paths.
The congruence between the original trajectory and pathway extracted from the elevation model was examined using a buffer distance equal to the original contour interval (5 m) as a boundary confidence level along the path as shown in Figure 15.The reason behind using this distance is due to the abovementioned reasons which lead to obtaining branched paths of the hydrographic network in different directions, but close to the real flow path.Table 4 shows the percentage probability of convergences between extracted and original stream paths in the selected regions.
As seen in Figure 16, it is found that the probability of congruence between the derived and original stream paths increases in the three locations as the vertical divergence between contour lines decreases.However, the weakness of probability ratios for the hilly sample compared to the undulated and plain samples can be explained by the misallocation of the start point at which the flow path begins (perigee point), in addition to the multiplicity of maximum slopes directions in this terrain pattern as a result of the convergence of the intermediate contour lines after densification process, as well as the effect of horizontal spatial resolution of the grid model.

Discussion
As mentioned in Section 2, shortcomings appear frequently when using contour lines to generate DEMs especially in areas with complicated topographical changes (Ardiansyah and Yokoyama 2002).The used  techniques cannot reflect the topography of mountaintops and depressions with sufficient resolution for most purposes.By regenerating the contour lines of the terrain model to examine their congruence with the original contours, deviation and departure from the parent lines were found in most cases.This paper proposed an advanced technique based on densification of contour lines to overcome the drawback of DEM generation from contours.Contour lines of 5 m contour interval were extracted from 1:25,000 topographic map as a part of the study area of different terrain patterns varies between hilly, undulated, and plain are used to prove the validity of our methods.After describing the methods proposed and evaluating the outcomes of the experiments, it is important to outline its contribution compared to other state-of-the-art approaches.
1) The results shown in Table 1 and Figure 5 proved that IAOV approach is the best approach to perform the orthogonal vectors between contour line.The false flat areas were reduced from 15.42% in the original model generated from the original contours only to less than 8% in all the enhanced models, this highlights the fact that our suggested algorithm has considerable advantages in dealing with false flat areas.
2) Abovementioned methods for DEM generation from contour lines (Rui et al. 2016;Taud, Parrot, and Alvarez 1999;Wang, Gao, and Wu 2005) show drawback and low performance on complex topography.They started to interpolate grid points among contours regardless of the position of the interpolated point i.e. does it fall on the steepest slope or not?In comparison, the proposed approach finds the steepest   slopes first, represented by orthogonal rays, and the interpolation is made along these rays using suitable contour interval value.On the other hand, the new technique is able to interpolate inside closed contour lines after defining the peak points.
3) The results provide satisfactory elevation precision of TIN and GRID models as depicted in Figures 9 and 10, it is also found that the precisions of the gridded models are close to those of the original TIN models from which they were derived.Furthermore, noticeable enhancement was found in the convergence between original contour lines and those generated after applying the proposed methodology.This indicates the ability of the proposed method to reflect geomorphic features very closely.
4) Based on the results of the quantitative study of the convergence between extracted and original stream paths, a reasonable percentage of matching was found in the three selected locations.However, in hilly location some weakness was found in the ratio of matching which might come from the fact that in this location, it is difficult to identify the start point where the flow path begins besides the complexity to detect the slope direction in this terrain pattern.
However, it is too difficult to obtain a complete congruence between surface runoff pathways found in maps and those derived from DEMs based on grid analysis due to many reasons.On one hand, the runoff pathways seen in maps represent real-world features that have been surveyed and drawn and not lines derived from topographic models or contours.In addition, in medium scales, according to the map resolution, the pathway is subjected to cartographic generalization and appeared as a line however in fact it is a water body bed defined by two banks.In other words, the pathway has a natural terrain sanctuary of a width that reaches a few meters.Thus, these lines in maps are more of a spatial symbolization than a true representation of the morphological shape of the watershed.On the other hand, the automatic extraction of pathways based on the digital elevation models is heavily affected by many factors, including the slope analysis algorithm adopted to calculate the values of the maximum slopes and their directions, the lack of correct determination of peak and perigee points, in addition to the effect of resolution on the generalization of the elevation values of the topography of the Earth's surface.

Conclusions
Digital elevation models based solely on original contour lines will not be able to correctly extract the hydrographic network paths as the increase of associated errors and decrease of reliability due to the heterogeneity and discontinuity of the spatial data used in the terrain modeling.However, enhancing the terrain model with stream paths from external sources will not be sufficient to obtain a reliable hydrologically appropriate DEM.In the methods developed, a proposed algorithm is integrated and performed to overcome the discontinuity of modeling data through densification of contour lines.Applying contours intensification algorithm to improve the topographic modeling has shown a significant promise in extracting accurate elevation products such as grid DEMs and hydrographic network paths for determining the basin boundaries and predicting the areas of inundation and flooding.The convergence between original contour lines and the generated ones is affected by the number of intermediate contour lines added and the topographical pattern considered.A special quantitative study was adopted to examine the degree of matching between stream paths derived from grid models based on the proposed algorithm and the original stream paths.Based on the analysis of the obtained results, the algorithm has a significant advantage in enhancement terrain modeling based on contour lines especially in the undulated and plain samples.This approach can be applied when there is no high-resolution DEM for the given area, but there are topographic maps when there is a need to compare a contemporary DEM with a historic topography for an area, and to fix the errors of freely-available global DEMs.Future research may also focus on developing an algorithm to find out the location of perigee points to overcome the weakness found in modeling hilly locations.
Figure 2(b).The point P x 0 ; y 0 ; z 0 ð Þ lies on the surface S, the vertical plane passes through P in the direction of u intersect S in a curve C, and the slope of the tangent line T to C at the point P is the rate of change of z in the direction of u.Then the directional derivative of f x 0 ; y 0 ð Þ in the direction of a unit vector u ¼ a; b ð Þ is:

Figure 2 .
Figure 2. (a) Slopes of tangent lines on the surface; (b) the vertical plane passes through the direction of unit vector (James 2017).

Figure 3 .
Figure 3. Gradient vector and tangent to the level curve are orthogonal vectors (James 2017).

Figure 4 .
Figure 4. Flow chart of the proposed algorithm.

Figure 5 .
Figure 5. False flat area appeared due to the wrong distribution of triangles sides.

Figure 6 .
Figure 6.(a) The Orthogonal vectors based on four adopted methods (CAOV, CDOV, IAOV, IDOV) from left to right, (b) The intermediate contour lines extracted by linear interpolation along orthogonal vectors, (c) The new contour lines generated from TIN Models and the area in red color of false flat area, and (d) the symmetric differences between intermediate contour lines and newly generated contour lines.

Figure 7 .
Figure 7. (a) The densification of level curves 0.5 m vertical interval, causing the generation of TIN model with 1.34% of false flat area (red regions), and (b) The densification of level curves 0.25 m vertical interval reduces the false flat area to 0.9%.

Figure 8 .
Figure 8.The examination samples of the study area according to three characteristic terrain locations, scale 1:25,000.

Figure 9 .
Figure 9.The precisions of TIN models as the increase of the intermediate contour lines.

Figure 10
Figure10presents the slope map for all locations.Curvature parameter tool calculates the second derivative value of the input surface on a cell-by-cell basis, i.e. it describes the shape of a slope.The statistics in Table3along with the spatial distribution shown in Figure11prove the superiority of DEM4 than TOPO raster DEM especially in undulated and hilly regions.Figure11(a,c,e) are clear and reflect characteristic elements of the topography more than Figure11(b,d,f) that based on TOPO raster DEM.Aspect identifies the downslope direction of the maximum rate of change in value from each cell to its neighbors.values indicate the Figure10presents the slope map for all locations.Curvature parameter tool calculates the second derivative value of the input surface on a cell-by-cell basis, i.e. it describes the shape of a slope.The statistics in Table3along with the spatial distribution shown in Figure11prove the superiority of DEM4 than TOPO raster DEM especially in undulated and hilly regions.Figure11(a,c,e) are clear and reflect characteristic elements of the topography more than Figure11(b,d,f) that based on TOPO raster DEM.Aspect identifies the downslope direction of the maximum rate of change in value from each cell to its neighbors.values indicate the

Figure 13 .
Figure 13.Spatial distribution of elevational changes between DEM4 and TOPO raster DEM (a) Hilly location; (b) Undulated location; and (c) plain location.

Figure 14 .
Figure 14.Comparison between original contour lines (brown color) and the generated ones (red color) from the extracted DEMs of the three types of terrain, line of maximum slope (blue color).

Figure 15 .
Figure 15.The convergence between extracted and original stream paths.

Figure 16 .
Figure 16.The convergence probability between extracted and original stream paths as the increase of the intermediate contour lines.

Table 1 .
Comparison between the statistics of tested methods.

Table 2 .
The precisions of TIN models.

Table 3 .
The geomorphometric parameters of DEMs and their statistics.

Table 4 .
Percentage probability of convergences between extracted and original stream paths.