I recently looked into this topic and here's what I found.
The surfnorm algorithm uses bicubic interpolation as defined in the 1981 Keys paper below, but it relies primarily on calculating the derivatives rather than the full interpolation. It computes finite differences using a 3x3 kernel around each node on the surface, which effectively considers the immediate neighbors. The derivative can be derived from the bicubic interpolant formula after eq 5 of the paper.
The diagonal vectors are essentially the tangential vectors derived from these finite differences. By computing finite differences in both the x and y directions, you get two vectors that lie in the tangent plane of the surface at each point. Crossing these vectors gives you the normal vector.
- Keys, R. (1981). Cubic convolution interpolation for digital image processing. IEEE transactions on acoustics, speech, and signal processing, 29(6), 1153-1160. [link]
It's worth noting that the algorithm in surfnorm has remained quite consistent over time, but as with any software, there's always potential for updates. So, it's a good idea to keep an eye on the documentation and the code itself for the most current information.
