Beach-ball like differential operators of a two-dimensional function

Question

I'm looking for references, known names of, and other useful pointers and insight about (pairs of) differential operators that are "beach-ball like" because they sample a 2-dimensional function in these infinitesimally-sized regular patterns with the indicated alternating polarities:

Figure 1. Pairs of differential operators and a beach ball.

These include a hexagonal and a decagonal pattern. The pairs of operators can be formed by:

$$\begin{gather}\lim_{h\to 0}\frac{\sum_{N=0}^{4N + 1} (-1)^n f\bigg(x + h\cos\left(\frac{2\pi n}{4N + 2}\right), y + h\sin\left(\frac{2\pi n}{4N + 2}\right)\bigg)}{h^{2N + 1}},\\ \lim_{h\to 0}\frac{\sum_{N=0}^{4N + 1} (-1)^n f\bigg(x + h\sin\left(\frac{2\pi n}{4N + 2}\right), y + h\cos\left(\frac{2\pi n}{4N + 2}\right)\bigg)}{h^{2N + 1}},\end{gather}\tag{1}$$

although I'm not certain about the normalization factor $h^{-(2N+1)}$, which at least does not collapse the operator to zero or blow it up to infinity for $N=0$, which is simply a coefficient times differentiation:

$$\begin{gather}N=0:\\ 2\frac{d}{dx}f(x, y),\\ 2\frac{d}{dy}f(x, y),\end{gather}\tag{2}$$

or for $N=1$, which I think is:

$$\begin{gather}N=1:\\ \frac{1}{4}\left(\frac{d}{dx}\right)^3f(x,y)-\frac{3}{4}\frac{d}{dx}\left(\frac{d}{dy}\right)^2f(x, y),\\ \frac{1}{4}\left(\frac{d}{dy}\right)^3f(x,y)-\frac{3}{4}\frac{d}{dy}\left(\frac{d}{dx}\right)^2f(x, y).\end{gather}\tag{3}$$

Applying these to a 2-d Gaussian function and plotting:

Figure 2. Color-mapped 1:1 scale (pixel:unit) plots of, in order: A 2-d Gaussian function with standard deviation $\sigma = 16$, derivative of the Gaussian function with respect to horizontal coordinate $x$, differential operator $\frac{1}{4}\big(\frac{d}{dx}\big)^3-\frac{3}{4}\frac{d}{dx}\big(\frac{d}{dy}\big)^2$ applied to the Gaussian function. Color key: blue: minimum, white: zero, red: maximum.

Python source for Fig. 2:

import matplotlib.pyplot as plt
import numpy as np
import scipy.ndimage

sig = 16  # Standard deviation
N = 161   # Image width
x = np.zeros([N, N])
x[N//2, N//2] = 1
h = scipy.ndimage.gaussian_filter(x, sigma=[sig, sig], order=[0, 0], truncate=(N//2)/sig)
ddx = scipy.ndimage.gaussian_filter(x, sigma=[sig, sig], order=[0, 1], truncate=(N//2)/sig)
h1x = scipy.ndimage.gaussian_filter(x, sigma=[sig, sig], order=[0, 3], truncate=(N//2)/sig) - 3*scipy.ndimage.gaussian_filter(x, sigma=[sig, sig], order=[2, 1], truncate=(N//2)/sig)
plt.imsave('h.png', plt.cm.bwr(plt.Normalize(vmin=-h.max(), vmax=h.max())(h)))
plt.imsave('ddx.png', plt.cm.bwr(plt.Normalize(vmin=-ddx.max(), vmax=ddx.max())(ddx)))
plt.imsave('h1x.png', plt.cm.bwr(plt.Normalize(vmin=-h1x.max(), vmax=h1x.max())(h1x)))
plt.imsave('gaussiankey.png', plt.cm.bwr(np.repeat([(np.arange(N)/(N-1))], 16, 0)))

I found some literature about an application which includes additionally those similar differential operators that would have an even number of samples on each half of the circle in Fig. 1, and (my interpretation aided by @KBDave's answer) represents each pair of operators as the real and the imaginary part of the operator. With those, the rotation between the real and imaginary parts would be such that the imaginary part has samples on the circle (see Fig. 1) half-way between the samples of the real part.

Figure 3. Top: real part, bottom: imaginary part, of a function and complex differential operators applied to it. From Pietro Perona, "Deformable kernels for early vision", Technical Report MIT-LIDS-P-2039, October 1991, also published in April 1995, IEEE Transactions on Pattern Analysis and Machine Intelligence 17(5):222-227.

This gives a more complete set of differential operators to look at.

Answer 1

Working heuristically, suppose that $f(x,y)=\mathrm{e}^{ax+by}$. Then

$$\begin{split}\sum_{n=0}^{4N + 1} (-1)^n f\bigg(x + h\cos\left(\tfrac{2\pi n}{4N + 2}\right), y + h\sin\left(\tfrac{2\pi n}{4N + 2}\right)\bigg)&=\sum_{n=0}^{4N + 1} (-1)^n \mathrm{e}^{ah\cos\left(\tfrac{2\pi n}{4N + 2}\right)+bh\sin\left(\tfrac{2\pi n}{4N + 2}\right)}\\ &=\sum_{n=0}^{4N + 1} (-1)^n\mathrm{e}^{h\Re (c\zeta^{-n})} \end{split}$$ where $c=a+\mathrm{i}b$ and where $\zeta$ is a primitive $(4N+2)$-th root of unity. But $$\sum_{n=0}^{4N + 1} (-1)^n\mathrm{e}^{h\Re (c\zeta^{-n})}=\frac{h^{2N+1}\Re c^{2N+1}}{2^{2N-1}(2N)!} +o(h^{2N+1})\text{,}$$

a result that follows from the Cauchy residue theorem, series expansion in $h$, and

$$\sum_{n=0}^{2N}\frac{1}{z-\cos(\theta+\tfrac{2\pi n}{2N+1})}=\frac{T'_{2N+1}(z)}{T_{2N+1}(z)-\cos(2N+1)\theta}$$ where $T$ denotes a Chebyshev polynomial.

Therefore for sufficiently "nice" $f$ we have

$$\sum_{n=0}^{4N + 1} (-1)^n f\bigg(x + h\cos\left(\tfrac{2\pi n}{4N + 2}\right), y + h\sin\left(\tfrac{2\pi n}{4N + 2}\right)\bigg)=\frac{h^{2N+1}\Re \left((\partial_x+\mathrm{i}\partial_y)^{2N+1}\right)}{2^{2N-1}(2N)!}f +o(h^{2N+1})$$

and it is a matter of analysis to determine the function space in which this argument is precise.