# Piecewise Bivariate Hermite Interpolations for Large Sets of Scattered Data

Authors: Renzhong Feng; Yanan Zhang
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

The requirements for interpolation of scattered data are high accuracy and high efficiency. In this paper, a piecewise bivariate Hermite interpolant satisfying these requirements
is proposed. We firstly construct a triangulation mesh using the given scattered point set. Based on this mesh, the computational point (x,y) is divided into two types: interior point and exterior point. The value of Hermite interpolation polynomial on a triangle will be used as the approximate value if point (x,y) is an interior point, while the value of a Hermite interpolation function with the form of weighted combination will be used if it is an exterior point. Hermite interpolation needs the first-order derivatives of the interpolated function which is not directly given in scatted data, so this paper also gives the approximate derivative at every scatted point using local radial basis function interpolation. And numerical tests indicate that the proposed piecewise bivariate Hermite interpolations are economic and have good approximation capacity.


---

## Body

## 1. Introduction

The approximation to higher dimensional scattered data is one of the hot and difficult problems in approximation theory field. Due to the characteristics of scattered data, such as the large amount, irregularity, and high dimensionality, it is difficult to construct the approximation methods for them. However, the approximation method to scattered data has been widely applied in many fields, that is, estimating the parameter of input nonlinear system [1–3], solving partial differential equations, surface reconstruction in reverse engineering, data visualization, and so on. The most frequently-used approximation methods include interpolation by spline, interpolation by radial basis function, and the least square approximation. The interpolations by spline have two types, which are global interpolation and local interpolation, respectively. The global spline interpolation [4] is not able to deal with large scale of scattered data, while the local spline interpolation [5–7] needs smooth joining between piece and piece. The radial basis function interpolation [8, 9] requires solving a large scale of linear system, and the least square approximation [10, 11] also requires solving a certain scale of linear system. This is a problem that doesnot easily work out in the computation. It is well known that Hermite interpolation has higher approximation accuracy than interpolation that only interpolates function values. This is because it interpolates not only the function value, but also the derivative value. However, how to construct multivariate Hermite interpolation for large sets of scattered data is an important research consideration which is also a difficult task. At present, some results have been published [12–14]; however, these methods all require solving a certain scale of linear system. Therefore, they are not suitable for approximating to large scale of scattered data.This paper proposes two piecewise bivariate Hermite interpolation methods for large sets of scattered data. One of them uses exact derivative in Hermite interpolation, and the other one does the approximate derivative. Both methods are economic, free us from solving any linear systems, and have better approximation capacity. Therefore, they are especially suitable for the approximation to large sets of scattered data. These interpolation methods firstly use Delaunay triangulation method [15] to make a triangulation mesh based on the given scattered point set and then divide the computational point (x,y) into two types of interior point and exterior outer point. If point (x,y) is an interior point, then the value of the Hermite interpolation function on the triangle which point (x,y) lies on is regarded as the approximate value of the approximated function, while if point (x,y) is an exterior point, then the value of a Hermite interpolation with a form of weighted combination is regarded as the approximate value. Hermite interpolation needs the first-order derivatives of the interpolated function; however, the derivative information is not directly given at scatted data, so this paper gives the approximate derivative at every scatted point by using local radial basis function interpolation.The structure of this paper is organised as follows. The Hermite interpolation on triangle is presented in Section2. The construction of piecewise bivariate Hermite interpolation is presented in Section 3. The estimation of partial derivative at every scattered point by local radial basis function interpolation is presented in Section 4. Some numerical tests and comparisons between constructed schemes and other interpolation schemes in accuracy and efficiency are carried out in Section 5. This paper ends with a section of brief conclusion.


## 2. Hermite Interpolation Based on Triangle

Given the values and partial derivatives of a bivariate functionf(x,y) at three vertices (xi,yi),i=1,2,3 of a triangle T123 (their vertices 1, 2, 3 are arranged counter to clockwise), that is, fi,(fx)i,(fy)i,  i=1,2,3 (see Figure 1). A bivariate Hermite interpolation polynomial H123(x,y) is to be constructed satisfying the following interpolation conditions:
(1)H123(xi,yi)=fi,∂H123∂x(xi,yi)=(fx)i,∂H123∂y(xi,yi)=(fy)i,i=1,2,3.Figure 1
Interpolation triangleT123.In order to infer the expression of Hermite interpolation from interpolation conditions (1), the barycenter coordinate of a triangle needs to be introduced first. It is assumed that P(x,y) is a point on triangle T123 and its barycentric coordinates is denoted as (l1,l2,l3). Then the relation between its cartesian coordinates and barycentric coordinates is
(2)l1=(x2y3-x3y2)+(y2-y3)x+(x3-x2)y2S,l2=(x3y1-x1y3)+(y3-y1)x+(x1-x3)y2S,l3=(x1y2-x2y1)+(y1-y2)x+(x2-x1)y2S,l1+l2+l3=1,l1,l2,l3≥0,x=x1l1+x2l2+x3l3,y=y1l1+y2l2+y3l3,
where S is the area of the triangle T123 and
(3)S=|1x1y11x2y21x3y3|.
Making use of transformation (2), the interpolation conditions (1) can be changed into the following form under barycentric coordinates:
(4)H123(1,0)=f1,H123(0,1)=f2,H123(0,0)=f3,∂H123∂l1(1,0)=(x1-x3)(fx)1+(y1-y3)(fy)1,∂H123∂l1(0,1)=(x1-x3)(fx)2+(y1-y3)(fy)2,∂H123∂l1(0,0)=(x1-x3)(fx)3+(y1-y3)(fy)3,∂H123∂l2(1,0)=(x2-x3)(fx)1+(y2-y3)(fy)1,∂H123∂l2(0,1)=(x2-x3)(fx)2+(y2-y3)(fy)2,∂H123∂l2(0,0)=(x2-x3)(fx)3+(y2-y3)(fy)3.A complete bivariate polynomial of three degreesp3(x,y)=∑i+j=03‍cijxiyj has ten undetermined coefficients. Under barycentric coordinate its expression is
(5)p3(x,y)=c1l13+c2l23+c3l33+c4l12l2+c5l12l3+c6l22l1+c7l22l3+c8l32l1+c9l32l2+c10l1l2l3.Since the Hermite interpolation conditions (1) or (4) only provide 9 conditions, we construct a Hermite interpolation polynomial with nine coefficients
(6)Hf(l1,l2)=c1l13+c2l23+c3l33+c4l12l2+c5l12l3+c6l22l1+c7l22l3+c8l32l1+c9l32l2+0.5(c7+c8+c9)l1l2l3,
which meets interpolation condition (4). Then the values of the coefficients in formula (6) can be obtained:
(7)c1=f1,c2=f2,c3=f3,c4=3f1-∂H123∂l1(1,0)+∂H123∂l2(1,0),c5=3f1-∂H123∂l1(1,0),c6=3f2-∂H123∂l2(0,1)+∂H123∂l1(0,1),c7=3f2-∂H123∂l2(0,1),c8=3f3+∂H123∂l1(0,0),c9=3f3+∂H123∂l2(0,0).Combining (4), (6), and (7), the Hermite interpolation polynomial satisfying (4) can be written as
(8)H123(x,y)=∑i=13‍[αi(x,y)fi+βi(x,y)(fx)i+γi(x,y)(fy)i],
where
(9)α1(x,y)=l13+3l2l12+3l3l12+2l1l2l3,α2(x,y)=l23+3l1l22+3l3l22+2l1l2l3,α3(x,y)=l33+3l12l32+3l2l32+2l1l2l3,β1=(x2-x1)(l12l2+0.5l1l2l3)+(x3-x1)(l12l3+0.5l1l2l3),β2=(x1-x2)(l22l1+0.5l1l2l3)+(x3-x2)(l22l3+0.5l1l2l3),β3=(x1-x3)(l32l1+0.5l1l2l3)+(x2-x3)(l32l2+0.5l1l2l3),γ1=(y2-y1)(l12l2+0.5l1l2l3)+(y3-y1)(l12l3+0.5l1l2l3),γ2=(y1-y2)(l22l1+0.5l1l2l3)+(y3-y2)(l22l3+0.5l1l2l3),γ3=(y1-y3)(l32l1+0.5l1l2l3)+(y2-y3)(l32l2+0.5l1l2l3).


## 3. Bivariate Piecewise Hermite Interpolation for Large Sets of Scattered Data

Given a set of scattered data{(xj,yj,fj)}j=1N⊂R3, which are assumed to be sampled from a functionf(x,y),(x,y)∈Ω. The partial derivatives of f(x,y) at the set of points {(xj,yj)}j=1N⊂Ω⊂R2 are given as (fx)j,(fy)j,j=1,…,N. Taking the N scattered points as nodes, a triangulation mesh, T, is constructed in domain Ω using Delaunay triangulation method. If point (x,y)∈Ω⊂R2 lies on a triangle of T, it is called as an interior point, otherwise, an exterior point. For example, points Q1,Q2 in Figure 2 are interior points, while point Q3 is an exterior point. Next we will give the approximate value Qf(x,y) of f(x,y) according to the category of point (x,y).Figure 2
The exhibition of interior and exterior points.


### 3.1. Interior Point

Suppose that point(x,y) lies on a triangle Tjkl⊂T whose vertices are (xj,yj),(xk,yk),(xl,yl). We now construct a Hermite interpolation polynomial as (8) on triangle Tjkl satisfying the following conditions:
(10)Hjkl(xj,yj)=fj,∂Hjkl∂x(xj,yj)=(fx)j,∂Hjkl∂y(xj,yj)=(fy)j,Hjkl(xk,yk)=fk,∂Hjkl∂x(xk,yk)=(fx)k,∂Hjkl∂y(xk,yk)=(fy)k,Hjkl(xl,yl)=fl,∂Hjkl∂x(xl,yl)=(fx)l,∂Hjkl∂y(xl,yl)=(fy)l.
Then let Qf(x,y)=Hjkl(x,y), so that Qf(x,y)≈f(x,y).


### 3.2. Exterior Point

When point(x,y) doesnot belong to any triangle within T, we construct the following form of Hermite interpolation:
(11)Qf(x,y)=∑j=1N‍Wj(x,y)Hf,j(x,y),
satisfying interpolation conditions
(12)Qf(xj,yj)=fj,∂Qf∂x(xj,yj)=(fx)j,∂Qf∂y(xj,yj)=(fy)j,j=1,2,…,N.In this paper, we takeHf,j(x,y) as a Hermite interpolation polynomial with the same form of (8) on triangle Tjkl with the least area among all triangles with vertex (xj,yj), satisfying interpolation conditions (10). Hf,j(x,y) is called as node basis function and function Wj(x,y) is called as weighted function. In order to satisfy interpolation conditions (12), these weighted functions are required for satisfying conditions
(13)Wj(xi,yi)=δij,i,j=1,…,N,(14)∑j=1N‍Wj(x,y)≡1,Wj(x,y)≥0,j=1,…,N,(15)∂Wj∂x(xi,yi)=∂Wj∂y(xi,yi)=0,i,j=1,…,N.In the paper we take the weighted functions in (11) as those in the literature [16]:
(16)Wj(x,y)=(1/ρj2)∑k=1N(1/ρk2),
where
(17)1ρk=(Rk-dk)+Rkdk,(Rk-dk)+={Rk-dkRk-dk≥0,0Rk-dk<0,dk=dk(x,y)=(x-xk)2+(y-yk)2.The value ofRk in (17) can be selected by user properly. It can be seen that when (x,y)→(xj,yj), (1/ρj2)→+∝, thereby Wj(x,y)→1; when (x,y)→(xk,yk),k≠j and (1/ρk2)→+∝, thereby Wj(x,y)→0. Thus Wj satisfies condition (13). Wj in formula (16) satisfies condition (14) obviously, and the proof of which can be found in literature [16]. If 1/ρj is determined by (17), then when dj(x,y)>Rj,Wj(x,y)=0. So in case of Rj being appropriately selected (Rj is also called as influence radius of node basis function), there are only several nonzero terms in sum ∑k=1N‍Wk(x,y). Then it can be found that interpolation function (11) which takes functions (16) and (17) as its weight functions is a local Hermite interpolation. When applying (11), we can use different influence radius at every node, so we can use the same radius R. The following computational formulation of uniform influence radius R in (18) is from [16]:
(18)R=D2NwN,D=maxi,jdi(xj,yj).NW is the number of points used inside a circle of radius R.


### 3.3. Bivariate Piecewise Hermite Interpolation

Summing up Sections3.1 and 3.2, a bivariate piecewise Hermite interpolation function is constructed on Ω to approximate the function f(x,y)(19)Qf(x,y)={Hjkl(x,y),if(x,y)∈Γ,∑j=1N‍Wj(x,y)Hf,j(x,y),if(x,y)∈Ω-Γ,
where Tjkl is a triangle within T with node  (xj,yj),(xk,yk),(xl,yl) as its vertices and Hjkl is a Hermite interpolation polynomial on triangle Tjkl.We now summarize the description of the constructed interpolation function, written as Algorithm1.Algorithm 1.
Consider the following.
Step  1. Generate a triangle mesh T in domain Ω using Delaunay triangulation method based on given scattered point set {(xj,yj)}j=1N⊂Ω.
Step  2. Judge the category of point (x,y): interior point or exterior point.
Step  3. If point (x,y) is an interior point, then find out the triangle Tjkl which includes point (x,y), compute the value of the interpolation function Hjkl at point (x,y) using formulas (2), (8), and (10), and take the value as the value of Qf(x,y).
Step  4. If point (x,y) is an exterior point, then compute D=maxi,jdi(xj,yj) and select NW to define R=(D/2)NW/N. The default value of NW is set to 9, and it responds to an uniform radius. Find out the points (xjl,yjl),l=1,…,k from scattered point set {(xj,yj)}j=1N, which belongs to the circle with radius R and center (x,y), compute the values Wjl(x,y),l=1,…,k and Hf,jl(x,y),l=1,…,k, give the sum ∑l=1k‍Wjl(x,y)Hf,jl(x,y), and assign it to Qf(x,y).


## 3.1. Interior Point

Suppose that point(x,y) lies on a triangle Tjkl⊂T whose vertices are (xj,yj),(xk,yk),(xl,yl). We now construct a Hermite interpolation polynomial as (8) on triangle Tjkl satisfying the following conditions:
(10)Hjkl(xj,yj)=fj,∂Hjkl∂x(xj,yj)=(fx)j,∂Hjkl∂y(xj,yj)=(fy)j,Hjkl(xk,yk)=fk,∂Hjkl∂x(xk,yk)=(fx)k,∂Hjkl∂y(xk,yk)=(fy)k,Hjkl(xl,yl)=fl,∂Hjkl∂x(xl,yl)=(fx)l,∂Hjkl∂y(xl,yl)=(fy)l.
Then let Qf(x,y)=Hjkl(x,y), so that Qf(x,y)≈f(x,y).


## 3.2. Exterior Point

When point(x,y) doesnot belong to any triangle within T, we construct the following form of Hermite interpolation:
(11)Qf(x,y)=∑j=1N‍Wj(x,y)Hf,j(x,y),
satisfying interpolation conditions
(12)Qf(xj,yj)=fj,∂Qf∂x(xj,yj)=(fx)j,∂Qf∂y(xj,yj)=(fy)j,j=1,2,…,N.In this paper, we takeHf,j(x,y) as a Hermite interpolation polynomial with the same form of (8) on triangle Tjkl with the least area among all triangles with vertex (xj,yj), satisfying interpolation conditions (10). Hf,j(x,y) is called as node basis function and function Wj(x,y) is called as weighted function. In order to satisfy interpolation conditions (12), these weighted functions are required for satisfying conditions
(13)Wj(xi,yi)=δij,i,j=1,…,N,(14)∑j=1N‍Wj(x,y)≡1,Wj(x,y)≥0,j=1,…,N,(15)∂Wj∂x(xi,yi)=∂Wj∂y(xi,yi)=0,i,j=1,…,N.In the paper we take the weighted functions in (11) as those in the literature [16]:
(16)Wj(x,y)=(1/ρj2)∑k=1N(1/ρk2),
where
(17)1ρk=(Rk-dk)+Rkdk,(Rk-dk)+={Rk-dkRk-dk≥0,0Rk-dk<0,dk=dk(x,y)=(x-xk)2+(y-yk)2.The value ofRk in (17) can be selected by user properly. It can be seen that when (x,y)→(xj,yj), (1/ρj2)→+∝, thereby Wj(x,y)→1; when (x,y)→(xk,yk),k≠j and (1/ρk2)→+∝, thereby Wj(x,y)→0. Thus Wj satisfies condition (13). Wj in formula (16) satisfies condition (14) obviously, and the proof of which can be found in literature [16]. If 1/ρj is determined by (17), then when dj(x,y)>Rj,Wj(x,y)=0. So in case of Rj being appropriately selected (Rj is also called as influence radius of node basis function), there are only several nonzero terms in sum ∑k=1N‍Wk(x,y). Then it can be found that interpolation function (11) which takes functions (16) and (17) as its weight functions is a local Hermite interpolation. When applying (11), we can use different influence radius at every node, so we can use the same radius R. The following computational formulation of uniform influence radius R in (18) is from [16]:
(18)R=D2NwN,D=maxi,jdi(xj,yj).NW is the number of points used inside a circle of radius R.


## 3.3. Bivariate Piecewise Hermite Interpolation

Summing up Sections3.1 and 3.2, a bivariate piecewise Hermite interpolation function is constructed on Ω to approximate the function f(x,y)(19)Qf(x,y)={Hjkl(x,y),if(x,y)∈Γ,∑j=1N‍Wj(x,y)Hf,j(x,y),if(x,y)∈Ω-Γ,
where Tjkl is a triangle within T with node  (xj,yj),(xk,yk),(xl,yl) as its vertices and Hjkl is a Hermite interpolation polynomial on triangle Tjkl.We now summarize the description of the constructed interpolation function, written as Algorithm1.Algorithm 1.
Consider the following.
Step  1. Generate a triangle mesh T in domain Ω using Delaunay triangulation method based on given scattered point set {(xj,yj)}j=1N⊂Ω.
Step  2. Judge the category of point (x,y): interior point or exterior point.
Step  3. If point (x,y) is an interior point, then find out the triangle Tjkl which includes point (x,y), compute the value of the interpolation function Hjkl at point (x,y) using formulas (2), (8), and (10), and take the value as the value of Qf(x,y).
Step  4. If point (x,y) is an exterior point, then compute D=maxi,jdi(xj,yj) and select NW to define R=(D/2)NW/N. The default value of NW is set to 9, and it responds to an uniform radius. Find out the points (xjl,yjl),l=1,…,k from scattered point set {(xj,yj)}j=1N, which belongs to the circle with radius R and center (x,y), compute the values Wjl(x,y),l=1,…,k and Hf,jl(x,y),l=1,…,k, give the sum ∑l=1k‍Wjl(x,y)Hf,jl(x,y), and assign it to Qf(x,y).


## 4. The Estimation of Partial Derivative

In the process of using bivariate piecewise Hermite interpolationQf, it demands the first-order derivatives of the interpolated function f(x,y) at every point, but the scattered data {(xj,yj,fj)}j=1N  doesnot provide the derivatives. Thus we use the first-order derivative of a local radial basis interpolation to approximate the derivative of the interpolated function f(x,y) and replace this exact derivative in (19) with the approximate derivative.The Delaunay triangulationT is a convex subdomain of Ω. Every scattered point Xj=(xj,yj) is either a boundary point or an interior node point of T. We put the vertices of all triangles with point (xj,yj) as into a set, written as Aj=(xj,yj)∪{(xjl,yjl)}l=1m, see Figure 3.The local radial basis interpolation point setAj around point Xj, (a) Xj is a boundary point of T and (b) Xj is an interior node point of T.
(a)
(b)We take the point setAj as local radial basis function interpolation point set and multiquadric radial basis functions
(20)ϕj=(x-xj)2+(y-yj)+c2,ϕj1=(x-xj1)2+(y-yj1)2+c2,⋮ϕjm=(x-xjm)2+(y-yjm)2+c2
as local interpolation basis functions.By solving the linear systemGC=F,
(21)G=(ϕj(xj,yj)ϕj1(xj,yj)…ϕjm(xj,yj)ϕj(xj1,yj1)ϕj1(xj1,yj1)…ϕjm(xj1,yj1)…⋯……ϕj(xjm,yjm)ϕj1(xjm,yjm)…ϕjm(xjm,yjm)),CT=(c0c1c2…cm),FT=(fjfj1fj2…fjm),
we obtained a local radial basis interpolation function pj(x,y)=c0ϕj(x,y)+c1ϕj1(x,y)+c2ϕj2(x,y)+⋯+cmϕjm(x,y), and then compute the first order derivatives of pj(x,y) at Xj=(xj,yj), that is,
(22)∂pj∂x(xj,yj)=c1(xj-xj1)ϕj1(xj,yj)+c2(xj-xj2)ϕj2(xj,yj)+⋯+cm(xj-xjm)ϕjm(xj,yj),∂pj∂y(xj,yj)=c1(yj-yj1)ϕj1(xj,yj)+c2(yj-yj2)ϕj2(xj,yj)+⋯+cm(yj-yjm)ϕjm(xj,yj).The bivariate piecewise interpolation functionQf using the previous approximate derivatives is written as
(23)Q~f(x,y)={H~jkl(x,y),if(x,y)∈Γ,∑j=1N‍Wj(x,y)H~f,j(x,y),if(x,y)∈Ω-Γ.
Algorithm 2 summarizes the computational process of interpolation Q~f(x,y).Algorithm 2.
Consider the following.
Step  1. Generate a triangle mesh T in Ω using Delaunay triangulation method based on the given scattered point set {(xj,yj)}j=1N⊂Ω.
Step  2. Find out the local radial basis interpolation point set Aj at every point (xj,yj), and compute the approximate derivatives (∂pj/∂x)(xj,yj),(∂pj/∂y)(xj,yj) at every point (xj,yj) using interpolation formulation (21) and (22).
Step  3. Judge the category of point (x,y): interior point or exterior point.
Step  4. If point (x,y) is an interior point, then find out the triangle Tjkl which point (x,y) falls in, compute the value of the interpolation function Hjkl at point (x,y) using formulas (2), (8), and (10), and assign the value to the function Q~f(x,y).
Step  5. If point (x,y) is an exterior point, then compute D=maxi,jdi(xj,yj) and select NW to define R=(D/2)NW/N. The default value of NW is set to 9, it responds to an uniform radius. Find out the points (xjl,yjl),l=1,…,k from scattered point set {(xj,yj)}j=1N, which belongs to the circle with radius R and center (x,y), compute the values Wjl(x,y),l=1,…,k and H~f,jl(x,y),l=1,…,k, give the sum ∑l=1k‍Wjl(x,y)H~f,jl(x,y), and assign it to Q~f(x,y).


## 5. Numerical Test

In this section,Qf(x,y) and Q~f(x,y) are used to approximate Franke function
(24)f(x,y)=0.75exp(-(9x-2)2+(9y-2)24)×0.75exp(-(9x+1)210-(9y+1)10)-0.2exp(-(9x-4)2-(9y-7)2)+0.5exp(-(9x-7)2+(9y-3)24),0≤x,y≤1.
Firstly, a scattered data set {(xj,yj,fj)}j=1N is sampled from f(x,y); then, the interpolation functions Qf(x,y) and Q~f(x,y) are generated using the sampled data. Aiming at different sampled data sets, mean square error (MSE) and maximum error (MME) of the two approximate functions are computed and their approximate capacities are also compared. Meanwhile, the approximate accuracy and computational efficiency of Qf(x,y) and Q~f(x,y) are compared with Shepard interpolation [17] and Weighted Local RBF interpolant (abbr. WLRBF) [9]. Finally, the CPU time of these four methods is compared which are devised by MATLAB (7.12.0 R2011a) installed in a computer with Processor: Intel(R) Q8400 2.66 GHz, RAM: 4 GB. The computational results show that the methods introduced in this paper need less CPU time and have higher accuracy.Figure4 shows the 300 scattered points which are randomly selected from [0,1]×[0,1], and Figure 5 is Delaunay triangulation based on Figure 4. Figure 6 presents the graph of Franke function. Figures 7, 8, 9, and 10 describe the interpolation function graphs for the four methods of Qf(x,y), Q~f(x,y), WLRBF, and Shepard which are also based on the 300 scattered points shown in Figure 4. It can be seen that Qf(x,y), Q~f(x,y), and WLRBF can accurately approach Franke function and are better than Shepard method. Figures 11, 12, 13, and 14 display the error graphs of the four methods based on 50*50 test points. Comparing Figure 11 with Figure 13, we can see that Qf(x,y) is much better than WLRBF method. Comparing Figures 12 and 13 for the interior points, Q~f(x,y) is obviously better than WLRBF method. Besides, with the number of scattered data increasing, it is found that approximation capability of Q~f(x,y) is closer to that of Qf(x,y). In case of large set of scattered points, Q~f(x,y) can be hence absolutely superior to WLRBF.Figure 4
300 scattered points randomly selected in[0,1]×[0,1].Figure 5
Delaunay triangulation based on Figure4.Figure 6
The graph of Franke function.Figure 7
The graph ofQf(x,y) based on Figure 4.Figure 8
The graph ofQ~f(x,y) based on Figure 4.Figure 9
The graph of WLRBF based on Figure4.Figure 10
The graph of Shepard interpolation based on Figure4.Figure 11
The error graph ofQf(x,y).Figure 12
The error graph ofQ~f(x,y).Figure 13
The error graph of WLRBF.Figure 14
The error graph of Shepard interpolation.It can be seen from Table1 that the approximation accuracy of Qf(x,y) is higher than that of Q~f(x,y) at the same scattered points. Meanwhile, MSE and MME of the two methods decrease with the increase of the point number. Besides, the approximating accuracies to f(x,y) of Qf(x,y) and Q~f(x,y) get closer at the same time. These results indicate that the local radial function interpolation method can be used to estimate the partial derivative when the number of data points is large enough. Therefore, Q~f(x,y) can be used to approximate f(x,y) when the partial derivative information of the scattered data is unknown.Table 1
The comparison of MSE and MME ofQf(x, y) and Q~f(x, y) for the same scattered points.
Number of scatter data
N
W
Influence radiusR
Q
f
(
x
,
y
)
-
f
(
x
,
y
)
Q
~
f
(
x
,
y
)
-
f
(
x
,
y
)
MSE
MME
MSE
MME
100
9
0.1818
6.3064
e
-
005
0.0578
1.9664
e
-
004
0.0951
300
9
0.1179
1.2890
e
-
006
0.0110
3.8058
e
-
006
0.0162
500
9
0.0918
1.0176
e
-
007
0.0030
3.0822
e
-
007
0.0043
800
9
0.0705
2.0574
e
-
008
0.0012
1.0705
e
-
007
0.0023
1000
9
0.0654
1.2458
e
-
008
0.0011
5.2834
e
-
008
0.0022It can be seen from Table2 that the errors of approximation functions Qf(x,y) and Q~f(x,y) increase with the radius of influence decreasing.Table 2
The influence of influence radius onQf(x,y) and Q~f(x,y).
Number of scatter data
N
W
Influence radiusR
Q
f
(
x
,
y
)
-
f
(
x
,
y
)
Q
~
f
(
x
,
y
)
-
f
(
x
,
y
)
MME
MME
300
4
0.0786
0.0489
0.0906
300
9
0.1179
0.0352
0.0747
300
12
0.1361
0.0325
0.0504
300
16
0.1572
0.0307
0.0481
300
25
0.1965
0.0293
0.0447Tables1 and 3 indicate that the approximating accuracy to f(x,y) using Qf(x,y) and Q~f(x,y) is much more accurate than that using WLRBF method and Shepard method with the same number of scattered points.Table 3
The comparison of MSE and MME of WLRBF method and Shepard method for the same scattered points.
Number of scatter data
N
W
Influence radiusR
WLRBF
Shepard’s Method
MSE
MME
MSE
MME
100
9
0.1818
1.9381
e
-
004
0.0923
1.9297
e
-
003
0.01842
300
9
0.1179
3.6377
e
-
006
0.0125
6.8555
e
-
004
0.1840
500
9
0.0918
9.6671
e
-
007
0.0088
4.2839
e
-
004
0.1967
800
9
0.0705
2.8995
e
-
007
0.0046
2.3199
e
-
004
0.0800
1000
9
0.0654
1.6362
e
-
007
0.0033
1.5849
e
-
008
0.0645The data in Table4 indicates that the computation efficiencies of Qf(x,y) and Q~f(x,y) are much higher than those of WLRBF and Shepard methods.Table 4
The computing time of the four methods on the50×50 grid points in the unit square.
Number of scatter data
Q
f
(
x
,
y
)
Q
~
f
(
x
,
y
)
WLRBF
Shepard’s method
100
5.6
15.46
185.64
52.54
300
28.59
30.35
399.38
113.99
500
42.31
13.77
586.63
180.94
800
64.21
68.04
779.32
225.97
1000
80.28
80.66
1102.64
271.22


## 6. Conclusion

In this paper, Delaunay triangulation based on planar point set is used to obtain the piecewise bivariate Hermite interpolation function in order to approximate three-dimensional scattered data sets. When point(x,y) falls on a triangle of the triangulation, the Hermite interpolation on the triangle is used for approximate calculation. If point (x,y) falls outside the triangulation area, the node basis function value weighted sum of the points that are closer to (x,y) is used as the approximate value of f(x,y). Since the constructed interpolation function needs the first derivative value of the approximated function, however, the given scattered data set does not provide such information, our interpolation scheme uses local radial basis interpolation function to estimate the first derivative of each scattered point. Numerical experiments show that our methods have strong approximation ability to scattered data and the estimation of derivatives by local radial basis interpolation has high accuracy. Owing to no demand for solving linear system and the weight functions with local support, our methods are easily implemented and have high computational efficiency, so it is better than B-Spline least square fitting which needs solving a big enough linear system. The use of local radial basis function interpolation to estimate the derivatives is still consuming time. Therefore, how to construct a simple and high-accurate numerical differential formula is one of our future works. Meanwhile, for nonuniform distributed scattered data, how to more reasonably approximate them is another task in further work. In addition, in the process of numerical experiments, we found that the triangulation based on scattered point set has a greater impact on the test results. Hence, how to construct a triangulation more suitable for approximate schemes also needs to be considered in the future work.



---

*Source: 239703-2013-04-08.xml*