Let $S$ denote the unit sphere in ${ℝ}^3$ centered at the origin

with $\parallel \cdot \parallel$ representing the Euclidean norm in ${ℝ}^3$. Consider $𝒳\subset S$ a set of sample points ${𝐱}_i,i=1,\mathrm{\dots },n$, lying on the unit sphere $S$, for which the values $f_i,i=1,\mathrm{\dots },n$, of a function $f:S\to ℝ$ are known. For $𝐱\in S$, the spherical Shepard operator is given by

where the cardinal basis functions $A_{i,\mu }$ are defined in the barycentric form as

for $\mu \in {ℝ}_{+}$ a control parameter. Here, the metric $g$ that is used represents the geodesic distance between $𝐱$ and ${𝐱}_i$ and it is defined as

with ${}_{}{}^{\prime \prime }\cdot _{}^{\prime \prime }$ denoting the Euclidean inner product of the two vectors.

The cardinal basis functions $A_{i,\mu }$ have the property that their partial derivatives vanish at the nodes, up to a certain degree.

The original Shepard method has been introduced in [37]. Over time, numerous improvements were made to overcome the drawbacks of the initial method. For instance, the function values $f_i$ have been replaced by the values of another operator that provide better approximation results, and a higher degree of exactness (see, e.g, [2, 4, 5, 6, 7, 23, 26]). Other improvements include the consideration of other forms of basis functions (see, e.g., [27, 28, 33, 34, 41]). Here we restrict our attention to the case of the spherical Shepard-Bernoulli operator. The Bernoulli operator or its combination with the Shepard interpolant has been extensively studied in the univariate and bivariate cases (see, e.g., [2, 5, 16, 17, 23]). The problem of interpolating large sets of scattered data on the sphere using Shepard-like methods has started to be studied recently by several authors (see, e.g, [1, 8, 10, 11, 12, 13, 25, 39, 40]).

We construct the new spherical operators based on two types of Shepard-type basis functions. The combination between these interpolants and the Bernoulli operators will be performed on triangular domains, after computing the Delaunay triangulation of the sphere $S$. Numerical tests on various functions and sets of nodes demonstrate the efficiency of the proposed methods. Furthermore, we showcase two practical applications for modeling real-life phenomena, specifically topographic and temperature prediction problems.

The triangular Shepard operator has been introduced by Little in [30] and further developed, for instance, in [14, 22, 24]. In [14], Cavoretto, De Rossi, Dell’Accio, and Di Tommaso discussed the implementation of the 2D-triangular Shepard method. They proposed an algorithm that successfully performs a compact triangulation of the domain, i.e., one that allows triangles to overlap or be disjoint. This is done based on a searching technique for determining the closest neighboring points in the interpolation scheme and on the construction of a block-based partitioning structure. In [15], they discussed the 3D analog of the previously mentioned method. The algorithm used to obtain a fast computation of the tetrahedral Shepard method involves finding a searching procedure that efficiently returns the closest neighbors of each interpolation node. Similar to the previous case, this is done by constructing a structure that partitions the domain in cubic blocks and by using a limited number of these blocks to determine the nearest neighbors. In both cases, they achieved an improvement in terms of computational efficiency.

Another type of triangulation, which requires additional conditions, is the Delaunay triangulation. The problem of data interpolation on the sphere, with the aid of triangulation methods, has been addressed and solved, for example, in [31, 32]. Renka proposed in [35] an algorithm for constructing the Delaunay triangulation of a sphere, which we will follow in the sequel. We recall some definitions from [35].

In this section, we recall some notions regarding the triangular Bernoulli operators in ${ℝ}^2$, that were studied in [23]. We shall consider in this part that $𝐱=(x,y)\in {ℝ}^2$.

First, let ${𝐱}_i,{𝐱}_j,{𝐱}_k$ be three nodes specified in counterclockwise order, forming the vertices of a triangle $T(i)$, with the referring vertex ${𝐱}_i$. To the ease of notations, we shall consider ${𝐱}_i:={𝐱}_1$, ${𝐱}_j:={𝐱}_2$, ${𝐱}_k:={𝐱}_3$ inside the triangle $T(i)$.

For a point $𝐱$, let ${\lambda }_1,{\lambda }_2,{\lambda }_3$ be its barycentric coordinates relative to $T(i)$. They are defined as

where $𝒜(𝐱,𝐲,𝐳)$ represents the signed area of the triangle with vertices $𝐱,𝐲,𝐳$.

The directional derivatives of a differentiable function $f$ along the sides of the triangle $T(i)$ are given as

The composition of the directional derivatives is expressed as

for $({\gamma }_1,{\gamma }_2)\in \mathbb{N}\times \mathbb{N}$.

The Bernoulli polynomials $B_n(x)$, $n\in \mathbb{N}$, $x\in ℝ$, can be defined recursively as

this being one of the many possible definitions, as it is shown in [19]. For $x=0$, one obtains the Bernoulli numbers $B_n$, i.e., $B_n=B_n(0)$.

There are also some other related polynomial expansions to Bernoulli polynomials in the triangle, such as expansions using Lidstone and complementary Lidstone polynomials, as described in [3], [18], [20], [21]. These kind of operators we intend to generalize for scattered data sets on the sphere in some future studies.

To obtain the new combined Shepard operators of Bernoulli type, we will use the spherical coordinates $(\varphi ,\theta )$ corresponding to a point $𝐱\in S$, given in cartesian coordinates $(x,y,z)$.

If $(\varphi ,\theta )$ represent the latitude and the longitude, respectively, with $\varphi \in [0,\pi ],\theta \in [0,2\pi ]$, then, for $𝐱=(x,y,z)\in S$, we have the following relation between the two coordinates representations:

Similarly to the planar case (3.2), one can obtain the directional derivatives of a function $f$ with respect to $\varphi$ and $\theta$ (see, e.g., [31]), so, $B_m^{T(i)}[\tilde{f}](\varphi ,\theta )$ can be written similarly to (3.1), considering $\tilde{f}(\varphi ,\theta )=f(x,y,z)$.

As mentioned in [31], the drawback that appears when using latitude and longitude derivatives is the fact that at the poles they are not well defined. A solution that was proposed to overcome this drawback is to slightly move the poles when one of them is a data point.

Moreover, in this case, the barycentric coordinates (3.1) are computed based on the signed area of a spherical triangle of vertices $𝐱,𝐲,𝐳\in S$, that is

with ${}_{}{}^{\prime \prime }\cdot _{}^{\prime \prime }$ denoting the Euclidean inner product and ${}_{}{}^{\prime \prime }\times _{}^{\prime \prime }$ the vector cross product.

Let $T$ be the Delaunay triangulation of a set $𝒳$ containing $n$ data samples ${𝐱}_i\in S,i=1,\mathrm{\dots },n,$ and ${𝒯}_i$ the set of all triangles that have a vertex in ${𝐱}_i$, for each $i=1,\mathrm{\dots },n$. Choosing the representative triangle $T(i)\subset {𝒯}_i$, on which the operator $B_m^{T(i)}[\tilde{f}]$ is constructed, such that the approximation error is minimum, we introduce the spherical Shepard-Bernoulli operator, defined as

with $B_m^{T(i)}[\tilde{f}](\varphi ,\theta )$ given in (3.1), $A_{i,\mu }$ in (1.2), and $\tilde{f}(\varphi ,\theta )=f(x,y,z)$.

We obtain the following result.

In the sequel, we consider another approach, proposed in [39] and [40], that consists of constructing a Bernoulli operator on each triangle from the triangulation $T$ of $S$. As mentioned in [39], the triangulation can be general, with overlapping or disjoint triangles, but here we restrict our attention to the case of Delaunay triangulation $T$. Let $N$ be the number of triangles that form the triangulation. In this part, ${𝒯}_i,i=1,\mathrm{\dots },N,$ will denote a triangle from the triangulation, with vertices ${𝐱}_{i_1},{𝐱}_{i_2},{𝐱}_{i_3}$, such that each node ${𝐱}_i,i=1,\mathrm{\dots },n,$ is a vertex of at least one triangle, i.e.,

Each basis function ${\mathrm{\Phi }}_{i,\mu }$ corresponding to the triangle ${𝒯}_i$, $i=1,\mathrm{\dots },N,$ will have the form [39]

for $\mu \in {ℝ}_{+}$ a control parameter.

We have that ${\mathrm{\Phi }}_{i,\mu },$ $i=1,\mathrm{\dots },N,$ are non-negative, satisfy the cardinality relations and form a partition of unity (see, e.g., [39]), i.e.,

We can now define the modified spherical Shepard-Bernoulli operator as follows,

with $B_m^{{𝒯}_i}[\tilde{f}](\varphi ,\theta )$ given in (3.1), ${\mathrm{\Phi }}_{i,\mu }$ in (4.2), and $\tilde{f}(\varphi ,\theta )=f(x,y,z)$.

Furthermore, we list some properties of $S_{B_m}^2$.

To illustrate the theoretical results, we compute the approximation errors for both interpolants introduced, $S_{B_m}^1$ and $S_{B_m}^2$, considering the orders $m=1$ and $m=2$ for the Bernoulli operators. We construct them considering two cases for choosing the nodes and the points. We use the following test functions (see, e.g., [27], [28], [33], [34])

For the first case, we evaluate the operators on 1000 random points on the unit sphere $S$, using sets of 500, 1000, and 2000 spiral nodes, respectively. These nodes are constructed based on a method proposed in [36]. To describe them, we use their spherical coordinates

with

Tables 1, 2 and 3 present the maximum approximation errors $e_{max}$, the mean approximation errors $e_{mean}$ and the root mean squared errors $e_{rms}$ of $S_{B_m}^1[f_i]$ and $S_{B_m}^2[f_i],$ $i=1,\mathrm{\dots },5$, for $m=1,2$. We also post the computational times (CPU) for the methods introduced above. Comparing the results, we observe that the modified forms $S_{B_m}^2$, for $m=1,2$, produce better approximation results than the classical operator, and they also have reduced computational costs.

In the second case, we perform the evaluation on $600$ spiral points, using sets of $100$, $500$, and $1000$ Halton nodes, respectively. The latter are spherical pseudorandom points, generated with the method proposed in [38]. We recall their definition from [38].

With the expansion using a prime base $p$ of a natural number $k$

one can define the function

obtaining in the case $p=2$, the sequence ${\alpha }_2(k),$ for $k=0,1,2,\mathrm{\dots }$, that is called the Van der Corput sequence.

To generate spherical Halton points, two $p-$adic Van der Corput sequences with different prime bases should be used to construct the following mappings

As mentioned in [38], the first mapping is a linear scaling to a cylindrical domain, $\varphi \in [0,2\pi ),t\in [-1,1]$, while the second one is a $z-$preserving radial projection from the unit cylinder $C=\{𝐱=(x,y,z):x^2+y^2=1,|z|\le 1\}$ to the unit sphere $S$.