Discrete spherical means of directional derivatives and Veronese maps

# Discrete spherical means of directional derivatives and Veronese maps

Alexander Belyaev  Boris Khesin  Serge Tabachnikov  Electrical, Electronic & Computer Engineering, Heriot-Watt University, Edinburgh, EH14 4AS, UK Department of Mathematics, University of Toronto, Toronto, ON M5S 2E4, Canada  Department of Mathematics, Pennsylvania State University, University Park, PA 16801, USA
###### Abstract

We describe and study geometric properties of discrete circular and spherical means of directional derivatives of functions, as well as discrete approximations of higher order differential operators. For an arbitrary dimension we present a general construction for obtaining discrete spherical means of directional derivatives. The construction is based on using the Minkowski’s existence theorem and Veronese maps. Approximating the directional derivatives by appropriate finite differences allows one to obtain finite difference operators with good rotation invariance properties. In particular, we use discrete circular and spherical means to derive discrete approximations of various linear and nonlinear first- and second-order differential operators, including discrete Laplacians. A practical potential of our approach is demonstrated by considering applications to nonlinear filtering of digital images and surface curvature estimation.

## Introduction

Mean value properties of functions play crucial role in analysis of many partial differential equations Friedman and Littman (1962), numerical interpolation and integration Stroud (1971), computer tomography Natterer (2001), and many other areas of mathematics and engineering. In this paper we undertake a study of discrete circular and spherical means of directional derivatives of functions in view of their applications to the finite difference methods. We present a general construction for obtaining discrete spherical means of directional derivatives, analyze general properties of such means, and use them to derive finite difference approximations with good rotation-invariant properties for linear and nonlinear first- and second-order partial differential operators. The two principal components of our approach are based on the use of the Minkowski’s existence theorem Schneider (1993) and Veronese maps Harris (1992).

For the sake of simplicity, we focus on finite difference approximations of very basic differential operators, the gradient and Laplacian. We start from an elementary analysis of discrete circular means of derivatives and demonstrate how such discrete means can be used for practical estimation of surface curvatures (Section 1). Then, in any dimension, we obtain Minkowski-type formulas for general and regular grids by using the Veronese map (Section 2). Finally, we analyze rotation invariance properties of discrete 2D gradient, Laplacian, and related nonlinear operators defined on a square grid (Section 3).

Various finite difference approximations of the gradient and Laplacian are widely used in fluid mechanics, electromagnetics, finance modeling, image processing, and many other areas. Our interest in constructing reliable discrete approximations of these operators is threefold. Firstly, reliable and consistent approximations of the gradient and Laplacian are needed in various diffusion-type hydrodynamical models, such as the convection-diffusion equation

 ut+v⋅∇u=ϵΔu,

where the density is convected by a velocity field and dissipates at different time scales. Its attractors exhibit very peculiar behavior as and (see Arnold and Khesin (1998)) and we hope that the approximation scheme developed in this paper could complement the analytic treatment of the attractors. The Laplacian case is also particularly useful in computations related to the fast dynamo problem, where the proposed approximation formulas may simplify computations of the magnetic diffusion for iterations of a seed magnetic field, see Childress (1992).

Secondly, linear and nonlinear diffusion-type equations are widely used in modern image processing for a variety of tasks including edge detection, image restoration, and image decomposition into structure and and texture components Weickert (1998); Aubert and Kornprobst (2006). As demonstrated in Weickert and Scharr (2002), the use of discrete nonlinear diffusion operators with good rotation invariance properties is highly beneficial for such tasks. In addition, nonlinear image diffusion provides us with a useful testbed for the theoretical considerations below. Thirdly, accurate and reliable approximations of the surface gradient and Laplacian are key ingredients in many shape analysis studies. Below we discuss applications of the approximation formulas both in image processing and geometry. Finally we note that the approach via the Veronese map described in the present paper can also be used for higher order differential operators, such as bi-Laplacian.

## 1 Circular means of derivatives and their applications

### Continuous and discrete circular means

A simple way to demonstrate the rotational invariance of the 2D Laplacian

 Δ≡∂2∂x2+∂2∂y2

consists of representing it as the circular mean (respectively, the spherical mean in 3D) of the second-order directional derivatives

 12Δf=1π∫π0∂2f∂{e}2φdφ,where {e}φ=(cosφ,sinφ),l∂∂{e}φ=cosφ∂∂x+sinφ∂∂y. (1)

The rotational invariance of the squared gradient can be demonstrated in a similar way as

 12|∇f|2=1π∫π0(∂f∂{e}φ)2dφ. (2)

More generally, the mean value representation for a directional derivative of can be obtained from

 12[a∂f∂x+b∂f∂y]=1π∫π0[acosφ+bsinφ]∂f∂{e}φdφ, (3)

where and are constants. Notice that (3) with and implies (1) and setting and yields the equality (2).

Let us now derive discrete counterparts of the above mean value representations (1), (2), and (3). Consider a set of directions , . Then a weighted sum of the second-order directional derivatives can be expressed as

 ∑wk∂2f∂{e}2φk=∂2f∂x2∑wkcos2φk+∂2f∂y2∑wksin2φkMMM+2∂2f∂x∂y∑wkcosφksinφk=12Δf∑wkMMM+12(∂2f∂x2−∂2f∂y2)Re(∑wke2iφk)+∂2f∂x∂yIm(∑wke2iφk). (4)

Finding weights such that

 ∑wke2iφk=0 (5)

yields, after a simple normalization a discrete mean value representation for :

 12Δf=∑wk∂2f∂{e}2φk, (6)

and similar representations of the squared gradient

 12|∇f|2=∑wk(∂f∂{e}φk)2 (7)

and directional derivative

 12[a∂f∂x+b∂f∂y]=∑wk[acosφk+bsinφk]∂f∂{e}φk, (8)

respectively. Obviously the same set of weights can be used to obtain a discrete mean value representation for a quasi-Laplacian

 12∇⋅(a(x,y)∇f)=∑wk∂∂{e}φk(a(x,y)∂f∂{e}φk). (9)

Note that the appearance of double angles in the formula (5) exactly corresponds to the Veronese map , considered in Example 2.4 of the next section.

### Discrete mean value formulas and Minkowski problem for polygons

The key problem of finding normalized weights in the relation (see (5) above) has a simple geometric meaning and is directly related to the Minkowski existence theorem (also called Minkowski’s Problem, see e.g. Schneider (1993)) for polyhedra.

###### Theorem 1.1 (Minkowski’s existence theorem).

Consider unit vectors spanning and a set of positive weights . Then the equation

 w1n1+w2n2+⋯+wknk=0

is a necessary and sufficient condition for existence of a convex polyhedron with facet unit normals and corresponding facet areas . Furthermore, this polyhedron is unique up to translation.

The necessity part of this theorem is, of course, trivial and follows immediately from the Gauss divergence theorem.

Let us note that in (1), (2), and (3) one deals with the directions , , which can be parameterized by the unit vector . Now consider a closed polygon such that the vectors are the outward normals of the polygon’s edges. Then, according to the Gauss divergence theorem, the edge lengths provide us with the desired set of weights . Vice versa, given a set of unit vectors and positive weights satisfying the first condition in (5), Minkowski’s existense theorem guarantees an existence of a convex polygon with edge lengths and outward normals . The left image of Figure 1 illustrates this geometric interpretation of (5).

In particular, if we consider a polygon circumscribed around the unit circle, as seen in the middle image of Figure 1, the lengths of the polygon sides are given by

 tanβj−1+tanβj,βk=φk+1−φk,

and, therefore, the normalized weights are

 wj=tanβj−1+tanβj∑(tanβk−1+tanβk). (10)

This set of weights was used in Langer et al. (2007) for estimating the mean curvature on surfaces approximated by dense triangle meshes.

In Section 2 we show how to solve this key problem in any dimension and, in particular, get a discretization of the corresponding Laplacian.

### Application to curvature estimation

In the previous examples, we have considered function data defined on a regular grid. Here we deal with less organized data, namely with triangle meshes approximating smooth surfaces embedded into the Euclidean space .

Reliable estimation of curvature characteristics of polygonal surfaces is very important for many computer graphics and geometric modeling applications. Starting from a seminal paper of Taubin Taubin (1995) integral-based approaches are frequently used for curvature estimation purposes Langer et al. (2007); Pottmann et al. (2009); Mérigot et al. (2009) (see also references therein). Our approach to curvature estimation is based on using discrete circular means of directional surface derivatives. Given a smooth surface approximated by a dense triangle mesh, we employ such discrete circular means for estimating the mean curvature and the so-called curvedness Koenderink and Van Doorn (1992) , where and stand for the principal curvatures. (While the curvedness is not a classical curvature measure, it is widely used in numerous applications including protein interaction analysis Bradford and Westhead (2005), heart physiology Zhong et al. (2009), and human visual perception Valenti et al. (2009).)

Let be a surface orientation normal. Denote by and the principal directions corresponding to the principal curvatures and , respectively. For a point on the surface, consider a unit tangent vector , which makes angle with , and the corresponding directional curvature . Similar to (1) and (2) we have

 H≡12(kmax+kmin)=12π∫2π0k(φ)dφ (11)

and

 R2≡12(k2max+k2min)=12π∫2π0(∂n∂t(φ))2dφ, (12)

respectively. Here (11) immediately follows from Euler’s formula

 k(φ)=kmaxcosφ2+kminsinφ2

and (12) is obtained from a similar representation

 (∂n∂t(φ))2=k2maxcos2φ+k2minsin2φ,

which, in turn, can be easily derived from the Rodrigues formulas

 ∂n∂tmax=−kmaxtmax,∂n∂tmin=−kmintmin.

Similar to (6) and (7), let us consider discrete counterparts of (11) and (12)

 H=∑wjk(φj)andR2=∑wj(∂n∂t(φj))2, (13)

where weights satisfy (5) and are normalized by . For the sake of simplicity, we use the set of weights (10) proposed in Langer et al. (2007).

Consider now a triangle mesh approximating the surface and assume that each mesh vertex is equipped with a unit normal which delivers a fairly good approximation of the corresponding surface normal. Given a mesh vertex and its 1-ring neighboring vertices , as seen in the left image of Figure 2, simple approximations of tangential vectors at and the corresponding directional curvature and directional derivative of are given by

 tj≈→PQj∥P−Qj∥,k(φj)≈2→PQj⋅n∥P−Qj∥2,∂n∂tj≈n(Qj)−n(P)∥Qj−P∥, (14)

respectively. We refer to the right image of Figure 2 for a geometric explanation of the directional curvature approximation. Finally, (13) and (14) provide us with discrete approximations of the mean curvatures and curvedness .

Figure 3 presents color encodings for the mean curvature and curvedness detected on a geometrical model with many small-scale details (see the electronic version of this paper for a color version of this figure).

A visual comparison of (13) and (14) with the curvature estimation methods developed in Taubin (1995) and Pottmann et al. (2009) indicates that our discretization is less robust to noise but better reflects small-scale surface details. One can improve the approximation by using more sophisticated approximations for the tangent directions, directional curvatures, and directional derivatives of the normal instead of common-sense estimates (14). Here our goal was to reveal potential applications of our approach to curvature-based shape analysis.

## 2 Discrete spherical means of derivatives

### General construction and Veronese maps

In this section we describe a generalization of the discrete mean value formulas of the preceding section from the circle case to arbitrary higher dimensions. In particular, we prove the following

###### Theorem 2.1.

Given a degree and a collection of points , in general position on the unit sphere in , with , there exist weights (not all equal to zero), which can be found constructively, such that for every harmonic homogeneous polynomial in of degree one has . If , the same statement holds for any homogeneous polynomial in of degree .

For , we have , and this theorem can be viewed as a direct corollary of the Minkowski theorem. Note that for a possible application to computations of the mean curvature of a manifold, the function can be the sum of a constant () and a homogeneous quadratic term (), and one needs to find the weights for obtaining the mean value of such a function from its values at the given set of points. (Of course, all nonzero harmonics should give us the zero mean.) The cases of computations for the discrete Laplacian and square gradient operators are similar.

For instance, we use Theorem 2.1 to construct a discrete Laplacian as follows.

###### Definition 2.2.

Let be a vertex in an arbitrary irregular grid with neighbouring vertices . Define the discrete Laplacian of a function at a point to be

 12Δf(O):=∑wi∂2f∂{e}2Pi, (15)

where points are normalized from to lie on the unit sphere centered at , the weights are given by the above theorem for , and are discrete approximations of the second derivatives of at in the direction .

Definition (15) is inspired by formula (6) and, as we will see later, generalizes the standard discrete Laplacians defined on regular grids.

###### Remark 2.3.

The result of Theorem 2.1 is very close to constructions of so called “spherical designs.” The latter is a finite set of points on such that the average value of any polynomial of degree less than or equal to on this set equals the average value of the polynomial on the sphere. The concept of a spherical design was introduced in Delsarte et al. (1977). The existence and structure of spherical designs in a circle was studied in Hong (1982). In Seymour and Zaslavsky (1984) the existence of designs of all sufficiently large sizes was proved: there is a number depending on and , such that a design exists for any . Many examples and recent results can be found in the survey Bannai and Bannai (2009).

The statement of Theorem 2.1 is somewhat different: we do not require the weights to be equal. In many applications, however, one uses the regular latices, and the weights turn out to be equal, as we discuss below.

### Proof–Construction

First, we prove Theorem 2.1 for a linear functional (i.e. for ). In this case the solution is given by the Minkowski theorem 1.1. Namely, consider the tangent hyperplanes to the sphere at points ; they form a polyhedron. Here we use the fact that , are in general position, so that no hyperplanes are parallel and no more than intersect at a point (note that the polyhedron is not necessarily circumscribed about the sphere, but the hyperplanes containing its facets are tangent to the sphere). Let be the signed codimension-one volumes of the faces of this polyhedron. Then according to the Minkowski theorem (where we understand as vectors in and where the weights may be negative), and hence for every linear function.

Next, consider harmonic polynomials of degree in . Their restrictions on the sphere form the space of spherical harmonics of degree . Consider the harmonic Veronese map , where , that is, the map defined by taking a basis of harmonic polynomials in of degree . The dimension is given by the formula: . Then where is a linear function on . Let be the constructed above weights for the collection of points in . This gives us the required formula . Note that the general position condition on now means that the points in , after normalization, satisfy the above condition that the tangent planes to the sphere form a polyhedron.

Finally consider all homogeneous polynomials in of degree . Their restrictions on the sphere form the space . The dimension of this space is , and we repeat the same construction involving the Veronese map for each homogeneous component as before.

###### Example 2.4.

Let us demonstrate this approach in the case of a circle and quadratic functions. The harmonic Veronese map for quadratic polynomials sends the coordinate functions to the space of quadratic harmonic polynomials (whose restrictions on the sphere give us spherical harmonics). For the circle this harmonic Veronese map is where , or , see a detailed analysis in the preceding section.

Note that if we map to the full 3-dimensional space of quadratic polynomials , then the image of the circle will lie in the plane (due to the relation ), so the point images will not be in general position to yield a polyhedron circumscribed around a 2-sphere in .

###### Example 2.5.

For the case of 2-sphere and quadratic functions the space of harmonic polynomials is 5-dimensional, and our harmonic Veronese map sends to this . The point images lie in general position, and sufficiently many points (at least 6) would yield a polyhedron, which allows one to find the corresponding weights .

###### Remark 2.6.

Since the points should lie on the sphere , rather than be generic vectors in , one can project them to the sphere by rescaling. Namely, let be a (sufficiently generic) collection of non-zero vectors in . Assign the following weights to these vectors. Consider the collection of points on the unit sphere , draw the tangent hyperplanes to the sphere, and let be the signed -dimensional volumes of the faces of the resulting polyhedron.

###### Proposition 2.7.

For every linear function on , one has: , where the weights are .

Proof. The Minkowski formula says that , hence . QED

### Weights for regular lattices

Now, from the general consideration of points in general position we move to regular ones, and consider the set of all “-mid-points,” i.e., the set of points whose coordinates are all possible combinations of nonzero coordinates equal to among spots, starting with to , which are scaled to belong to the unit sphere . Here is any fixed integer between 1 and . We also confine to the case of quadratic harmonics, in the general problem above.

###### Example 2.8.

In dimension 3, for we get consisting of 6 midpoints of the cube faces, which are also the vertices of an octahedron: and . The set consists of 12 scaled mid-points of the cube edges: , and . The set contains exactly 8 vertices of the cube: .

###### Theorem 2.9.

Let be a homogeneous quadratic function, be the set of -mid-points for any fixed , . Then

 1vol(Sn)∫SnF dS=1#Mm∑x∈MmF(x),

where is the volume of -dimensional sphere and is the number of points in .

One can see that in all these “regular cases” one counts the points of with equal weights. This theorem generalizes the corresponding 3D consideration (see e.g. Patra and Karttunen (2005)), which gives three basis stencils used in the Laplacian computations.

Proof. Write , where and a harmonic quadratic polynomial. For a constant the claim is obvious. For one has , so we need to show that for every harmonic quadratic polynomial .

Consider the full Veronese map , given by

 (x1,...,xn+1)↦(x21,...,x2n+1;x1x2,...,xnxn+1),

with . Let be the space of linear functions on that vanish on the vector with first nonzero components (we separate them by semicolon). Then, for each quadratic harmonic , one has for some . Indeed, harmonic polynomials form a hyperplane in the space of linear functions on being given by one linear relation (coming from the sphere equation ), and the polynomials are harmonic and equal to 0 on this vector . Thus we need to check that for all and any .

For instance, for the set , which consists of -type points, the image consists of vectors : with 1’s in one of the first coordinates. In the case of with “diagonal” points

 (0,...,0,±1,0,...,0,±1,0,...,0)/√2,

the image consists of vectors :

 (0,...,0,1,0,...0,1,0,...,0;0,...,0,±1,0,...,0)/2

with all possible pairs of 1’s among first coordinates and the at the respective, -place. Similarly, for any one can see that in each of these cases we have . Hence for all , as claimed. QED.

Thus, while Theorem 2.1 allows one to find the weights for a generic set of vectors , Theorem 2.9 gives the explicit weights in the case when is a homogeneous quadratic polynomial (i.e., ) and vectors correspond to the -mid-points of the -dimensional cube, .

###### Remark 2.10.

For computations we are often interested mostly in the case . Here is the explicit statement for adjusted to this case: let be a homogeneous quadratic function, be the set of 12 points: , , and . Then

 14π∫S2F dS=112∑P∈MF(P).

Thus in all of the above cases we need to use the same weights on these points. This equality of weights for a given set can also be obtained from the symmetry consideration. One can also take combination of the above lattices with with different weights for different (provided that the total sum of the weights is equal to 1), e.g., by considering arbitrary combinations of the 3 lattices in : mid-edges, vertices of the cube, and at the middle of faces. By imposing additional requirements, one can optimize these coefficients: for instance, to achieve an isotropic approximation, or with the next approximation smallest in an appropriate sense.

###### Remark 2.11.

Points of , evidently, form the coordinate -simplex in the space spanned by first coordinates, which is a slim subspace in . In the case of , the points of form the standard -simplex in a hyperplane in . Indeed, in the latter case, we see that for all one has for and , i.e., the vectors have the same length and same angles between them.

###### Remark 2.12.

Above we considered two discretization procedures, both using the Veronese map: a generic set of points and the mid-points sets. While for a generic set in the coefficients are given by the Minkowski formula for the corresponding polyhedron in the image, it is not the case in the mid-point cases. Indeed, due to a very regular structure of the image points , the corresponding tangent planes go to infinity, and the corresponding volumes of the faces are infinite.

In a sense, if the image points are linearly dependent in , and the corresponding polyhedron has faces going to infinity (like in the standard basis case), we have to quotient out these “directions to infinity,” and consider the volume of faces in a polyhedron of smaller dimension.

###### Example 2.13.

For and for the standard basis we did not get a polygon on the plane . Instead, it reduces to a pair of points in . It is not enough to start with two generic points in this case, but if we take three generic points on the circle, we obtain an interpolation formula. This is a manifestation of the general phenomenon described in the following proposition.

###### Proposition 2.14.

There is no discrete mean value formula for quadratic harmonics and a given generic collection of points on .

Proof. Let , and let with be the full Veronese map:

 (x1,...,xn)↦(x21,...,x2n+1;x1x2,...,xixj,...).

Let . Let .

We want to find constants such that . This is equivalent to the following equality:

 Diag(wk)=(A∗A)−1,

where is the matrix made of the vectors , as the following linear algebra shows. Indeed, the equality is equivalent to

 ∑iwixk,ixl,i=δk,l,

which is equivalent to , the identity matrix. This is equivalent to , as claimed.

Of course, if is an orthonormal frame then is orthogonal, and all , as in the theorem above. The matrix is diagonal if and only if the vectors are pairwise orthogonal. This is the case for the set discussed above, which consists of the basis points, but, of course, not true for generic collection of points on . This implies that one cannot have an interpolation formula with generic points. QED

Note that on the other hand, if one has sufficiently many points, so that their Veronese images are linearly dependent in , then any such linear relation provides the weights that yield an interpolation formula.

### Higher-dimensional Laplacian approximation formulas

The results of the present section allow one to easily construct various discrete Laplacians and squared gradient similar to representations (6) and (7), respectively. In the beginning of the section we discussed the construction of the discrete Laplacian for an arbitrary grid, see formula (15):

 12Δf(O)=∑wi∂2f∂{e}2Pi.

To obtain the weights for a point with neigboring vertices we first normalize the latter to get lying on the unit sphere centered at . Then we consider the Veronese map and obtain points . Finally we obtain the weights by finding the signed -dimensional volumes of the faces of the polyhedron in the -space, see Proposition 2.7. An analog of the formula (7) for the square gradient instead of the Laplacian is completely analogous.

For a regular grid, Theorem 2.9 delivers similar formulas for a discrete Laplace operator where all weights are equal:

 Δf(O)=Cm∑∂2f∂{e}2Pi, (16)

where points belong to the grid of -midpoints, , and the constant normalizes the sum, , by comparing the number of neighbours in with that in the coordinate grid , the most straightforward definition of the discrete Laplacian.

Furthermore, one can consider linear combination of the Laplace approximations (16) with different coefficients satisfying . Any such linear combination delivers a different discrete Laplacian and one can run a separate optimization problem on these coefficients to find a discrete Laplacian with better rotation-invariant properties. We discuss the work in this direction for the 2D Laplacian in the next section.

## 3 Discrete circular means of derivatives and isotropic finite differences

In this section we show how to use discrete circular means of Section 1 for approximation of various differential operators in 2D by finite differences with good rotation invariance properties.

The simplest way to construct such approximations consists of substituting in (6), (7), (8), and (9) three-point finite differences instead of the first- and second-order derivatives. For example, constructing a discrete Laplacian corresponding to (6) can be done as follows. Consider a point inside a convex bounded domain . Let the straight line determined by and direction intersect the boundary in two points and . Denote by and the distances from to and , respectively. The second-order -directional derivative of can be approximated by

 ∂2f∂{e}2φ≈2ρ1ρ2([f(q1)ρ1+f(q2)ρ2]/[1ρ1+1ρ2]−f(p)), (17)

where and are assumed to be small. Now a discrete approximation of the Laplacian is obtained by summing up (17) over the set of directions , .

Discrete Laplacian (6), (17) can be considered as a special case of a general construction studied in Wardetzky et al. (2007) (see also references therein).

It is interesting to note that integration of (17) with respect to

 Δf(p)≈2π∫0[(f(q1)ρ1+f(q2)ρ2)/(1ρ1+1ρ2)]dφρ1ρ2/2π∫0dφρ1ρ2−f(p) (18)

leads to an extension of a pseudo-harmonic interpolation scheme proposed in Gordon and Wixom (1974). It can be shown that, given defined on , (18) and its generalizations deliver the harmonic interpolation if is a circle (a ball in the multidimensional case) Belyaev (2006).

Unfortunately combining circular means with directional finite differences does not automatically guarantee rotation invariance properties of the approximating differential operators. The reason is that the finite differences introduce directional errors. However, as seen below, for a regular grid, such directional errors can be uniformly distributed over the directions by choosing appropriate weights in discrete circular means (6), (7), (8), and (9).

### Gradient and Laplacian stencils for square grids

Now let us focus on analyzing rotation invariance properties of regular grid stencils corresponding to (6) and (7).

Consider a two-dimensional square grid with step-size . Each grid vertex has eight nearest neighbors: two horizontal, two vertical, and four diagonal. Thus, according to our geometric interpretation of (6), (7), and (8), we are given four directions and corresponding unit normals

 {e2iφ},φ=0,π/4,π/2,3π/4

which define a rectangle aligned with the coordinate axes.

Consider such a rectangle with the edge lengths and , as shown in the right image of Figure 1, and replace the first- and second-order directional derivatives in (6), (7), and (8) by corresponding central differences. We arrive at the following stencils for the -derivative and Laplacian

 (19)

and

 Δ≈LΔ≡1h2(w+2)⎡⎢⎣1w1w−4(w+1)w1w1⎤⎥⎦, (20)

respectively. The stencil for is obtained from by -rotation. Here and everywhere below the matrices are understood as finite difference operators acting on functions defined on a square grid with step-size .

Formulas (19) and (20) provide us with consistent parameterizations of the stencils for the gradient and Laplacian.

###### Remark 3.1.

The standard central difference and five-point Laplacian is obtained for . Setting leads to the so-called Prewitt masks for the 1st-order derivatives and Laplacian Pratt (2001); Gonzalez and Woods (2008). The case is commonly used in image processing applications and yields the standard Sobel mask for the derivative Pratt (2001); Gonzalez and Woods (2008) and a 9-point discrete Laplacian, which possesses good isotropic properties Kamgar-Parsi et al. (1999). The case was analyzed in Bickley (1948), where it was shown that

 Dx|w=4 = 18h⎡⎢⎣−101−404−101⎤⎥⎦=∂∂x+h212Δ∂∂x+O(h4), (21) LΔ∣∣w=4 = 16h2⎡⎢⎣1414−204141⎤⎥⎦=Δ+h212Δ2+O(h4), (22)

as . Since the Laplacian is isotropic, the right hand sides of (21) and (22) deliver asymptotically optimal asymptotically rotation-equivariant stencils for the -derivative and Laplacian, respectively, as . In particular, this explains why (20) with is often used for a numerical solution of the Poisson equation (Kantorovich and Krylov, 1956, Chapter 3, §1), (Birkhoff and Lynch, 1984, Chapter 3, §10).

Discrete nine-point Laplacian (20) can be represented as a linear combination of two basis five-point Laplacians

 LΔ=αL++βL× with α=ww+2 and β=2w+2, where
 L+=1h2⎡⎢⎣0101−41010⎤⎥⎦andL×=12h2⎡⎢⎣1010−40101⎤⎥⎦.

Setting yields and .

In Figure 4, we use a Gaussian bump function

 g(x,y)=exp(x2+y2) (23)

to test isotropic properties of four discrete stencils

 LΔ∣∣w=∞≡L+,LΔ∣∣w=0≡L×,LΔ∣∣w=2, and LΔ∣∣w=4.

Namely, the figure displays contour lines of . As (22) suggests, with demonstrates the best performance if the grid step-size is sufficiently small.

Similar results are obtained for the -derivative stencils (19) and corresponding -derivative stencil. Analogously to the case of the Laplacian, expansion (21) implies that (19) with is rotationally optimal for sufficiently small .

###### Remark 3.2.

It is easy to increase the approximation accuracy of (21) and (22). For example, one can rewrite (21) as

 ∂∂x=(1+1h2Δ)−1Dx|w=4+O(h4)=(1−1h2Δ)Dx|w=4+O(h4)

and use stencil (22) for approximating .

### Quasi-Laplacian stencil with good rotation invariance properties

Consider now a quasi-Laplacian operator

 La≡∇⋅[a(x,y)∇]=aΔ+∇a⋅∇ (24)

where is given. In view of (22), (21), and (9), the optimal stencils for are obtained from (20) and (19) when .

In practice, the diffusion coefficient in (24) is defined at the staggered grid points which lie on halfway between the standard grid points. Following (Aubert and Kornprobst, 2006, A.3.2) let us approximate (24) using a linear combination of two basis stencils

 αLa++βLa×

with

 La+[f]≡1h2[a+0fi+1,j+a−0fi−1,j+a0+fi,j+1+a0−fi,j−1MMMMMMMMMMMMMMMMMMMMMM−(a+0+a−0+a0++a0−)fi,j]

and

 La×[f]≡12h2[a++fi+1,j+1+a−+fi−1,j+1+a+−fi+1,j−1+a−−fi−1,j−1MMMMMMMMMMMMMMMMMMMMMM−(a+++a−++a+−+a−−)fi,j],

where , , and denote the values of at , , and , respectively.

Assume . Since the same weights are used in the discrete directional mean value representations (6) and (9) we can expect that the optimal rotation-invariant stencil for quasi-Laplacian is obtained when and