Spline quasi-interpolation in the Bernstein basis and its application to digital elevation models

Received: 30 November 2021 Revised: 18 April 2022 Accepted: 13 July 2022 DOI: 10.1002/mma.8602 RESEARCH ARTICLE Spline quasi-interpolation in the Bernstein basis and its application to digital elevation models Francisco J. Ariza-López1Domingo Barrera2Salah Eddargani2 María José Ibáñez2Juan F. Reinoso3 1Department of Cartographic, Geodesic and Photogrammetry Engineering, University of Jaén, Jaén, Spain 2Department of Applied Mathematics, University of Granada, Granada, Spain 3Department of Architectural and Engineering Graphic Expression, University of Granada, Granada, Spain Correspondence Salah Eddargani, Department of Applied Mathematics, University of Granada, Granada, Spain. Email: [email protected] Communicated by: J. Vigo-Aguiar Funding information Spanish State Research Agency, Grant/Award Number: PID2019-106195RB-I00; Universidad de Granada/CBUA A nonstandard low-cost spline approximation method for approximating bivariate functions is constructed. It is applied for Digital Elevation approximation and then its accuracy in the downscaling process is studied. KEYWORDS altimetric error, Bernstein–Bézier coefficients, DEM, upscaling, downscaling, tensor product, quasi-interpolation 1INTRODUCTION Digital Terrain Elevation Models (DETM) are cartographic products with multiple applications in fields such as civil engineering, hydrology, agriculture, environment, and geology.1The quality of the results achieved in each field will largely depend on the data capture method and the transformation of the captured data into a type of DETM. The two most important capture methods currently are photogrammetry and light detection and ranging (LiDAR) for their ability to cover large areas of land that characterizes cartography. To achieve these large territorial areas, data capture is carried out with airborne sensors. The result of the data capture is a point cloud that after processing can be transformed into two types of DETM: An irregular triangle network (TIN) or a regular mesh (DEM). The TIN would be defined by a point number and its corresponding 3D coordinates together with the relationship of adjacent points, which it forms edges of the corresponding triangles with. On the other hand, the DEM would be defined by a matrix in which the value of each cell corresponds to the altitude and the planimetric coordinates are easily calculated knowing the cell size in addition to the coordinates of the first row and column of the matrix. We will focus on DEMs as one of the products delivered by most of the national and regional cartographic agencies of the states. DEMs, although they represent a territorial area, have a discrete character, typical of the matrices through which they are represented. Thus, if it is desired to obtain the slope of the ground, it could be only known in the coordinates associated with the matrix of cell centers. Other variables derived from the DEM could be the orientation, curvature, or flow accumulation, which determines the drainage network, and its computation presents the same problem of spatial uncertainty as the above-mentioned slope. This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made. © 2022 The Authors. Mathematical Methods in the Applied Sciences published by John Wiley & Sons, Ltd. Math Meth Appl Sci. 2023;46:1687–1698. wileyonlinelibrary.com/journal/mma 1687 ARIZA-LÓPEZ ET AL. The discontinuity of DEMs poses a problem when it is necessary to compare variables derived from two DEMs that represent the same territorial space, but each one having a different cell size. Upscaling2and downscaling3–6 are commonly used techniques with regard to this problem, and different approximation methods have been employed for this purpose, such as nearest neighbor, bilinear7or bicubic8interpolation. In this paper, we propose a novel method of function approximation that will allow us to perform both upscaling and downscaling, and we will study its accuracy in the downscaling process. For the function underlying a DEM, we will construct a nonstandard quasi-interpolant that will provide a C2bicubic piecewise surface from which it will be possible to estimate each of the elements required in practice. Quasi-interpolation is a well-known technique for constructing approximants directly from the available information of the function to be approximated, which can be reduced to the values it takes on a set of points. Typically, a quasi-interpolant for a given function will be a linear combination of elements of a set of appropriate non-negative and compactly supported functions, which form a partition of unity. It can be defined to satisfy some required properties.9–13 The nonstandard feature of the quasi-interpolant constructed here comes from the fact that it will be a tensor product quasi-interpolant defined from a univariate non-standard quasi-interpolant, which is constructed on each sub-interval by providing the coefficients in the Bernstein basis. Therefore, on each square of the quadrangulation associated with the DEM the coefficients of the representation of the quasi-interpolant in terms of the corresponding Bernstein polynomials are directly defined from the values at the points in a neighborhood.14–17 In Section 2 a nonstandard univariate quasi-interpolant is constructed for functions defined on the real line, endowed with a uniform partition. Some numerical tests are presented to show the performance of this kind of spline quasi-interpolant. Next, it is appropriately modified to approximate functions defined over an interval. The nonstandard bivariate quasi-interpolant is obtained as the tensor product of the univariate scheme with itself. Numerical tests are presented. In Section 3 the approximation method proposed is applied to a DEM and the altimetric error will be analyzed by comparing altitudes and the planimetric error using the automatic contour line algorithm.18,19 2QUASI-INTERPOLATION IN THE BERNSTEIN BASIS 2.1 The univariate case Suppose that for a real function 𝑓the values 𝑓(vi)and 𝑓(ei),i∈Z, are known, where vi∶= a+ih and ei∶= vi+1 2h,with h>0 the size of the partition Δ∶={vi,i∈Z},anda∈R. We want to construct a cubic spline Q𝑓of 𝑓defined on Δ.Q𝑓reduces on each interval Ii∶= [vi,vi+1],i∈Z,toacubic polynomial, so that it can be expressed in Bernstein's basis relative to Ii,thatis, Q𝑓|Ii(x)= 3 ∑ k=0 bk,iBk(x−vi h),x∈Ii,(1) for some coefficients bk,i∈R,where Bk(t)∶= (3 k)tk(1−t)3−k,t∈[0,1], are the cubic Bernstein polynomials relative to the interval [0,1]. The Bernstein–Bézier (BB) coefficients bk,iof Q𝑓on Ii are linked to the domain points vi+k 3h,k=0,…,3. The union without repetition of all of them gives rise to the set ∶= {a+ih 3,i∈Z}. The BB-coefficients of the quasi-interpolant will be determined to get C2continuity and exactness on the space P2of quadratic polynomials. Since the partition is uniform, one strategy is to start from a partition of by subsets iand appropriately determine the BB-coefficients corresponding to the points in i. Specifically, it is fulfilled that =⋃i∈Zi, with i∶= {ui,vi,wi},whereui∶= vi−h 3and wi∶= vi+h 3. In Figure 1, the above construction is illustrated. The upper subfigure illustrates the structure of the set . In the lower subfigure, the subsets icorresponding to various vertices are shown. According to the notation introduced above, the quasi-interpolant defined in (1) is rewritten as Q𝑓|Ii(x)=c(vi)B0(x−vi h)+c(wi)B1(x−vi h)+c(ui+1)B2(x−vi h)+c(vi+1)B3(x−vi h),x∈Ii,(2) 1688 10991476, 2023, 2, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/mma.8602 by Universidad De Granada, Wiley Online Library on [03/12/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License ARIZA-LÓPEZ ET AL. FIGURE 1 (Top)Subsetofdomain points in for some subintervals induced by Δ. (Bottom) Labeling domain points of various subintervals [Colour figure can be viewed at wileyonlinelibrary.com] where c(p)stands for the BB-coefficient of the domain point p. Figure 1 shows the four domain points whose associated BB-coefficients are involved in the above expression. We propose to define the BB-coefficients of Q𝑓in each interval Iishown in (2) as linear combinations of the values of 𝑓at knots v𝓁and midpoints e𝓁in a small neighborhood of that interval. Problem 1. Find masks 𝛼∶= (𝛼0,𝛼 1,𝛼 2,𝛼 3,𝛼 4),𝛽∶= (𝛽0,𝛽 1,𝛽 2,𝛽 3,𝛽 4)and 𝛾∶= (𝛾0,𝛾 1,𝛾 2,𝛾 3,𝛾 4)such that the quasi-interpolant Q𝑓defined by (2) on each interval I𝓁with BB-coefficients c(u𝓁)∶=𝛼0𝑓(v𝓁−1)+𝛼1𝑓(e𝓁−1)+𝛼2𝑓(v𝓁)+𝛼3𝑓(e𝓁)+𝛼4𝑓(v𝓁+1), c(v𝓁)∶=𝛽0𝑓(v𝓁−1)+𝛽1𝑓(e𝓁−1)+𝛽2𝑓(v𝓁)+𝛽3𝑓(e𝓁)+𝛽4𝑓(v𝓁+1), c(w𝓁)∶=𝛾0𝑓(v𝓁−1)+𝛾1𝑓(e𝓁−1)+𝛾2𝑓(v𝓁)+𝛾3𝑓(e𝓁)+𝛾4𝑓(v𝓁+1) (3) is C2continuous and Q𝑓=𝑓for all 𝑓∈P2. The fact that the partition is uniform leads to the imposition that the masks do not depend on the subinterval in which the quasi-interpolant is calculated. Proposition 2. The quasi-interpolant Q𝑓defined locally by (2) and (3) is C1-continuous if and only if 𝛼k−2𝛽k+𝛾k=0,k∈{0,1,2,3,4}.(4) Proof. Q𝑓is C1-continuous if and only if c(w𝓁)−2c(v𝓁)+c(u𝓁)=0forall𝓁∈Z. Replacing the BB-coefficients c(w𝓁),c(v𝓁),andc(u𝓁)with their expressions given in (3), it holds c(w𝓁)−2c(v𝓁)+c(u𝓁)=E0,1𝑓(v𝓁−1)+E1,1𝑓(e𝓁−1)+E2,1𝑓(v𝓁)+E3,1𝑓(e𝓁)+E4,1𝑓(v𝓁+1), where Ek,1∶= 𝛼k−2𝛽k+𝛾k. This expression will be zero for any function 𝑓if and only if the conditions in the statement are satisfied. Proposition 3. Let us suppose that the quasi-interpolant Q𝑓defined locally by (2) and (3) is C1-regular. Then, it is C2-continuous if and only if 2𝛼0−2𝛾0−𝛾2=0,2𝛼1−2𝛾1−𝛾3=0,𝛼 0+2𝛼2−2𝛾2−𝛾4=0,𝛼 1+2𝛼3−2𝛾3=0,𝛼 2+2𝛼4−2𝛾4=0.(5) Proof. Being Q𝑓C1-continuous, C2-continuity is achieved if and only if c(u𝓁+1)+2c(u𝓁)−2c(w𝓁)−c(w𝓁−1)=0for all 𝓁∈Z. Replacing c(u𝓁+1),c(u𝓁),c(w𝓁),andc(w𝓁−1)with their expressions given in (3), and taking into account that 𝛽kdepends on 𝛼kand 𝛾k,then c(u𝓁+1)+2c(u𝓁)−2c(w𝓁)−c(w𝓁−1)=E0,2𝑓(v𝓁−1)+E1,2𝑓(e𝓁−1)+E2,2𝑓(v𝓁)+E3,2𝑓(e𝓁)+E4,2𝑓(v𝓁+1), where E0,2∶= 2𝛼0−2𝛾0−𝛾2,E1,2∶= 2𝛼1−2𝛾1−𝛾3,E2,2∶= 𝛼0+2𝛼2−2𝛾2−𝛾4,E3,2∶= 𝛼1+2𝛼3−2𝛾3 and E4,2∶= 𝛼2+2𝛼4−2𝛾4. Therefore, Q𝑓is C2-continuous if and only if Ek,2=0, k=0,…,4, and the proof is complete. The following result is easily obtained. Lemma 4. The BB-coefficients in the interval I𝓁of monomials mk(x)∶= (x−vi h)k,k =0,1,2,are(1,1,1,1), (0,1∕3,2∕3,1)and (0,0,1∕3,1), respectively. 1689 10991476, 2023, 2, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/mma.8602 by Universidad De Granada, Wiley Online Library on [03/12/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License ARIZA-LÓPEZ ET AL. We are now in a position to state the main result of this sub-section. Proposition 5. Problem 1 has a unique solution, provided by the masks 𝛼=(0,8 15 ,2 5,4 15 ,−1 5),𝛽=(−1 10 ,2 5,2 5,2 5,−1 10 ) and 𝛾=(−1 5,4 15 ,2 5,8 15 ,0). Proof. By (2), the BB-coefficients relative to interval I𝓁of the quasi-interpolant Q𝑓are (c(v𝓁),c(w𝓁),c(u𝓁+1),c(v𝓁+1)) and are given by (3). The BB-coefficients satisfy conditions (4) and (5) to ensure C2-continuity of Q𝑓. Moreover, exactness on P2is required, so that BB-coefficients of Qmk on I𝓁must be equal to those of mk. They are given in Lemma 4 and those of Qmkare (1,1,1,1), (1 2(−2𝛾0−𝛾1+𝛾3+2𝛾4),1 2(−2𝛽0−𝛽1+𝛽3+2𝛽4),1 2(𝛼1+2𝛼2+3𝛼3+4𝛼4),1 2(𝛾1+2𝛾1+3𝛾3+4𝛾4))and (1 4(4𝛾0+𝛾1+𝛾3+4𝛾4),1 4(4𝛽0+𝛽1+𝛽3+4𝛽4),1 4(𝛼1+4𝛼2+9𝛼3+16𝛼4),1 4(𝛾0+4𝛾1+9𝛾3+16𝛾4)). Therefore, there are 10 linear equations that guarantee C2-continuity and 12 other that produce the required exactness. In total, there are 22 equations in 15 unknowns. A symbolic calculus system makes it possible to prove that such a system of equations has a unique solution, which gives rise to the masks indicated in the statement. For 𝑓∈C3(I𝓁)and for the quasi-interpolation operator Q, given by the unique solution to Problem 1, there exists a constant C𝓁independent of hsuch that ‖𝑓−Q𝑓‖∞,I𝓁≤C𝓁h3‖ ‖ ‖𝑓(3)‖ ‖ ‖∞,I𝓁 . 2.2 Numerical tests in the univariate case To test the performance of the quasi-interpolation scheme defined above, we consider the following three test functions, defined on [0,1], although they can be extended to a larger interval: 1. 𝑓1(x)∶=3 4e−2(9x−2)2−1 5e−(9x−7)2−(9x−4)2+1 2e−(9x−7)2−1 4(9x−3)2+3 4e1 10 (−9x−1)−1 49 (9x+1)2. 2. 𝑓2(x)∶=1 2xcos4(4(x2+x−1)). 3. 𝑓3(x)∶=(3+cos (2𝜋x))−1log (1 x3+1)sin (7𝜋x). Their plots appear in Figure 2. The tests are carried out for a sequence of uniform partition Δnassociated with the vertices vi=ih,i=0,…,n, where h∶= 1 n. FIGURE 2 Plots of univariate functions 𝑓1,𝑓2,and𝑓3[Colour figure can be viewed at wileyonlinelibrary.com] 1690 10991476, 2023, 2, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/mma.8602 by Universidad De Granada, Wiley Online Library on [03/12/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License ARIZA-LÓPEZ ET AL. TABLE 1 Error estimates and NCOs for test functions 𝑓1,𝑓2,and𝑓3 𝒇1𝒇2𝒇3 nError NCO Error NCO Error NCO 16 0.04662822716 −0.03579025054 −0.008119172533 − 32 0.004072919573 3.517 0.005506927133 2.700 0.0007660623721 3.406 64 0.000327887317 3.635 0.0004341112452 3.665 0.00007860374495 3.285 128 0.000033860570 3.276 0.00004277811450 3.343 8.9665358890173 ×10−63.132 256 4.072326029745 ×10−63.056 5.5514176444525 ×10−62.946 1.0827003465402 ×10−63.050 The quasi-interpolation error is estimated as max 𝓁=1,…,200 |𝑓(x𝓁)−𝑓(x𝓁)|, where x𝓁,𝓁=1,…,200, are equally spaced points in [0,1]. The numerical convergence order (NCO) is given by the rate NCO ∶= log2(E(2n) E(n)), where E(m)stands for the estimated error associated with Δm. Table 1 shows the numerical errors and the numerical convergence orders (NCOs) for 𝑓1,𝑓2,and𝑓3. 2.3 Quasi-interpolation of functions defined on an interval When the function 𝑓is defined on an interval [a,b]endowed with a uniform partition with knots xi∶= a+ih,i= 0,…,n,h∶= b−a n, then the above numerical scheme is not applicable since the values 𝑓(a−h),𝑓(a−h 2),𝑓(b+h 2) and 𝑓(b+h), which should appear in (3), are not defined. Next, it is shown how to define those values. Lemma 6. Let us suppose that 𝑓(a−h)and 𝑓(a−h 2)are approximated by 𝜙0𝑓(a)+𝜙1𝑓(a+h 2)+𝜙2𝑓(a+h)and 𝜑0𝑓(a)+𝜑1𝑓(a+h 2)+𝜑2𝑓(a+h). Then, the approximation errors are null when the data come from a quadratic polynomial if and only if 𝜙0=6,𝜙1=−8,𝜙2=3,𝜑0=3,𝜑1=−3,and𝜑2=1. Similarly, if 𝑓(b+h 2)and 𝑓(b+h)are approximated by 𝜑0𝑓(b)+𝜑1𝑓(b−h 2)+𝜑2𝑓(b−h)and 𝜙0𝑓(b)+𝜙1𝑓(b−h 2)+𝜙2𝑓(b−h), then the approximation errors satisfy the same condition regarding exactness if and only if 𝜑0=1,𝜑1=−3,𝜑2=3,𝜙0=3,𝜙1=−8,and𝜙2=6. Proof. The results are obtained by imposing that the errors are zero when 𝑓=1,x,x2, and solving the resulting systems of equations. Substituting in the BB-coefficients corresponding to 𝓁=0and𝓁=n−1 given by (3) and in Proposition 5 the approximations obtained to 𝑓(a−h),𝑓(a−h 2),𝑓(b+h 2)and 𝑓(b+h), we obtain the BB-coefficients in the subintervals [a,a+h]and [b−h,b]. They are given next. Proposition 6. The BB-coefficients of the quasi-interpolant Q𝑓to 𝑓∈C([a,b])in [a,a+h]are 𝑓(a), 4 3𝑓(a+h 2)−1 3𝑓(a+h), 8 15𝑓(a+h 2)+2 5𝑓(a+h)+4 15𝑓(a+3 2h)−1 5𝑓(a+2h), −1 10𝑓(a)+2 5𝑓(a+h 2)+2 5𝑓(a+h)+2 5𝑓(a+3 2h)−1 10𝑓(a+2h). 1691 10991476, 2023, 2, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/mma.8602 by Universidad De Granada, Wiley Online Library on [03/12/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License ARIZA-LÓPEZ ET AL. Those in [b−h,b]are −1 10𝑓(b−2h)+2 5𝑓(b−3 2h)+2 5𝑓(b−h)+2 5𝑓(b−h 2)−1 10𝑓(b), −1 5𝑓(b−2h)+4 15𝑓(b−3 2h)+2 5𝑓(b−h)+8 15𝑓(b−h 2), −1 3𝑓(b−h)+4 3𝑓(b−h 2), 𝑓(b). 2.4 Nonstandard quasi-interpolation of bivariate functions The aim of this subsection is to define a quasi-interpolant for a bivariate function 𝑓(x,𝑦)defined on a square [a,b]× [a,b]from the nonstandard quasi-interpolation operator Qconstructed from the general masks in Proposition 5 and the boundary masks in Proposition 6. Low computational cost is an essential feature of the bivariate numerical approximation scheme to be used in modeling DEMs which, in general, will correspond to large-area terrains, so that Qwill perform well a priori. The 2D quasi-interpolant is defined as the tensor product of Qwith itself. For 𝑦∈[a,b]and x∈Ii,i=1,…,n−2, the quasi-interpolant Q𝑓(·,𝑦 )of 𝑓as a function of the first variable is written as g(𝑦)∶= Bi,0(x)( 𝛽0𝑓(vi−1,𝑦 )+𝛽1𝑓(ei−1,𝑦 )+𝛽2𝑓(vi,𝑦 )+𝛽3𝑓(ei,𝑦 )+𝛽4𝑓(vi+1,𝑦 )) +Bi,1(x)( 𝛾0𝑓(vi−1,𝑦 )+𝛾1𝑓(ei−1,𝑦 )+𝛾2𝑓(vi,𝑦 )+𝛾3𝑓(ei,𝑦 )+𝛾4𝑓(vi+1,𝑦 )) +Bi,2(x)( 𝛼0𝑓(vi,𝑦 )+𝛼1𝑓(ei,𝑦 )+𝛼2𝑓(vi+1,𝑦 )+𝛼3𝑓(ei+1,𝑦 )+𝛼4𝑓(vi+2,𝑦 )) +Bi,3(x)( 𝛽0𝑓(vi,𝑦 )+𝛽1𝑓(ei,𝑦 )+𝛽2𝑓(vi+1,𝑦 )+𝛽3𝑓(ei+1,𝑦 )+𝛽4𝑓(vi+2,𝑦 )) , where Bi,k(x)∶= Bk(x−vi h),k=0,1,2,3. Therefore, for 𝑗=1,…,n−2, the tensor product quasi-interpolant Ti,𝑗 𝑓of 𝑓(x,𝑦 )on Ii×I𝑗is given by Ti,𝑗 𝑓(x,𝑦 )∶= B𝑗,0(𝑦)(𝛽0g(v𝑗−1)+𝛽1g(e𝑗−1)+𝛽2g(v𝑗)+𝛽3g(e𝑗)+𝛽4g(v𝑗+1)) +B𝑗,1(𝑦)(𝛾0g(v𝑗−1)+𝛾1g(e𝑗−1)+𝛾2g(v𝑗)+𝛾3g(e𝑗)+𝛾4g(v𝑗+1)) +B𝑗,2(𝑦)(𝛼0g(v𝑗)+𝛼1g(e𝑗)+𝛼2g(v𝑗+1)+𝛼3g(e𝑗+1)+𝛼4g(v𝑗+2)) +B𝑗,3(𝑦)(𝛽0g(v𝑗)+𝛽1g(e𝑗)+𝛽2g(v𝑗+1)+𝛽3g(e𝑗+1)+𝛽4g(v𝑗+2)), with B𝑗,k(𝑦)∶= Bk(𝑦−v𝑗 h),k=0,1,2,3. After some calculations, Ti,𝑗 is found to be a linear combination of the Bernstein polynomials associated with the square Ii×I𝑗, whose coefficients are determined from the values of f at vertices and midpoints of a neighborhood of that square. Figure 3 shows how the Bernstein polynomials and the corresponding coefficients are arranged in a rectangular structure. FIGURE 3 Arrangement of Bernstein polynomials relative to square Ii×I𝑗and its associated BB-coefficients 1692 10991476, 2023, 2, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/mma.8602 by Universidad De Granada, Wiley Online Library on [03/12/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License ARIZA-LÓPEZ ET AL. More precisely, Ti,𝑗 𝑓(x,𝑦 )= 3 ∑ p=0 3 ∑ q=0 𝜙p,qBp(x−vi h)Bq(𝑦−v𝑗 h),(x,𝑦 )∈Ii×I𝑗,i,𝑗 ∈{1,…,n−2}, with 𝜙0,0=𝛽2 0𝑓(vi−1,v𝑗−1)+𝛽0𝛽1(𝑓(ei−1,v𝑗−1)+𝑓(vi−1,e𝑗−1))+𝛽2 1𝑓(ei−1,e𝑗−1)+𝛽0𝛽2(𝑓(vi−1,v𝑗)+𝑓(vi,v𝑗−1)) +𝛽1𝛽2(𝑓(ei−1,v𝑗)+𝑓(vi,e𝑗−1))+𝛽2 2𝑓(vi,v𝑗)+𝛽0𝛽3(𝑓(ei,v𝑗−1)+𝑓(vi−1,e𝑗)) +𝛽1𝛽3(𝑓(ei−1,e𝑗)+𝑓(ei,e𝑗−1))+𝛽2𝛽3(𝑓(ei,v𝑗)+𝑓(vi,e𝑗))+𝛽2 3𝑓(ei,e𝑗)+𝛽0𝛽4(𝑓(vi−1,v𝑗+1)+𝑓(vi+1,v𝑗−1)) +𝛽1𝛽4(𝑓(ei−1,v𝑗+1)+𝑓(vi+1,e𝑗−1))+𝛽2𝛽4(𝑓(vi,v𝑗+1)+𝑓(vi+1,v𝑗))+𝛽3𝛽4(𝑓(ei,v𝑗+1)+𝑓(vi+1,e𝑗)) +𝛽2 4𝑓(vi+1,v𝑗+1), 𝜙0,1=𝛽0𝛾0𝑓(vi−1,v𝑗−1)+𝛽1𝛾0𝑓(ei−1,v𝑗−1)+𝛽2𝛾0𝑓(vi,v𝑗−1)+𝛽3𝛾0𝑓(ei,v𝑗−1)+𝛽4𝛾0𝑓(vi+1,v𝑗−1)+𝛽0𝛾1𝑓(vi−1,e𝑗−1) +𝛽1𝛾1𝑓(ei−1,e𝑗−1)+𝛽2𝛾1𝑓(vi,e𝑗−1)+𝛽3𝛾1𝑓(ei,e𝑗−1)+𝛽4𝛾1𝑓(vi+1,e𝑗−1)+𝛽0𝛾2𝑓(vi−1,v𝑗)+𝛽1𝛾2𝑓(ei−1,v𝑗) +𝛽2𝛾2𝑓(vi,v𝑗)+𝛽3𝛾2𝑓(ei,v𝑗)+𝛽4𝛾2𝑓(vi+1,v𝑗)+𝛽0𝛾3𝑓(vi−1,e𝑗)+𝛽1𝛾3𝑓(ei−1,e𝑗)+𝛽2𝛾3𝑓(vi,e𝑗) +𝛽3𝛾3𝑓(ei,e𝑗)+𝛽4𝛾3𝑓(vi+1,e𝑗)+𝛽0𝛾4𝑓(vi−1,v𝑗+1)+𝛽1𝛾4𝑓(ei−1,v𝑗+1)+𝛽2𝛾4𝑓(vi,v𝑗+1) +𝛽3𝛾4𝑓(ei,v𝑗+1)+𝛽4𝛾4𝑓(vi+1,v𝑗+1), 𝜙0,2=𝛼0𝛽0𝑓(vi−1,v𝑗)+𝛼1𝛽0𝑓(vi−1,e𝑗)+𝛼2𝛽0𝑓(vi−1,v𝑗+1)+𝛼3𝛽0𝑓(vi−1,e𝑗+1)+𝛼4𝛽0𝑓(vi−1,v𝑗+2) +𝛼0𝛽1𝑓(vi−1,v𝑗+2)+𝛼1𝛽1𝑓(ei−1,e𝑗) +𝛼2𝛽1𝑓(ei−1,v𝑗+1)+𝛼3𝛽1𝑓(ei−1,e𝑗+1)+𝛼4𝛽1𝑓(ei−1,v𝑗+2)+𝛼0𝛽2𝑓(vi,v𝑗)+𝛼1𝛽2𝑓(vi,e𝑗)+𝛼2𝛽2𝑓(vi,v𝑗+1) +𝛼3𝛽2𝑓(vi,e𝑗+1)+𝛼4𝛽2𝑓(vi,v𝑗+2)+𝛼0𝛽3𝑓(ei,v𝑗)+𝛼1𝛽3𝑓(ei,e𝑗)+𝛼2𝛽3𝑓(ei,v𝑗+1)+𝛼3𝛽3𝑓(ei,e𝑗+1) +𝛼4𝛽3𝑓(ei,v𝑗+2)+𝛼0𝛽4𝑓(vi+1,v𝑗)+𝛼1𝛽4𝑓(vi+1,e𝑗)+𝛼2𝛽4𝑓(vi+1,v𝑗+1)+𝛼3𝛽4𝑓(vi+1,e𝑗+1)+𝛼4𝛽4𝑓(vi+1,v𝑗+2), 𝜙0,3=𝛽2 0𝑓(vi−1,v𝑗)+𝛽0𝛽1(𝑓(ei−1,v𝑗)+𝑓(vi−1,e𝑗))+𝛽2 1𝑓(ei−1,e𝑗)+𝛽0𝛽2(𝑓(vi−1,v𝑗+1)+𝑓(vi,v𝑗)) +𝛽1𝛽2(𝑓(ei−1,v𝑗+1)+𝑓(vi,e𝑗)) +𝛽2 2𝑓(vi,v𝑗+1)+𝛽0𝛽3(𝑓(ei,v𝑗)+𝑓(vi−1,e𝑗+1))+𝛽1𝛽3(𝑓(ei−1,e𝑗+1)+𝑓(ei,e𝑗)) +𝛽2𝛽3(𝑓(ei,v𝑗+1)+𝑓(vi,e𝑗+1))+𝛽2 3𝑓(ei,e𝑗+1) +𝛽0𝛽4(𝑓(vi−1,v𝑗+2)+𝑓(vi+1,v𝑗))+𝛽1𝛽4(𝑓(ei−1,v𝑗+2)+𝑓(vi+1,e𝑗))+𝛽2𝛽4(𝑓(vi,v𝑗+2)+𝑓(vi+1,v𝑗+1)) +𝛽3𝛽4(𝑓(ei,v𝑗+2)+𝑓(vi+1,e𝑗+1))+𝛽2 4𝑓(vi+1,v𝑗+2). Coefficients 𝜙i,𝑗 ,i=1,2,3and𝑗=0,1,2,3, have similar expressions. These values of the coefficients correspond to a square having no side on the boundary. For the remainder, a similar procedure using the masks of the univariate quasi-interpolant attached to the boundary leads to similar expressions. Next, the performance of this nonstandard quasi-interpolation operator is illustrated by considering the Franke and Nielson test functions: g1(x,𝑦)∶=1 2exp (−((9x−7)2+1 4(9𝑦−3)2))+3 4exp (−1 49(9x+1)2−1 10(9𝑦+1)) −1 5exp (−(9x−4)2−(9𝑦−7)2)+3 4exp (−((9x−2)2+(9𝑦−2)2)), g2(x,𝑦)∶=1 2𝑦cos4(4(x2+𝑦−1)). Figure 4 shows the plots of Franke and Nielson functions. Table 2 shows the errors and NCOs for h=2−r,r=4,…,8. They are in good agreement with the theoretical results. 1693 10991476, 2023, 2, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/mma.8602 by Universidad De Granada, Wiley Online Library on [03/12/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License ARIZA-LÓPEZ ET AL. FIGURE 4 From left to right, plots of functions g1and g2[Colour figure can be viewed at wileyonlinelibrary.com] TABLE 2 Error estimates and NCOs for test functions g1and g2g1g2 n error NCO error NCO 16 0.3625040 −0.257841 − 32 0.0700742 2.37104 0.0489511 2.39707 64 0.0103237 2.76293 0.00712815 2.77974 128 0.00118445 3.12366 0.000912965 2.96490 256 0.000134193 3.14184 0.000112483 3.02085 FIGURE 5 Georeferenced area of our DEM (gray area represents the hillshade map covered by the DEM) [Colour figure can be viewed at wileyonlinelibrary.com] 3MODELING OF A DEM FROM TENSOR PRODUCT QUASI-INTERPOLATION IN THE BERNSTEIN BASIS 3.1 Material and methodology Our DEM reference (DEMref2×2)willbetheDigital Terrain Model MDT02 produced by the Instituto Geográfico Nacional of Spain with a cell size equal to 2m ×2m. The coordinate system uses the geodetic reference system ETRS89 and UTM projection UTM zone 30N (see Figure 5). The DEMref2×2size is 4096m ×4096m centered at the UTM coordinates (583023m,4718590m). The goodness of the proposed approximation algorithm based on quasi-interpolation in the Bernstein basis is assessed producing a DEM (DEMder) of larger cell sizes than the DEMref2×2. Then a surface is adjusted and it is generated a resampled DEM (DEMres) with the same cell size than DEMref2×2cell. The planimetric and altimetric discrepancies between DEMder and DEMres are calculated separately. The process is divided into the following phases: 1694 10991476, 2023, 2, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/mma.8602 by Universidad De Granada, Wiley Online Library on [03/12/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License ARIZA-LÓPEZ ET AL. 1. The DEMref2×2is resampled by the closest neighbor method at 3 different resolution levels (4m ×4m, 8m ×8m and 16m ×16m), which produces three different DEMresp which we call, respectively, DEMresp4×4,DEM resp8×8,and DEMresp16×16. They can be represented generically as DEMrespX×X. 2. The quasi-interpolation operator fits a surface on each of the DEMrespX×X, denoting such a surface as SrespX×X. 3. Since the definition domain SrespX×Xis the same as DEMref2×2,S respX×Xcan be evaluated at the same points that define the DEMref2×2matrix, with which we obtain a homologous DEMder2×2_from_X×XofthesamesizeasDEM ref2×2. 4. The planimetric quality is assessed using the method based on contour lines introduced in Reinoso19:theresults given by quasi-interpolation are compared with those produced by the bicubic interpolation.8 5. The altimetric quality is assessed by calculating the mean discrepancy in absolute value between the altitudes of DEMref2×2and DEMder2×2_from_X×X. 3.2 Planimetric error assessment The planimetric error assessment is based on the area enclosed between homologous contours computed on both DEMs, i.e. DEMref2×2and DEMder2×2_from_X×X: 1. The contours of both DEMs are computed (Figure 6A,B). Homologous contours are matched automatically (curves C4aand C4bin Figure 6A,B). Homologous contours are those with the same height and placed in approximately the same place (see Reinoso19 to understand how the algorithm automatically eliminate the ambiguity). 2. Overlapping homologous contours (Figure 6C) lead to areas enclosed (gray area in Figure 6D). The planimetric error (Pei)referredtotheith pair of homologous contours (Cia,Cib) is computed as the surface enclosed between both curves (Si) divided by the mean length of those contours (Lmi=1 2(Lia+Lib).i.e.Pei=Si∕Lmi. 3. The mean planimetric error linked to the DEMder2×2_from_X×Xis denoted as PeDEMder2x2_from_X×Xand its value is the weighted average of the Peiof all the homologous contours. The weighting factor is the average length of homologous contours, divided by the total length of the average contours being the total length LTot =∑n i=1Lmi,thatis, PeDEMder2x2_from_X×X=1 LTot n ∑ i=1 Pei·Lmi. FIGURE 6 Homologous contours and area between them 1695 10991476, 2023, 2, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/mma.8602 by Universidad De Granada, Wiley Online Library on [03/12/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License