Heuristic and Numerical Geometrical Methods for Estimating the Elevation and Slope at Points Using Level Curves. Application for Embankments

: Both the calculation of ground slopes at points on the map and the elevation estimation for these points bear signiﬁcant importance and also have applications in various domains, such as civil engineering, road and railway design. The paper presents two methods that use level curves: one that is fast and approximate and another which is slower, but more precise. The running speed of the two proposed methods and their results are compared by performing 100 million experiments. The paper also presents how these methods can be applied to optimize embankments. An accurate method to calculate the horizontal plane of the excavation/ﬁlling when building a new house is also presented.


Introduction
The possibility of calculating the slope and elevation using level curves has practical applicability in several fields of activity, such as civil engineering [1], architecture [2,3], railway [4], road and bridge engineering [5][6][7][8], topography [9], etc. Although the results found can be used in various areas, this article is limited to solving two problems related to civil engineering and to railways, roads, and bridges, both ultimately related to the embankments optimization. The problem is well known, but it was not closely examined by design engineers, who preferred an intuitive, approximate solution. The current article presents a rigorous and mathematically accurate solution to find the optimum between excavations and fillings, resulting in a reduction of the embankments related costs. The article presents in detail the solutions of two concrete problems in two different fields of activity.
The first problem investigated is in the field of civil engineering and consists of determining the optimal elevation corresponding to the finished floor of the ground level, when the building has to be positioned on a slope [1][2][3]. The second problem investigated is in the design of a new railway or road when there are several fixed points on the topographic map through which the new road [10][11][12][13] or the new railway [4,11] will pass.
The designer has to know the slope and elevation in each of these points as accurately as possible. Using the methods proposed in this paper, the two parameters (slope and elevation) can be determined fast and with a high level of accuracy using the level curves of the map, without the need for direct measurement on the spot. Furthermore, when designing roads or railways, the risk of failure, the frictional resistance, the safety factor against sliding or failure, the translational fill failure, or the cut volume and the road fill per meter require knowing the slopes (side slope, cut slope, or fill slope).
The paper is organized as follows: an algorithm to calculate the ground slope at a given point, using level curves, is deduced in Section 2. This problem is reduced to Since a level curve is approximated with consecutive coplanar segments of straight lines, our problem can be reduced to finding the shortest length segment |P1P2| that connects two coplanar line segments located on two consecutive level curves. When the two 2D points P1(x1,y1) and P2(x2,y2) are determined, they are transformed into 3D points by adding the elevations z1 and, respectively, z2 of the level curves on which each of the points are located. Using the 3D coordinates of these 2 points, the slope m = tan(α) at point Q (see Figure 2) can be calculated as follows: Since a level curve is approximated with consecutive coplanar segments of straight lines, our problem can be reduced to finding the shortest length segment |P 1 P 2 | that connects two coplanar line segments located on two consecutive level curves. When the two 2D points P 1 (x 1 ,y 1 ) and P 2 (x 2 ,y 2 ) are determined, they are transformed into 3D points by adding the elevations z 1 and, respectively, z 2 of the level curves on which each of the points are located. Using the 3D coordinates of these 2 points, the slope m = tan(α) at point Q (see Figure 2) can be calculated as follows: The first problem (denoted IPS) intended to be solved is to identify the pair of segments (located on two level curves) to which the finding of the shortest length segment should be applied. Of course, for each of the two level curves, the candidates are the closest segments to point Q. Let us denote by Si the set of candidates (segments) for curve i (i = 1, 2). For each pair of segments s1 ∈ S1 ands2 ∈ S2, respectively, the corresponding shortest length segment [P1P2] between s1 and s2 is calculated. If s1 contains P1 and s2 contains P2, then the pair (s1, s2) is added to the set (denoted C) of solution candidates for IPS ( Figure 3). The set Si (i = 1, 2) can be built as follows: first, we find the point T belonging to the level curve i that is closest to Q. We consider then k ≥ 1 points (L1, L2, …, Lk) to the left of T on the level curve and k points (R1, R2, …, Rk) located to right (see Figure 4). Finally, the set Si consists of 2k segments of straight lines given by the consecutive points on the level curve i: Lk, Lk-1, …, L1, T, R1, R2, …, Rk, i.e.,: Now, let us see how k (giving the number of points in the vicinity of T) is considered. Usually (in real tests), k is enough to be 1. In few cases k is found to be greater. So, the calculations are started with k = 1 and, if necessary, it is increased.
The Algorithm 1 for solving IPS (denoted AIPS) is: Determining the slope at Q using P 1 and P 2 .
The first problem (denoted IPS) intended to be solved is to identify the pair of segments (located on two level curves) to which the finding of the shortest length segment should be applied. Of course, for each of the two level curves, the candidates are the closest segments to point Q. Let us denote by S i the set of candidates (segments) for curve i (i = 1, 2). For each pair of segments s 1 ∈ S 1 and s 2 ∈ S 2 , respectively, the corresponding shortest length segment [P 1 P 2 ] between s 1 and s 2 is calculated. If s 1 contains P 1 and s 2 contains P 2 , then the pair (s 1 , s 2 ) is added to the set (denoted C) of solution candidates for IPS ( Figure 3).  The first problem (denoted IPS) intended to be solved is to identify the pair of segments (located on two level curves) to which the finding of the shortest length segment should be applied. Of course, for each of the two level curves, the candidates are the closest segments to point Q. Let us denote by Si the set of candidates (segments) for curve i (i = 1, 2). For each pair of segments s1 ∈ S1 ands2 ∈ S2, respectively, the corresponding shortest length segment [P1P2] between s1 and s2 is calculated. If s1 contains P1 and s2 contains P2, then the pair (s1, s2) is added to the set (denoted C) of solution candidates for IPS ( Figure 3). The set Si (i = 1, 2) can be built as follows: first, we find the point T belonging to the level curve i that is closest to Q. We consider then k ≥ 1 points (L1, L2, …, Lk) to the left of T on the level curve and k points (R1, R2, …, Rk) located to right (see Figure 4). Finally, the set Si consists of 2k segments of straight lines given by the consecutive points on the level curve i: Lk, Lk-1, …, L1, T, R1, R2, …, Rk, i.e.,: Now, let us see how k (giving the number of points in the vicinity of T) is considered. Usually (in real tests), k is enough to be 1. In few cases k is found to be greater. So, the calculations are started with k = 1 and, if necessary, it is increased.
The Algorithm 1 for solving IPS (denoted AIPS) is: The set S i (i = 1, 2) can be built as follows: first, we find the point T belonging to the level curve i that is closest to Q. We consider then k ≥ 1 points (L 1 , L 2 , . . . , L k ) to the left of T on the level curve and k points (R 1 , R 2 , . . . , R k ) located to right (see Figure 4). Finally, the set S i consists of 2k segments of straight lines given by the consecutive points on the level curve i: L k , L k-1 , . . . , L 1 , T, R 1 , R 2 , . . . , R k , i.e.,: Now, let us see how k (giving the number of points in the vicinity of T) is considered. Usually (in real tests), k is enough to be 1. In few cases k is found to be greater. So, the calculations are started with k = 1 and, if necessary, it is increased.
The Algorithm 1 for solving IPS (denoted AIPS) is: The shortest line segment that connects both level curves is the shortest segment [P 1 P 2 ] from C.
After applying AIPS, the problem is reduced to calculating the shortest segment that passes through a given point Q and connects two coplanar lines d 1 and d 2 ( Figure 4). We shall denote this problem as SSTPC2L (shortest segment through a point that connects 2 lines). In Figure 4 the points A 11 and A 12 determine the line d 1 , and the points A 21 and A 22 determine the line d 2 .
Algorithm 1 Solving IPS (denoted AIPS) AIPS: Input: the points of the closest two level curves to Q Output: the shortest line segment [P1P2] that connects the given two level curves finished = false; While not finished do C = ϕ (empty set); Build the sets of line segments S1 and S2 (see Equation (2)); For each segment s1 from S1 do For each segment s2 from S2 do Calculate the shortest length segment P1P2 that connects s1 and s2 and passes through Q; If P1 ∈ s1 and P2 ∈ s2 then Add [P1P2] to C; End if; End for; End for; If C not is empty then finished = true; else k = k+1; End if; End while; The shortest line segment that connects both level curves is the shortest segmen [P1P2] from C.
After applying AIPS, the problem is reduced to calculating the shortest segment tha passes through a given point Q and connects two coplanar lines d1 and d2 ( Figure 4). We shall denote this problem as SSTPC2L (shortest segment through a point that connects 2 lines). In Figure 4 the points A11 and A12 determine the line d1, and the points A21 and A2 determine the line d2. We present two methods of solving SSTPC2L: an approximate method and an exac mathematical method. These two methods (presented in Sections 2.1 and 2.2) were partially described in [7]. We present two methods of solving SSTPC2L: an approximate method and an exact mathematical method. These two methods (presented in Sections 2.1 and 2.2) were partially described in [7].

Heuristic (Approximate) Method for SSTPC2L
From the topographic map, the point Q is known and two segments A 11 A 12 and A 21 A 22 are determined belonging to the level curves from Figure 1  and A 21 A 22 on curve 620). These two segments are found according to Figure 3. The method consists of drawing the perpendiculars from the given point Q on the two given lines d 1 and, respectively, d 2 ( Figure 5). We construct a line passing through Q which is parallel to the line determined by the feet of the two perpendiculars. The obtained line is denoted by d 3 , and its intersection with the lines d 1 and d 2 is the searched segment [P 1 P 2 ] (approximate solution of SSTPC2L).
Appl. Sci. 2021, 11, x FOR PEER REVIEW 5 of 22 From the topographic map, the point Q is known and two segments A11A12 and A21A2 are determined belonging to the level curves from Figure 1 (A11A12 on curve 610 and A21A2 on curve 620). These two segments are found according to Figure 3. The method consists of drawing the perpendiculars from the given point Q on the two given lines d1 and respectively, d2 ( Figure 5). We construct a line passing through Q which is parallel to the line determined by the feet of the two perpendiculars. The obtained line is denoted by d3 and its intersection with the lines d1 and d2 is the searched segment [P1P2] (approximate solution of SSTPC2L). The perpendicular from Q on di (i = 1, 2) is: The foot Q (x , y ) (i = 1, 2) of the perpendicular on di is obtained by solving the system of Equations (3) and (4).
The slope of the line Q1Q2 is: Since m is the slope of Q Q , then the parallel through the point Q(x , y ) is: Finally, the points Pi (i = 1, 2) are obtained by solving the system of Equations (3) and (6).
The following Algorithm 2 (denoted A1SSTPC2L) is obtained to solve the problem SSTPC2L: Algorithm 2 Solving the problem SSTPC2L A1SSTPC2L: Input: Q and the points of the closest two level curves to Q Output: The slope at point Q Apply AIPS to find A11, A12, A21, A22; Now, let's present how the points P 1 and P 2 are calculated. The equation of the line d i (i = 1, 2) starting from the points A i1 and A i2 is: The perpendicular from Q on d i (i = 1, 2) is: The foot Q i x Q i , y Q i (i = 1, 2) of the perpendicular on d i is obtained by solving the system of Equations (3) and (4).
The slope of the line Q 1 Q 2 is: Since m 12 is the slope of Q 1 Q 2 , then the parallel through the point Q x Q , y Q is: Finally, the points P i (i = 1, 2) are obtained by solving the system of Equations (3) and (6).
The following Algorithm 2 (denoted A1SSTPC2L) is obtained to solve the problem SSTPC2L: Algorithm 2 Solving the problem SSTPC2L A1SSTPC2L: Input: Q and the points of the closest two level curves to Q Output: The slope at point Q Apply AIPS to find A 11 , A 12 , A 21 , A 22 ; Calculate the point P i by solving system of Equations (3) and (6) (i = 1, 2). Calculate the slope at point Q using Equation (1).

The Exact (Mathematical) Method for SSTPC2L
For this method, in order to simplify calculus, we apply some initial transformations (rotations and translations) to the points A 11 , A 12 , A 21 , A 22 , and Q so that the bisector of the angle between the lines d 1 and d 2 becomes parallel to Oy axis and the point Q is in the origin O, i.e., Q = O (0,0). To do so, we first calculate the intersection I of d 1 and d 2 : by solving the system of linear equations: The following transformations to A 11 , A 12 , A 21 , A 22 and Q are applied: Rotation with sin(−u) = −sin(u) and cos(−u) = cos(u), where: 11 )/r sin(u) = (y 12 − y 11 )/r (8) After this rotation, the line d 1 is over Ox axis of the map's system of coordinates.

3.
Rotation with sin(v) and cos(v), where: After this rotation, the bisector of the angle given by lines d 1 and d 2 is over Oy.

Translation with
The following 2 systems of equations have to be solved in order to find the points P i (x i ,y i ) (i = 1,2) where the lines d i (i = 1,2) intersect the line d: Since A i1 , A i2 define the line d i (i = 1,2) in (10) we have: Solving the 2 systems of equations, the points P i (x i ,y i ) (i = 1,2) are found, where: The distance between P 1 and P 2 must be minimized, i.e.,: min a f(a) (15) where: By replacing (13) and (14) in (16) the following formula for f(a) is obtained: We make the following notations: It is easy to see that the function f is decreasing for a ≤ a L , and is increasing for a ≥ a R , since f is continuous, a L and a R are the slopes of the perpendiculars QQ 1 and QQ 2 from Q on the lines d 1 and, respectively, d 2 , and "a" is the slope of the line d 3 . Consequently, there is a minimum of the function f denoted a min in the interval [a L , a R ]. So, the value a min is one of the solutions of the equation: on the interval [a L , a R ], where: Remark 1. Since the bisector of the angle between d 1 and d 2 is parallel to Oy axis, it is not difficult to see that: In order to obtain the solution a min of the Equation (20), a numerical method [14] can be used such as bisection method [15], or tangent method. Since f (a L ) < 0 and After a min is calculated, the points P i (i = 1, 2) are obtained using (13) and (14) and the distance between these two points is: In order to go back to the initial system of coordinates, the inverse initial transformations must be applied in reverse order to the points P 1 and P 2 , i.e.,:
Rotation with sin(-v) = sin(v) and cos(-v) = cos(v) (calculated in (9) Find the solution a min of the equation f (a) = 0 using bisection method on the interval (a L , a R ) (see Equations (19) and (20)); Calculate the coordinates of the points P i (i = 1,2) using (13) and (14), where a = a min ; Apply the inverse initial transformations in reverse order to the points P 1 and P 2 ; Calculate the slope at point Q using (1).
We implemented the above algorithm in Visual C++ and in the bisection method we set the error err = 0.0001 for finding the value a min . In Figure 6, a graphical output of our program is presented.
Find the solution amin of the equation f′(a) = 0 using bisection method on the interval (a , a ) (see Equations (19) and (20)); Calculate the coordinates of the points Pi (i = 1,2) using (13) and (14), where a = amin; Apply the inverse initial transformations in reverse order to the points P1 and P2; Calculate the slope at point Q using (1).
We implemented the above algorithm in Visual C++ and in the bisection method we set the error err = 0.0001 for finding the value amin. In Figure 6, a graphical output of our program is presented.

Results
Besides estimation of the ground slope for a given point, other two possible problems can be solved since we have now obtained a method to calculate the points P1, and P2 on the closest two level curves to a given point Q so that P1, Q, and P2 are collinear and the distance between P1 and P2 is minimum.

Elevation Estimation at a Point
The challenge is to estimate as well as possible the elevation z of a point Q. Assuming that z1 < z2 (see Figure 3) and since QQ' is parallel to P2P'2, in the triangle P1P2P'2 we have: From (25) we obtain the elevation of Q: It is immediate that, without the assumption of z1 < z2, the following formula can be used to estimate the elevation of Q:

Results
Besides estimation of the ground slope for a given point, other two possible problems can be solved since we have now obtained a method to calculate the points P 1 , and P 2 on the closest two level curves to a given point Q so that P 1 , Q, and P 2 are collinear and the distance between P 1 and P 2 is minimum.

Elevation Estimation at a Point
The challenge is to estimate as well as possible the elevation z of a point Q. Assuming that z 1 < z 2 (see Figure 3) and since QQ' is parallel to P 2 P 2 , in the triangle P 1 P 2 P 2 we have: From (25) we obtain the elevation of Q: It is immediate that, without the assumption of z 1 < z 2 , the following formula can be used to estimate the elevation of Q:

The Slope between Two Points
The problem of estimating the slope between two points is significant, seeing that when designing roads, the slope between any two points on a road cannot exceed a given maximum slope, e.g., for highways this value has to be less than 6% or 7%. Therefore, the chosen points from the future road can be tested for this eligibility. To do that, we consider the elevations of the points Q 1 and Q 2 calculated using formula (26) and we can estimate the slope m(Q 1 ,Q 2 ) of the road passing through Q 1 (x 1 ,y 1 ,z 1 ) and Q 2 (x 2 ,y 2 ,z 2 ) using the following formula (see Figure 7): The problem of estimating the slope between two points is significant, seeing that when designing roads, the slope between any two points on a road cannot exceed a given maximum slope, e.g., for highways this value has to be less than 6% or 7%. Therefore, the chosen points from the future road can be tested for this eligibility. To do that, we consider the elevations of the points Q1 and Q2 calculated using formula (26) and we can estimate the slope m(Q1,Q2) of the road passing through Q1(x1,y1,z1) and Q2(x2,y2,z2) using the following formula (see Figure 7):

Optimum Excavation/Filling
When a house is built on sloping land and/or with elevations and depressions, the elevation of a horizontal plane must be calculated so that the excavations over the plane would fill the space below the plane resulting in a horizontal platform on which a house can be built. In such a manner, the optimal solution between excavations and fillings is calculated, resulting in a reduction in the costs of embankments. We shall follow up with a presentation of a method to accurately calculate the elevation of this plane. The wellknown "flood fill" algorithm [16] from computer graphics is adapted to deal with this problem. Since the algorithm runs in a discrete space, the coordinates will be transformed into their discrete counterparts.
We start with the level curves from the vicinity of the house (Figure 8). Then a discretization is applied (Figure 9), meaning that the plane is divided into equidistant 2D points, e.g., in Figure 9, the distance of 40 cm between points was considered. The points

Optimum Excavation/Filling
When a house is built on sloping land and/or with elevations and depressions, the elevation of a horizontal plane must be calculated so that the excavations over the plane would fill the space below the plane resulting in a horizontal platform on which a house can be built. In such a manner, the optimal solution between excavations and fillings is calculated, resulting in a reduction in the costs of embankments. We shall follow up with a presentation of a method to accurately calculate the elevation of this plane. The well-known "flood fill" algorithm [16] from computer graphics is adapted to deal with this problem. Since the algorithm runs in a discrete space, the coordinates will be transformed into their discrete counterparts.
We start with the level curves from the vicinity of the house (Figure 8). Then a discretization is applied (Figure 9), meaning that the plane is divided into equidistant 2D points, e.g., in Figure 9, the distance of 40 cm between points was considered. The points of the contour of the house are transformed into their discrete counterparts (for each such a point, the closest discrete point is determined). Using Bresenham's algorithm [17], the discrete points of the contour are found (green points in Figure 9). Using a flood fill algorithm [16] in the discrete plane of the house, the discrete points inside the contour are determined (red points from Figure 9). The elevation of each red or green point is calculated using the level curves. Using these elevations, the optimal elevation of the excavation/filling plane is computed. So, the following algorithm denoted AOEF is obtained.   Using the elevations zi, calculate the elevation of the optimal excavation/filling plane (see AHEFP).
Applying the elevations zi of the discrete points Qi inside the contour of the house, the elevation of the optimal excavation/filling plane can be determined as follows. The elevations are sorted ascendingly. A horizontal plane π(h) (h is the height/elevation of the plane) is placed consecutively starting with the first (the lowest) elevation, continuing with the second, and ending with the last one (the highest). The following algorithm is obtained: Algorithm 5 Height of Excavation/Filling Plane (AHEFP):  Find the set of points C i using Bresenham's algorithm from point P i to point P i+1 ; end for; Using Bresenham's algorithm from point P n to point P 1 , find the set of points C n ; Find the points Q i (i = 1, 2, . . . , m) inside the contour given by C 1 ∪ C 1 ∪ . . . ∪ C n ; for each point Q i do Apply A2SSTPC2L to find the elevation z i of Q i ; end for; Using the elevations z i , calculate the elevation of the optimal excavation/filling plane (see AHEFP).
Applying the elevations z i of the discrete points Q i inside the contour of the house, the elevation of the optimal excavation/filling plane can be determined as follows. The elevations are sorted ascendingly. A horizontal plane π(h) (h is the height/elevation of the plane) is placed consecutively starting with the first (the lowest) elevation, continuing with the second, and ending with the last one (the highest). The following algorithm is obtained: Algorithm 5 Height of Excavation/Filling Plane (AHEFP):

Algorithm 5 Height of Excavation/Filling Plane (AHEFP)
Input: vector of elevations z = (z i ) i = 1, 2, . . . , m ; Sort ascending the vector z; S = 0; for i = 2 to n do S = S + z i − z 1 ; end for; S 1 = 0; S 2 = S; for i = 2 to n − 1 do At each iteration of the last "for" loop from AFHEFP, the height of the plane is considered equal to z i , S 1 is the sum of the distances between the elevation of z i the plane and the elevations below the plane, and S 2 is the sum of the distances between the elevations over the plane and the elevation z i of the plane, i.e.,: when starting a new iteration, S 1 is the sum of i-1 components of z and was calculated in the previous iteration. The new sum S 1 can be obtained from the previous by adding (i − 1) · (z i − z i-1 ). At the beginning of each iteration, S 2 is the sum of n-i+1 components of z and was calculated in the previous iteration. Thus, using the value of the previous S 2 , the new sum can be obtained from the previous by subtracting (n − i + 1) · (z i − z i-1 ).
Since S 1 is increasing from 0 to S and S 2 is decreasing from S to 0, the optimum is reached when S 2 becomes less or equal to S 1 . More exactly, the optimum height of the plane is inside the interval [z opt-1 , z opt ]. If the values z opt-1 and z opt are close enough, i.e., z opt − z opt-1 < ε, where ε > 0 is the fixed maximum distance to the optimum, e.g., ε = 5 cm, then any of the values z opt-1 and z opt are good approximations of the optimum. If z opt − z opt-1 ≥ ε, then a divide and conquer strategy is applied to get closer to optimum. To do that, two initial planes π(z opt-1 ) and π(z opt ) are considered. A new plane is placed in the middle. If the distance between S 2 calculated for π((z opt-1 +z opt )/2) and S 1 for π(z opt-1 ) is less than the distance between S 2 for π(z opt ) and S 1 for π((z opt-1 +z opt )/2), then the optimum is further calculated in the interval [z opt-1 , (z opt-1 +z opt )/2]. Otherwise, the optimum is calculated in the interval [(z opt-1 +z opt )/2, z opt ] and so on.
Algorithm 6 divide and conquer for optimum excavation/filling plane (ADCOEFP): Algorithm 6 Divide and conquer for optimum excavation/filling plane (ADCOEFP) a = z opt-1 ; b = z opt ; S a,1 = S 1 − (opt − 1) · (z opt − z opt-1 ); S a,2 = S 2 + (n -opt + 1) · (z opt − z opt-1 ); S b,1 = S 1 ; S b,2 = S 2 ; while b − a ≥ ε do c = (a+b)/2; S c,1 = S a,1 + (opt − 1) · (c − a); S c,2 = S a,2 − (n − opt + 1) · (c − a); if |S c,2 − S a,1 | < |S b,2 − S c,1 | then a = c; S a,1 = S c,1 ; else b = c; S b,2 = S c,2 ; end if; end while; S 1 = S a,1 ; S 2 = S a,2 ; At the end of the while loop, any of the values, a or b are good approximations of the optimum. We considered a as the optimal solution. The approximate excavation volume is V e = S 2 · d 2 , and the filling volume is V f = S 1 · d 2 . Consequently, the difference between the two volumes is: So, the two volumes are very similar. Figure 10 shows a real topographic map for a plot of land on which a house is to be built. The topographic map of the land was drawn up by a topometric engineer, with specific high-precision equipment, called in specialized terms "station." For small areas of land, the station is successively placed and GPS 3D coordinates are measured relative to a fixed position (in our situation it is a terminal in the Black Sea area located at a distance of approximately 717 km). The first chosen points of the topographic survey are those that describe the contour of the plot of land. In our case, the terrain contour is described by the 15 points (numbered in black). These points are inventoried and noted in a table (coordinate inventory) located on the topographic map. The following points will describe the contour of the level elevations, the positioning of the access road, the position of the neighboring houses (they are marked on the map in red together with the elevation), etc. In the case of lands with a high slope, such as the one exemplified, the measured points must be thickened for a greater accuracy of tracing the level curves (in total we have over 100 such points).

Numerical Example
The designer receives in electronic format a topographic map similar to the one in Figure 10. On this topographic support, the designer will place the designed building on the scale, having to respect a series of rigors imposed by law, among which: the minimum distance from the property limit and/or from the road axis as in Figure 11. Once the construction is located, the new contour will appear on the topographic map described by a few new points (minimum 4, in the example, for the description of the contour, 6 points are used). A situation plan is drawn up in which, according to the legislation, the 2D coordinates of the corners of the construction must be specified (see Table 1) together with the elevation of the horizontal excavation/filling plane. In Figure 11, this elevation is denoted by CTA and is estimated at 707 m by the designer. A good determination of this elevation is not easy, especially if the terrain has a high slope. The designer receives in electronic format a topographic map similar to the one in Figure 10. On this topographic support, the designer will place the designed building on the scale, having to respect a series of rigors imposed by law, among which: the minimum distance from the property limit and/or from the road axis as in Figure 11. Once the construction is located, the new contour will appear on the topographic map described by a few new points (minimum 4, in the example, for the description of the contour, 6 points are used). A situation plan is drawn up in which, according to the legislation, the 2D coordinates of the corners of the construction must be specified (see Table 1) together with the elevation of the horizontal excavation/filling plane. In Figure 11, this elevation is denoted by CTA and is estimated at 707 m by the designer. A good determination of this elevation is not easy, especially if the terrain has a high slope.  A detail of Figure 11 is presented in Figure 12 illustrating the points A 2 11 , A 2 12 , A 2 21 , and A 2 22 found using the algorithm AIPS. Table 2 presents the coordinates of the points A k 11 , A k 12 , A k 21 , and A k 22 calculated for each point Q k (k = 1, 2, . . . , 6) by applying AIPS. Using A2SSTPC2L and (26), the elevation for each point Q k (k = 1, 2, . . . , 6) is calculated.
Usually, for the dimensions of the project to be in accordance with the real construction, before commencing all the construction works, a plot is made on the field of the future building. For this purpose, the ground tracing is performed. The corners of the house are very precisely marked (the Q points). In our example, the provided measured elevations of the 6 points are presented in Table 3. Appl. Sci. 2021, 11, x FOR PEER REVIEW 15 of 22 Figure 11. Situation plan.  Figure 11 is presented in Figure 12 illustrating the points A , A , A , and A found using the algorithm AIPS.    As seen in Table 3, the estimated elevation and the measured one for each Q point are very close. The difference between them is less than 2 cm.
Using the coordinates of the points Q k (k = 1, 2, . . . , 6) from Table 1, Bresenham, and, then flood fill algorithms were applied resulting in 803 points located inside the perimeter of the house. The elevations of these points were calculated using A2SSTPC2L and (26). The algorithm AHEFP was applied for these 803 elevations and the height of the horizontal excavation/filling plane was accurately estimated as 706.2571 m, instead of the value 707 m proposed by the designer. Out of the 803 points, 412 were located below the horizontal plane, and 391 had the elevation over this plane.

Discussion
In order to compare the two proposed methods (from Sections 2.1 and 2.2) for determining the points P 1 and P 2 100 million experiments were performed. For each experiment, two random lines d 1 and d 2 were selected and the point Q was also randomly chosen.
The tests were performed using an ASUS ROG GL752VW-T4015D laptop with Intel ® Core™ i7-6700HQ 2.60 GHz processor, and 8 GB of RAM.
First, we analyzed the speed of the proposed methods. The total running time of the first method was 5.82 s. The total running time of the second method was 40.14 s. So, the first method is almost 7 times faster than the second one (see Figure 13). The first method is considerably faster because it requires elementary calculations, while the second one runs an iterative algorithm to solve the Equation (18). We also calculated the accuracy of the first method against the second method since the error for the second method can be set as low as desired. An average error of 0.61% was obtained, the lowest error being 0 and the highest 1.17% (Table 4). So, we concluded that the accuracy of the approximate first method is good enough.

Errors (First Method against Second Method)
Lowest error 0% Average error 0.61% Highest error 1.17% The steps of positioning the house will be next discussed, and the utility of the methods proposed in this paper will be shown. The designer receives a topographic map in electronic format. On this topographic support, the designer will place the designed building on the scale, having to comply with a series of rigors imposed by law, among which: the minimum distance from the property limit and/or from the road axis, etc. Once the construction is located, the outline of the house described by several new points will appear on the topographic map (minimum 4). At these new points, it is necessary to specify: the coordinates on the three dimensions according to which the elevation is specified, and the distances of the construction location from the property limits. The determination of the coordinates in the plan does not raise any concerns, but there are difficulties in obtaining the real vertical elevation especially if the terrain has a high slope. The article proposes a method for determining the exact elevation of the land at the corners of the building, and in fact for the entire contour. Depending on the contour dimensions, the designer establishes a position of the horizontal plane, a plane that corresponds to the finished floor of the ground level. The decision to establish this quota should be contingent on the volumes of embankments that will be made. The optimum is obtained when the excavation volume is approximately equal to the filling volume. The article solves this problem in Section 3.3.
For the design of a road, a primary route of the road is made on a topographic map and road recognition is performed. A final design of the road follows both along its entire We also calculated the accuracy of the first method against the second method since the error for the second method can be set as low as desired. An average error of 0.61% was obtained, the lowest error being 0 and the highest 1.17% (Table 4). So, we concluded that the accuracy of the approximate first method is good enough. Table 4. Errors of the first method obtained from 100 million experiments.

Errors (First Method against Second Method)
Lowest error 0% Average error 0.61% Highest error 1.17% The steps of positioning the house will be next discussed, and the utility of the methods proposed in this paper will be shown. The designer receives a topographic map in electronic format. On this topographic support, the designer will place the designed building on the scale, having to comply with a series of rigors imposed by law, among which: the minimum distance from the property limit and/or from the road axis, etc. Once the construction is located, the outline of the house described by several new points will appear on the topographic map (minimum 4). At these new points, it is necessary to specify: the coordinates on the three dimensions according to which the elevation is specified, and the distances of the construction location from the property limits. The determination of the coordinates in the plan does not raise any concerns, but there are difficulties in obtaining the real vertical elevation especially if the terrain has a high slope. The article proposes a method for determining the exact elevation of the land at the corners of the building, and in fact for the entire contour. Depending on the contour dimensions, the designer establishes a position of the horizontal plane, a plane that corresponds to the finished floor of the ground level. The decision to establish this quota should be contingent on the volumes of embankments that will be made. The optimum is obtained when the excavation volume is approximately equal to the filling volume. The article solves this problem in Section 3.3.
For the design of a road, a primary route of the road is made on a topographic map and road recognition is performed. A final design of the road follows both along its entire length and on its cross-sections. Figure 14 shows with a red line the optimal trajectory of a road depending on the elevation of the natural land. Because the considered road is a highway, the rules impose stricter restrictions on the slopes. If the allure of the natural terrain in certain areas is followed, the maximum allowed inclination of the road will inevitably be exceeded. As it can be seen in Figure 14, in these areas, measures must be taken to straighten the slope, such as fillings, excavations or the provision of bridges and viaducts, tunnels respectively. Through detailed awareness of the elevations and the inclination of the land in the points near the road, an optimal vertical tracing of the road can be established so that we have minimum excavations and fillings for an imposed road slope, or we can optimize the lengths of the bridges, respectively tunnels [12].
Appl. Sci. 2021, 11, x FOR PEER REVIEW 19 of 22 inevitably be exceeded. As it can be seen in Figure 14, in these areas, measures must be taken to straighten the slope, such as fillings, excavations or the provision of bridges and viaducts, tunnels respectively. Through detailed awareness of the elevations and the inclination of the land in the points near the road, an optimal vertical tracing of the road can be established so that we have minimum excavations and fillings for an imposed road slope, or we can optimize the lengths of the bridges, respectively tunnels [12]. When designing a new road or railway, a number of initial design factors must be taken into account, such as: road category, design speed, width of a traffic lane, number of lanes, traffic frequency, equipment size, etc. Depending on the route imposed on the road, there are also other specific aspects that must be taken into account: 1. Avoid high erosion hazard sites, particularly where mass failure is a possibility. 2. Utilize natural terrain features such as stable benches, ridgetops, and low gradient slopes to minimize the area of road disturbance. 3. If necessary, include short road segments with steeper gradients to avoid problem areas or to utilize natural terrain features. 4. Avoid midslope locations on long, steep, or unstable slopes. 5. Locate roads on well-drained soils and rock formations which dip into slopes rather than areas characterized by seeps, highly plastic clays, concave slopes hummocky topography, cracked soil and rock strata dipping parallel to the slope. 6. For logging road, utilize natural log landing areas (flatter, benched, well-drained land) to reduce soil disturbance associated with log landings and skid roads. 7. Avoid undercutting unstable, moist toe slopes when locating roads in or near a valley bottom. 8. Roll or vary road grades where possible to dissipate flow in road drainage ditches and culverts and to reduce surface erosion. 9. Select drainage crossings to minimize channel disturbance during construction and to minimize approach cuts and fills. 10. Locate roads far enough above streams to provide an adequate buffer, or provide structure or objects to intercept sediment moving down slope below the road. 11. If an unstable area such as a headwall must be crossed, consider end hauling excavated material rather than using sidecast methods. Avoid deep fills and compact all fills to accepted engineering standards. Design for close culvert and cross drain spacing to effectively remove water from ditches and provide for adequate energy dissipators below culvert outlets. Horizontal drains or interceptor drains may be When designing a new road or railway, a number of initial design factors must be taken into account, such as: road category, design speed, width of a traffic lane, number of lanes, traffic frequency, equipment size, etc. Depending on the route imposed on the road, there are also other specific aspects that must be taken into account:

1.
Avoid high erosion hazard sites, particularly where mass failure is a possibility.

2.
Utilize natural terrain features such as stable benches, ridgetops, and low gradient slopes to minimize the area of road disturbance.

3.
If necessary, include short road segments with steeper gradients to avoid problem areas or to utilize natural terrain features.

5.
Locate roads on well-drained soils and rock formations which dip into slopes rather than areas characterized by seeps, highly plastic clays, concave slopes hummocky topography, cracked soil and rock strata dipping parallel to the slope. 6.
For logging road, utilize natural log landing areas (flatter, benched, well-drained land) to reduce soil disturbance associated with log landings and skid roads. 7.
Avoid undercutting unstable, moist toe slopes when locating roads in or near a valley bottom. 8.
Roll or vary road grades where possible to dissipate flow in road drainage ditches and culverts and to reduce surface erosion. 9.
Select drainage crossings to minimize channel disturbance during construction and to minimize approach cuts and fills. 10. Locate roads far enough above streams to provide an adequate buffer, or provide structure or objects to intercept sediment moving down slope below the road.
11. If an unstable area such as a headwall must be crossed, consider end hauling excavated material rather than using sidecast methods. Avoid deep fills and compact all fills to accepted engineering standards. Design for close culvert and cross drain spacing to effectively remove water from ditches and provide for adequate energy dissipators below culvert outlets. Horizontal drains or interceptor drains may be necessary to drain excess groundwater [13].
As seen in Figure 15, the angle of inclination of the terrain is an important parameter to know in order to avoid erosion and then landslides [18,19]. From the point of view of the cross-sections, the position of the road in relation to the terrain slope must take into account the above constraints. By the same token, knowing the slope of the land is essential for establishing the position both horizontally and especially vertically of the road, and for achieving a low cost related to the embankment works. This position is established according to the Full Bench Road Prism, which contains the elements in Figure 16. For very steep road areas such as the one in Figure 17, knowing the slope is imperative because substantial savings in embankment volumes can be made, resulting From the point of view of the cross-sections, the position of the road in relation to the terrain slope must take into account the above constraints. By the same token, knowing the slope of the land is essential for establishing the position both horizontally and especially vertically of the road, and for achieving a low cost related to the embankment works. This position is established according to the Full Bench Road Prism, which contains the elements in Figure 16. From the point of view of the cross-sections, the position of the road in relation to the terrain slope must take into account the above constraints. By the same token, knowing the slope of the land is essential for establishing the position both horizontally and especially vertically of the road, and for achieving a low cost related to the embankment works. This position is established according to the Full Bench Road Prism, which contains the elements in Figure 16. For very steep road areas such as the one in Figure 17, knowing the slope is imperative because substantial savings in embankment volumes can be made, resulting in decreasing the quantities of gabions or excavations [13]. For very steep road areas such as the one in Figure 17, knowing the slope is imperative because substantial savings in embankment volumes can be made, resulting in decreasing the quantities of gabions or excavations [13]. The problem related to determining the slope of natural land at a point on the topographic map is necessary both in the longitudinal section and in the cross-section of roads or railways and have the same solution described in Sections 3.1 and 3.2.
The proposed methods can also be applied in environment protection, and hydrotechnical engineering [18,19]. The evaluation of the slope is also necessary for river levees. Slope stability analyses are conducted with rising water levels until certain failure is reached. Discharge occurs at a water height of around 7 cm above the crown, which means the slope would actually fail before reaching the hydraulic heads used for slope stability calculation during the overflow [20].

Conclusions
To conclude, the correct estimation of elevation at points has a practical application to calculate the optimum excavation/filling plane when a new house is built. The article presented a method to estimate accurately the elevation of this plane which reduces the costs of embankments.
When designing a new road or railway, designers need to know the ground slope at some points of the future road or rail as exactly as possible, as well as the elevation of each considered point, and, also, the slope of each consecutive two points. Instead of sending workers on the field to perform measurements, these values can be calculated with very good accuracy using the methods described in Sections 2.1, and 2.2. Two algorithms to solve SSTPC2L are presented, a fast and approximate one and an exact but slower one. The error of the first method is low (see Table 1). It has the advantage of being very fast and accurate enough. Thus, if there are many points for which SSTPC2L is applied, the The problem related to determining the slope of natural land at a point on the topographic map is necessary both in the longitudinal section and in the cross-section of roads or railways and have the same solution described in Sections 3.1 and 3.2.
The proposed methods can also be applied in environment protection, and hydrotechnical engineering [18,19]. The evaluation of the slope is also necessary for river levees. Slope stability analyses are conducted with rising water levels until certain failure is reached. Discharge occurs at a water height of around 7 cm above the crown, which means the slope would actually fail before reaching the hydraulic heads used for slope stability calculation during the overflow [20].

Conclusions
To conclude, the correct estimation of elevation at points has a practical application to calculate the optimum excavation/filling plane when a new house is built. The article presented a method to estimate accurately the elevation of this plane which reduces the costs of embankments.
When designing a new road or railway, designers need to know the ground slope at some points of the future road or rail as exactly as possible, as well as the elevation of each considered point, and, also, the slope of each consecutive two points. Instead of sending workers on the field to perform measurements, these values can be calculated with very good accuracy using the methods described in Section 2.1, and Section 2.2. Two algorithms to solve SSTPC2L are presented, a fast and approximate one and an exact but slower one. The error of the first method is low (see Table 1). It has the advantage of being very fast and accurate enough. Thus, if there are many points for which SSTPC2L is applied, the first method can be preferred. However, if we look for a more accurate solution, the second method is more suitable.