Piecewise Bivariate Hermite Interpolations for Large Sets of Scattered Data

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 (𝑥,𝑦) 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 (𝑥,𝑦) 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.

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 (𝑥,𝑦) into two types of interior point and exterior outer point. If point (𝑥,𝑦) is an interior point, then the value of the Hermite interpolation function on the triangle which point (𝑥,𝑦) lies on is regarded as the approximate value of the approximated function, while if point (𝑥,𝑦) 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 Section 2. 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 function 𝑓(𝑥,𝑦) at three vertices (𝑥_𝑖,𝑦_𝑖),𝑖=1,2,3 of a triangle 𝑇₁₂₃ (their vertices 1, 2, 3 are arranged counter to clockwise), that is, 𝑓_𝑖,(𝑓_𝑥)_𝑖,(𝑓_𝑦)_𝑖,  𝑖=1,2,3 (see Figure 1). A bivariate Hermite interpolation polynomial 𝐻₁₂₃(𝑥,𝑦) is to be constructed satisfying the following interpolation conditions: 𝐻₁₂₃𝑥_𝑖,𝑦_𝑖=𝑓_𝑖,𝜕𝐻₁₂₃𝑥𝜕𝑥^𝑖,𝑦_𝑖=𝑓_𝑥𝑖,𝜕𝐻₁₂₃𝑥𝜕𝑦^𝑖,𝑦_𝑖=𝑓_𝑦𝑖,𝑖=1,2,3. (1)

Figure 1 

Interpolation triangle 𝑇123.

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 𝑃(𝑥,𝑦) is a point on triangle 𝑇₁₂₃ and its barycentric coordinates is denoted as (𝑙₁,𝑙₂,𝑙₃). Then the relation between its cartesian coordinates and barycentric coordinates is 𝑙₁=𝑥₂𝑦₃−𝑥₃𝑦₂+𝑦₂−𝑦₃𝑥+𝑥₃−𝑥₂𝑦,𝑙2𝑆₂=𝑥₃𝑦₁−𝑥₁𝑦₃+𝑦₃−𝑦₁𝑥+𝑥₁−𝑥₃𝑦,𝑙2𝑆₃=𝑥₁𝑦₂−𝑥₂𝑦₁+𝑦₁−𝑦₂𝑥+𝑥₂−𝑥₁𝑦,𝑙2𝑆₁+𝑙₂+𝑙₃=1,𝑙₁,𝑙₂,𝑙₃≥0,𝑥=𝑥₁𝑙₁+𝑥₂𝑙₂+𝑥₃𝑙₃,𝑦=𝑦₁𝑙₁+𝑦₂𝑙₂+𝑦₃𝑙₃, (2) where 𝑆 is the area of the triangle 𝑇₁₂₃ and 𝑆=1𝑥₁𝑦₁1𝑥₂𝑦₂1𝑥₃𝑦₃. (3) Making use of transformation (2), the interpolation conditions (1) can be changed into the following form under barycentric coordinates: 𝐻₁₂₃(1,0)=𝑓₁,𝐻₁₂₃(0,1)=𝑓₂,𝐻₁₂₃(0,0)=𝑓₃,𝜕𝐻₁₂₃𝜕𝑙₁(1,0)=𝑥₁−𝑥₃𝑓_𝑥1+𝑦₁−𝑦₃𝑓_𝑦1,𝜕𝐻₁₂₃𝜕𝑙₁(0,1)=𝑥₁−𝑥₃𝑓_𝑥2+𝑦₁−𝑦₃𝑓_𝑦2,𝜕𝐻₁₂₃𝜕𝑙₁(0,0)=𝑥₁−𝑥₃𝑓_𝑥3+𝑦₁−𝑦₃𝑓_𝑦3,𝜕𝐻₁₂₃𝜕𝑙₂(1,0)=𝑥₂−𝑥₃𝑓_𝑥1+𝑦₂−𝑦₃𝑓_𝑦1,𝜕𝐻₁₂₃𝜕𝑙₂(0,1)=𝑥₂−𝑥₃𝑓_𝑥2+𝑦₂−𝑦₃𝑓_𝑦2,𝜕𝐻₁₂₃𝜕𝑙₂(0,0)=𝑥₂−𝑥₃𝑓_𝑥3+𝑦₂−𝑦₃𝑓_𝑦3. (4)

A complete bivariate polynomial of three degrees 𝑝₃(𝑥,𝑦)=∑³_𝑖+𝑗=0𝑐_𝑖𝑗𝑥^𝑖𝑦^𝑗 has ten undetermined coefficients. Under barycentric coordinate its expression is 𝑝₃𝑥,𝑦=𝑐₁𝑙³₁+𝑐₂𝑙³₂+𝑐₃𝑙³₃+𝑐₄𝑙²₁𝑙₂+𝑐₅𝑙²₁𝑙₃+𝑐₆𝑙²₂𝑙₁+𝑐₇𝑙²₂𝑙₃+𝑐₈𝑙²₃𝑙₁+𝑐₉𝑙²₃𝑙₂+𝑐₁₀𝑙₁𝑙₂𝑙₃. (5)

Since the Hermite interpolation conditions (1) or (4) only provide 9 conditions, we construct a Hermite interpolation polynomial with nine coefficients 𝐻_𝑓𝑙₁,𝑙₂=𝑐₁𝑙³₁+𝑐₂𝑙³₂+𝑐₃𝑙³₃+𝑐₄𝑙²₁𝑙₂+𝑐₅𝑙²₁𝑙₃+𝑐₆𝑙²₂𝑙₁+𝑐₇𝑙²₂𝑙₃+𝑐₈𝑙²₃𝑙₁+𝑐₉𝑙²₃𝑙₂+0.5𝑐₇+𝑐₈+𝑐₉𝑙₁𝑙₂𝑙₃, (6) which meets interpolation condition (4). Then the values of the coefficients in formula (6) can be obtained: 𝑐₁=𝑓₁,𝑐₂=𝑓₂,𝑐₃=𝑓₃,𝑐₄=3𝑓₁−𝜕𝐻₁₂₃𝜕𝑙₁(1,0)+𝜕𝐻₁₂₃𝜕𝑙₂(),𝑐1,0₅=3𝑓₁−𝜕𝐻₁₂₃𝜕𝑙₁(),𝑐1,0₆=3𝑓₂−𝜕𝐻₁₂₃𝜕𝑙₂(0,1)+𝜕𝐻₁₂₃𝜕𝑙₁(),𝑐0,1₇=3𝑓₂−𝜕𝐻₁₂₃𝜕𝑙₂(),𝑐0,1₈=3𝑓₃+𝜕𝐻₁₂₃𝜕𝑙₁(),𝑐0,0₉=3𝑓₃+𝜕𝐻₁₂₃𝜕𝑙₂(0,0). (7)

Combining (4), (6), and (7), the Hermite interpolation polynomial satisfying (4) can be written as 𝐻₁₂₃=𝑥,𝑦_3𝑖=1𝛼_𝑖𝑥,𝑦𝑓_𝑖+𝛽_𝑖𝑥,𝑦𝑓_𝑥𝑖+𝛾_𝑖𝑥,𝑦𝑓_𝑦𝑖, (8) where 𝛼₁𝑥,𝑦=𝑙³₁+3𝑙₂𝑙²₁+3𝑙₃𝑙²₁+2𝑙₁𝑙₂𝑙₃,𝛼₂𝑥,𝑦=𝑙³₂+3𝑙₁𝑙²₂+3𝑙₃𝑙²₂+2𝑙₁𝑙₂𝑙₃,𝛼₃𝑥,𝑦=𝑙³₃+3𝑙₁2𝑙²₃+3𝑙₂𝑙²₃+2𝑙₁𝑙₂𝑙₃,𝛽₁=𝑥₂−𝑥₁𝑙²₁𝑙₂+0.5𝑙₁𝑙₂𝑙₃+𝑥₃−𝑥₁𝑙²₁𝑙₃+0.5𝑙₁𝑙₂𝑙₃,𝛽₂=𝑥₁−𝑥₂𝑙²₂𝑙₁+0.5𝑙₁𝑙₂𝑙₃+𝑥₃−𝑥₂𝑙²₂𝑙₃+0.5𝑙₁𝑙₂𝑙₃,𝛽₃=𝑥₁−𝑥₃𝑙²₃𝑙₁+0.5𝑙₁𝑙₂𝑙₃+𝑥₂−𝑥₃𝑙²₃𝑙₂+0.5𝑙₁𝑙₂𝑙₃,𝛾₁=𝑦₂−𝑦₁𝑙²₁𝑙₂+0.5𝑙₁𝑙₂𝑙₃+𝑦₃−𝑦₁𝑙²₁𝑙₃+0.5𝑙₁𝑙₂𝑙₃,𝛾₂=𝑦₁−𝑦₂𝑙²₂𝑙₁+0.5𝑙₁𝑙₂𝑙₃+𝑦₃−𝑦₂𝑙²₂𝑙₃+0.5𝑙₁𝑙₂𝑙₃,𝛾₃=𝑦₁−𝑦₃𝑙²₃𝑙₁+0.5𝑙₁𝑙₂𝑙₃+𝑦₂−𝑦₃𝑙²₃𝑙₂+0.5𝑙₁𝑙₂𝑙₃. (9)

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

Given a set of scattered data {(𝑥_𝑗,𝑦_𝑗,𝑓_𝑗)}^𝑁_𝑗=1⊂𝑅³, which are assumed to be sampled from a function𝑓(𝑥,𝑦),(𝑥,𝑦)∈Ω. The partial derivatives of 𝑓(𝑥,𝑦) at the set of points {(𝑥_𝑗,𝑦_𝑗)}^𝑁_𝑗=1⊂Ω⊂𝑅² are given as (𝑓_𝑥)_𝑗,(𝑓_𝑦)_𝑗,𝑗=1,…,𝑁. Taking the 𝑁 scattered points as nodes, a triangulation mesh, 𝑇, is constructed in domain Ω using Delaunay triangulation method. If point (𝑥,𝑦)∈Ω⊂𝑅² lies on a triangle of 𝑇, it is called as an interior point, otherwise, an exterior point. For example, points 𝑄₁,𝑄₂ in Figure 2 are interior points, while point 𝑄₃ is an exterior point. Next we will give the approximate value 𝑄_𝑓(𝑥,𝑦) of 𝑓(𝑥,𝑦) according to the category of point (𝑥,𝑦).

Figure 2 

The exhibition of interior and exterior points.

3.1. Interior Point

Suppose that point (𝑥,𝑦) lies on a triangle 𝑇_𝑗𝑘𝑙⊂𝑇 whose vertices are (𝑥_𝑗,𝑦_𝑗),(𝑥_𝑘,𝑦_𝑘),(𝑥_𝑙,𝑦_𝑙). We now construct a Hermite interpolation polynomial as (8) on triangle 𝑇_𝑗𝑘𝑙 satisfying the following conditions: 𝐻_𝑗𝑘𝑙𝑥_𝑗,𝑦_𝑗=𝑓_𝑗,𝜕𝐻_𝑗𝑘𝑙𝑥𝜕𝑥^𝑗,𝑦_𝑗=𝑓_𝑥𝑗,𝜕𝐻_𝑗𝑘𝑙𝑥𝜕𝑦^𝑗,𝑦_𝑗=𝑓_𝑦𝑗,𝐻_𝑗𝑘𝑙𝑥_𝑘,𝑦_𝑘=𝑓_𝑘,𝜕𝐻_𝑗𝑘𝑙𝑥𝜕𝑥^𝑘,𝑦_𝑘=𝑓_𝑥𝑘,𝜕𝐻_𝑗𝑘𝑙𝑥𝜕𝑦^𝑘,𝑦_𝑘=𝑓_𝑦𝑘,𝐻_𝑗𝑘𝑙𝑥_𝑙,𝑦_𝑙=𝑓_𝑙,𝜕𝐻_𝑗𝑘𝑙𝑥𝜕𝑥^𝑙,𝑦_𝑙=𝑓_𝑥𝑙,𝜕𝐻_𝑗𝑘𝑙𝑥𝜕𝑦^𝑙,𝑦_𝑙=𝑓_𝑦𝑙. (10) Then let 𝑄_𝑓(𝑥,𝑦)=𝐻_𝑗𝑘𝑙(𝑥,𝑦), so that 𝑄_𝑓(𝑥,𝑦)≈𝑓(𝑥,𝑦).

3.2. Exterior Point

When point (𝑥,𝑦) doesnot belong to any triangle within 𝑇, we construct the following form of Hermite interpolation: 𝑄_𝑓𝑥,𝑦=^𝑁_𝑗=1𝑊_𝑗𝑥,𝑦𝐻_𝑓,𝑗𝑥,𝑦, (11) satisfying interpolation conditions 𝑄_𝑓𝑥_𝑗,𝑦_𝑗=𝑓_𝑗,𝜕𝑄_𝑓𝑥𝜕𝑥^𝑗,𝑦_𝑗=𝑓_𝑥𝑗,𝜕𝑄_𝑓𝑥𝜕𝑦^𝑗,𝑦_𝑗=𝑓_𝑦𝑗,𝑗=1,2,…,𝑁. (12)

In this paper, we take 𝐻_𝑓,𝑗(𝑥,𝑦) as a Hermite interpolation polynomial with the same form of (8) on triangle 𝑇_𝑗𝑘𝑙 with the least area among all triangles with vertex (𝑥_𝑗,𝑦_𝑗), satisfying interpolation conditions (10). 𝐻_𝑓,𝑗(𝑥,𝑦) is called as node basis function and function 𝑊_𝑗(𝑥,𝑦) is called as weighted function. In order to satisfy interpolation conditions (12), these weighted functions are required for satisfying conditions 𝑊_𝑗𝑥_𝑖,𝑦_𝑖=𝛿_𝑖𝑗,𝑖,𝑗=1,…,𝑁,_𝑁𝑗=1𝑊_𝑗𝑥,𝑦≡1,𝑊_𝑗𝑥,𝑦≥0,𝑗=1,…,𝑁,𝜕𝑊_𝑗𝑥𝜕𝑥^𝑖,𝑦_𝑖=𝜕𝑊_𝑗𝑥𝜕𝑦^𝑖,𝑦_𝑖=0,𝑖,𝑗=1,…,𝑁. (13)(14)(15)

In the paper we take the weighted functions in (11) as those in the literature [16]: 𝑊_𝑗𝑥,𝑦=1/𝜌²_𝑗∑^𝑁_𝑘=11/𝜌²_𝑘, (16) where 1𝜌_𝑘=𝑅_𝑘−𝑑_𝑘+𝑅_𝑘𝑑_𝑘,𝑅_𝑘−𝑑_𝑘+=𝑅_𝑘−𝑑_𝑘𝑅_𝑘−𝑑_𝑘≥0,0𝑅_𝑘−𝑑_𝑘𝑑<0,_𝑘=𝑑_𝑘𝑥,𝑦=𝑥−𝑥_𝑘²+𝑦−𝑦_𝑘². (17)

The value of 𝑅_𝑘 in (17) can be selected by user properly. It can be seen that when (𝑥,𝑦)→(𝑥_𝑗,𝑦_𝑗), (1/𝜌²_𝑗)→+∝, thereby 𝑊_𝑗(𝑥,𝑦)→1; when (𝑥,𝑦)→(𝑥_𝑘,𝑦_𝑘),𝑘≠𝑗 and (1/𝜌²_𝑘)→+∝, thereby 𝑊_𝑗(𝑥,𝑦)→0. Thus 𝑊_𝑗 satisfies condition (13). 𝑊_𝑗 in formula (16) satisfies condition (14) obviously, and the proof of which can be found in literature [16]. If 1/𝜌_𝑗 is determined by (17), then when 𝑑_𝑗(𝑥,𝑦)>𝑅_𝑗,𝑊_𝑗(𝑥,𝑦)=0. So in case of 𝑅_𝑗 being appropriately selected (𝑅_𝑗 is also called as influence radius of node basis function), there are only several nonzero terms in sum ∑^𝑁_𝑘=1𝑊_𝑘(𝑥,𝑦). 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 𝑅. The following computational formulation of uniform influence radius 𝑅 in (18) is from [16]: 𝐷𝑅=2𝑁_𝑤𝑁,𝐷=max_𝑖,𝑗𝑑_𝑖𝑥_𝑗,𝑦_𝑗. (18)𝑁_𝑊 is the number of points used inside a circle of radius 𝑅.

3.3. Bivariate Piecewise Hermite Interpolation

Summing up Sections 3.1 and 3.2, a bivariate piecewise Hermite interpolation function is constructed on Ω to approximate the function 𝑓(𝑥,𝑦)𝑄_𝑓=𝐻𝑥,𝑦^𝑗𝑘𝑙𝑥,𝑦,if𝑥,𝑦∈Γ,_𝑁𝑗∑₌₁𝑊_𝑗𝑥,𝑦𝐻_𝑓,𝑗𝑥,𝑦,if𝑥,𝑦∈Ω−Γ, (19) where 𝑇_𝑗𝑘𝑙 is a triangle within 𝑇 with node  (𝑥_𝑗,𝑦_𝑗),(𝑥_𝑘,𝑦_𝑘),(𝑥_𝑙,𝑦_𝑙) as its vertices and 𝐻_𝑗𝑘𝑙 is a Hermite interpolation polynomial on triangle 𝑇_𝑗𝑘𝑙.

We now summarize the description of the constructed interpolation function, written as Algorithm 1.

Algorithm 1. Consider the following.
Step  1. Generate a triangle mesh T in domain Ω using Delaunay triangulation method based on given scattered point set {(𝑥_𝑗,𝑦_𝑗)}^𝑁_𝑗=1⊂Ω.
Step  2. Judge the category of point (𝑥,𝑦): interior point or exterior point.
Step  3. If point (𝑥,𝑦) is an interior point, then find out the triangle 𝑇_𝑗𝑘𝑙 which includes point (𝑥,𝑦), compute the value of the interpolation function 𝐻_𝑗𝑘𝑙 at point (𝑥,𝑦) using formulas (2), (8), and (10), and take the value as the value of 𝑄_𝑓(𝑥,𝑦).
Step  4. If point (𝑥,𝑦) is an exterior point, then compute 𝐷=max_𝑖,𝑗𝑑_𝑖(𝑥_𝑗,𝑦_𝑗) and select 𝑁_𝑊 to define 𝑅=(𝐷/2)𝑁_𝑊/𝑁. The default value of 𝑁_𝑊 is set to 9, and it responds to an uniform radius. Find out the points (𝑥_𝑗𝑙,𝑦_𝑗𝑙),𝑙=1,…,𝑘 from scattered point set {(𝑥_𝑗,𝑦_𝑗)}^𝑁_𝑗=1, which belongs to the circle with radius 𝑅 and center (𝑥,𝑦), compute the values 𝑊_𝑗𝑙(𝑥,𝑦),𝑙=1,…,𝑘 and 𝐻_{𝑓,𝑗𝑙}(𝑥,𝑦),𝑙=1,…,𝑘, give the sum ∑^𝑘_𝑙=1𝑊_𝑗𝑙(𝑥,𝑦)𝐻_{𝑓,𝑗𝑙}(𝑥,𝑦), and assign it to 𝑄_𝑓(𝑥,𝑦).

4. The Estimation of Partial Derivative

In the process of using bivariate piecewise Hermite interpolation 𝑄_𝑓, it demands the first-order derivatives of the interpolated function 𝑓(𝑥,𝑦) at every point, but the scattered data {(𝑥_𝑗,𝑦_𝑗,𝑓_𝑗)}^𝑁_𝑗=1  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 𝑓(𝑥,𝑦) and replace this exact derivative in (19) with the approximate derivative.

The Delaunay triangulation 𝑇 is a convex subdomain of Ω. Every scattered point 𝑋_𝑗=(𝑥_𝑗,𝑦_𝑗) is either a boundary point or an interior node point of 𝑇. We put the vertices of all triangles with point (𝑥_𝑗,𝑦_𝑗) as into a set, written as 𝐴_𝑗=(𝑥_𝑗,𝑦_𝑗)∪{(𝑥_𝑗𝑙,𝑦_𝑗𝑙)}^𝑚_𝑙=1, see Figure 3.

(a)

(b)

Figure 3 

The local radial basis interpolation point set 𝐴𝑗 around point 𝑋𝑗, (a) 𝑋𝑗 is a boundary point of 𝑇 and (b) 𝑋𝑗 is an interior node point of 𝑇.

We take the point set 𝐴_𝑗 as local radial basis function interpolation point set and multiquadric radial basis functions 𝜙_𝑗=𝑥−𝑥_𝑗²+𝑦−𝑦_𝑗+𝑐²,𝜙_𝑗1=𝑥−𝑥_𝑗1²+𝑦−𝑦_𝑗1²+𝑐²,⋮𝜙_𝑗𝑚=𝑥−𝑥_𝑗𝑚²+𝑦−𝑦_𝑗𝑚²+𝑐² (20) as local interpolation basis functions.

By solving the linear system 𝐺𝐶=𝐹, 𝜙𝐺=^𝑗𝑥_𝑗,𝑦_𝑗𝜙_𝑗1𝑥_𝑗,𝑦_𝑗…𝜙_𝑗𝑚𝑥_𝑗,𝑦_𝑗𝜙_𝑗𝑥_𝑗1,𝑦_𝑗1𝜙_𝑗1𝑥_𝑗1,𝑦_𝑗1…𝜙_𝑗𝑚𝑥_𝑗1,𝑦_𝑗1𝜙…⋯……_𝑗𝑥_𝑗𝑚,𝑦_𝑗𝑚𝜙_𝑗1𝑥_𝑗𝑚,𝑦_𝑗𝑚…𝜙_𝑗𝑚𝑥_𝑗𝑚,𝑦_𝑗𝑚,𝐶^𝑇=𝑐₀𝑐₁𝑐₂…𝑐_𝑚,𝐹^𝑇=𝑓_𝑗𝑓_𝑗1𝑓_𝑗2…𝑓_𝑗𝑚, (21) we obtained a local radial basis interpolation function 𝑝_𝑗(𝑥,𝑦)=𝑐₀𝜙_𝑗(𝑥,𝑦)+𝑐₁𝜙_𝑗1(𝑥,𝑦)+𝑐₂𝜙_𝑗2(𝑥,𝑦)+⋯+𝑐_𝑚𝜙_𝑗𝑚(𝑥,𝑦), and then compute the first order derivatives of 𝑝_𝑗(𝑥,𝑦) at 𝑋_𝑗=(𝑥_𝑗,𝑦_𝑗), that is, 𝜕𝑝_𝑗𝑥𝜕𝑥^𝑗,𝑦_𝑗=𝑐₁𝑥_𝑗−𝑥_𝑗1𝜙_𝑗1𝑥_𝑗,𝑦_𝑗+𝑐₂𝑥_𝑗−𝑥_𝑗2𝜙_𝑗2𝑥_𝑗,𝑦_𝑗𝑐+⋯+^𝑚𝑥_𝑗−𝑥_𝑗𝑚𝜙_𝑗𝑚𝑥_𝑗,𝑦_𝑗,𝜕𝑝_𝑗𝑥𝜕𝑦^𝑗,𝑦_𝑗=𝑐₁𝑦_𝑗−𝑦_𝑗1𝜙_𝑗1𝑥_𝑗,𝑦_𝑗+𝑐₂𝑦_𝑗−𝑦_𝑗2𝜙_𝑗2𝑥_𝑗,𝑦_𝑗𝑐+⋯+^𝑚𝑦_𝑗−𝑦_𝑗𝑚𝜙_𝑗𝑚𝑥_𝑗,𝑦_𝑗. (22)

The bivariate piecewise interpolation function 𝑄_𝑓 using the previous approximate derivatives is written as 𝑄^∼_𝑓𝑥,𝑦=^∼𝐻_𝑗𝑘𝑙𝑥,𝑦,if𝑥,𝑦∈Γ,_𝑁𝑗∑₌₁𝑊_𝑗𝑥,𝑦𝐻^∼_𝑓,𝑗𝑥,𝑦,if𝑥,𝑦∈Ω−Γ. (23) Algorithm 2 summarizes the computational process of interpolation 𝑄^∼_𝑓(𝑥,𝑦).

Algorithm 2. Consider the following.
Step  1. Generate a triangle mesh 𝑇 in Ω using Delaunay triangulation method based on the given scattered point set {(𝑥_𝑗,𝑦_𝑗)}^𝑁_𝑗=1⊂Ω.
Step  2. Find out the local radial basis interpolation point set 𝐴_𝑗 at every point (𝑥_𝑗,𝑦_𝑗), and compute the approximate derivatives (𝜕𝑝_𝑗/𝜕𝑥)(𝑥_𝑗,𝑦_𝑗),(𝜕𝑝_𝑗/𝜕𝑦)(𝑥_𝑗,𝑦_𝑗) at every point (𝑥_𝑗,𝑦_𝑗) using interpolation formulation (21) and (22).
Step  3. Judge the category of point (𝑥,𝑦): interior point or exterior point.
Step  4. If point (𝑥,𝑦) is an interior point, then find out the triangle 𝑇_𝑗𝑘𝑙 which point (𝑥,𝑦) falls in, compute the value of the interpolation function 𝐻_𝑗𝑘𝑙 at point (𝑥,𝑦) using formulas (2), (8), and (10), and assign the value to the function 𝑄^∼_𝑓(𝑥,𝑦).
Step  5. If point (𝑥,𝑦) is an exterior point, then compute 𝐷=max_𝑖,𝑗𝑑_𝑖(𝑥_𝑗,𝑦_𝑗) and select 𝑁_𝑊 to define 𝑅=(𝐷/2)𝑁_𝑊/𝑁. The default value of 𝑁_𝑊 is set to 9, it responds to an uniform radius. Find out the points (𝑥_𝑗𝑙,𝑦_𝑗𝑙),𝑙=1,…,𝑘 from scattered point set {(𝑥_𝑗,𝑦_𝑗)}^𝑁_𝑗=1, which belongs to the circle with radius 𝑅 and center (𝑥,𝑦), compute the values 𝑊_𝑗𝑙(𝑥,𝑦),𝑙=1,…,𝑘 and 𝐻^∼_{𝑓,𝑗𝑙}(𝑥,𝑦),𝑙=1,…,𝑘, give the sum ∑^𝑘_𝑙=1𝑊_𝑗𝑙(𝑥,𝑦)𝐻^∼_{𝑓,𝑗𝑙}(𝑥,𝑦), and assign it to 𝑄^∼_𝑓(𝑥,𝑦).

5. Numerical Test

In this section, 𝑄_𝑓(𝑥,𝑦) and 𝑄^∼_𝑓(𝑥,𝑦) are used to approximate Franke function 𝑓−()𝑥,𝑦=0.75exp9𝑥−2²+9𝑦−2²4−()×0.75exp9𝑥+1²−109𝑦+1−()10−0.2exp9𝑥−4²−9𝑦−7²−()+0.5exp9𝑥−7²+9𝑦−3²4,0≤𝑥,𝑦≤1. (24) Firstly, a scattered data set {(𝑥_𝑗,𝑦_𝑗,𝑓_𝑗)}^𝑁_𝑗=1 is sampled from 𝑓(𝑥,𝑦); then, the interpolation functions 𝑄_𝑓(𝑥,𝑦) and 𝑄^∼_𝑓(𝑥,𝑦) 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 𝑄_𝑓(𝑥,𝑦) and 𝑄^∼_𝑓(𝑥,𝑦) 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.

Figure 4 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 𝑄_𝑓(𝑥,𝑦), 𝑄^∼_𝑓(𝑥,𝑦), WLRBF, and Shepard which are also based on the 300 scattered points shown in Figure 4. It can be seen that 𝑄_𝑓(𝑥,𝑦), 𝑄^∼_𝑓(𝑥,𝑦), 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 𝑄_𝑓(𝑥,𝑦) is much better than WLRBF method. Comparing Figures 12 and 13 for the interior points, 𝑄^∼_𝑓(𝑥,𝑦) is obviously better than WLRBF method. Besides, with the number of scattered data increasing, it is found that approximation capability of 𝑄^∼_𝑓(𝑥,𝑦) is closer to that of 𝑄_𝑓(𝑥,𝑦). In case of large set of scattered points, 𝑄^∼_𝑓(𝑥,𝑦) 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 Figure 4.

Figure 6 

The graph of Franke function.

Figure 7 

The graph of 𝑄𝑓(𝑥,𝑦) based on Figure 4.

Figure 8 

The graph of 𝑄∼𝑓(𝑥,𝑦) based on Figure 4.

Figure 9 

The graph of WLRBF based on Figure 4.

Figure 10 

The graph of Shepard interpolation based on Figure 4.

Figure 11 

The error graph of 𝑄𝑓(𝑥,𝑦).

Figure 12 

The error graph of 𝑄∼𝑓(𝑥,𝑦).

Figure 13 

The error graph of WLRBF.

Figure 14 

The error graph of Shepard interpolation.

It can be seen from Table 1 that the approximation accuracy of 𝑄_𝑓(𝑥,𝑦) is higher than that of 𝑄^∼_𝑓(𝑥,𝑦) 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 𝑓(𝑥,𝑦) of 𝑄_𝑓(𝑥,𝑦) and 𝑄^∼_𝑓(𝑥,𝑦) 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, 𝑄^∼_𝑓(𝑥,𝑦) can be used to approximate 𝑓(𝑥,𝑦) when the partial derivative information of the scattered data is unknown.

Table 1 

The comparison of MSE and MME of 𝑄_𝑓(𝑥, 𝑦) and 𝑄^∼_𝑓(𝑥, 𝑦) for the same scattered points.

It can be seen from Table 2 that the errors of approximation functions 𝑄_𝑓(𝑥,𝑦) and 𝑄^∼_𝑓(𝑥,𝑦) increase with the radius of influence decreasing.

Table 2 

The influence of influence radius on 𝑄_𝑓(𝑥,𝑦) and 𝑄^∼_𝑓(𝑥,𝑦).

Tables 1 and 3 indicate that the approximating accuracy to 𝑓(𝑥,𝑦) using 𝑄_𝑓(𝑥,𝑦) and 𝑄^∼_𝑓(𝑥,𝑦) 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.

The data in Table 4 indicates that the computation efficiencies of 𝑄_𝑓(𝑥,𝑦) and 𝑄^∼_𝑓(𝑥,𝑦) are much higher than those of WLRBF and Shepard methods.

Table 4 

The computing time of the four methods on the 50×50 grid points in the unit square.

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 (𝑥,𝑦) falls on a triangle of the triangulation, the Hermite interpolation on the triangle is used for approximate calculation. If point (𝑥,𝑦) falls outside the triangulation area, the node basis function value weighted sum of the points that are closer to (𝑥,𝑦) is used as the approximate value of 𝑓(𝑥,𝑦). 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.

Acknowledgment

Supported by the National Natural Science Foundation of China (no. 11271041), National Natural Science Key Foundation of China (no. 10931004) and ISTCP of China grant (no. 2010DFR00700).