#2277 dgrid3d: Mathematical error in the metric computation plot3d.c:qnorm()

[Winfried Ithheld]

• Error: dgrid3d's qnorm interpolation kernels are, in general, not inverse distance weights
• Affected platforms/Gnuplot versions: all since about 1993 [1], here in particular Debian buster (stable)
• Affected files:src/plot3d.c https://sourceforge.net/p/gnuplot/gnuplot-main/ci/master/tree/src/plot3d.c#l603
• Steps to reproduce: Plot scattered data with dgrid3d interpolation with qnorm kernel.
• Expected result: qnorm should be the inverse power distance function for all exponents.
• Actual result: Mathematical meaning of formulas heavily depends on exponent; artificial (lower) symmetry of interpolation kernel in the generic case
• Visual example: Leftmost figure of [2]: fourfold star-shaped halos around points, should be isotropic

Details:

As far as I understand the intentions of dgrid3d, it is meant to implement the Shepard inverse distance weighting for interpolation. But it does so only in some cases for the exponent dgrid3d_norm_value, namely 2, 4, 8, 16.

According to Shepard's original publication [3], and general consensus about Euclidean norm, the inverse power distance weight must read:
qnorm(dist_x,dist_y,q)=(dist_x*dist_x + dist_y*dist_y)**(-0.5*q) (gnuplot function syntax) for all non-zero q. The quadratic exponents for the distance norm (inner parenthesis) are strictly required (polar coordinate conversion). Square roots can only be avoided for even q.

In the qnorm() function in src/plot3d.c, the current code (cases default and 1) implements a function that is mathematically inequivalent to the above universal formula, in its q-dependent exponents in the sum over x and y contributions (these are not simple powers of the Euclidean distance).
One obviously different property is breaking of continuous rotation symmetry. This is easily visible when expressed in polar coordinates: at q=1 it reads 1/(x+y)=1/(r*sin(phi)+r*cos(phi)) and has fourfold (Manhattan-like) symmetry, while the correct 1/sqrt(x*x+y*y)=1/r would be invariant under arbitrary rotations.

The documentation (help dgrid3d) is correct in the sense that it describes the current discontinuous behaviour of the code. But it is imprecise in not stating that (and why) only for some powers of two, the grid interpolation actually represents a textbook inverse distance weighting, and uses some less well studied, anisotropic, fourfold-symmetric kernel in all other cases.

Arguably, the special treatment of the powers of two might also be left to the internal optimisation routines of modern compilers.

• Proposed resolution: Reimplement qnorm() based on the above formula, easily possible by using the existing function pythag() that is found directly below it in plot3d.c.
• Impact: Any dgrid3d-based interpolation for exponents 1, 3, 5-7, 9-15, and above 17 would produce different (correct) results for all end users. Structural behaviour will not change. No side effects expected, except attentive users astonished that the appearance of their data changes.
• Severity/Priority: High from the point of view of mathematical rigor, but limited real-world implications (low implementation risk, little change in output data, probably small userbase)

[1] According to the changelog, Gershon Elber added the code along the "dgrid3d" introduction in v3.5 https://sourceforge.net/p/gnuplot/gnuplot-main/ci/3.5/tree/command.c#l1368
[2] https://en.wikipedia.org/wiki/File:Shepard_interpolation_2.png
[3] https://dl.acm.org/doi/10.1145/800186.810616

[Ethan Merritt]

Do you have a data set for testing?
I coded the change, but the differences are almost too subtle to notice on the demos in the current gnuplot collection. It would help to have a test case.

[Winfried Ithheld]

See attached. In pm3d map view, the Manhattan ridges are clearly visible for the non-binary-exponent cases 15, 17, in case of the square as ridges wrongly radiating out of the central domain, and, more dramatically with the triangle, where domain boundaries are aligned onto a Cartesian grid. The cases with qnorm exponent 16 are isotropic as they should be.

You might also plot the circular plateau with qnorm 1 in normal 3d mode to see how the central tip looks like a pyramid with four sharp edges along the coordinate axes, instead of a cone.

[Ethan Merritt]

Thanks!
Corrected code is now commited to the main branch.

[Winfried Ithheld]

As far as I could check it (withouth comparison to analytical calculations), it now looks qualitatively right (isotropic). Compared to http://www.gnuplot.info/demo/dgrid3d.3.png the gnu-valley demo also lost its spurious groove in the middle.
Difficult to compare with analytical results, but let's face it, interpolation is a dirty business anyway (-:

Thank you!