Numerical approximation of principal curvature

I have a surface given by z-values on an xy-grid (a 2D-array of values).
To calculate surface tension, I need to calculate mean curvature in every point.
To calculate mean curvature, I need to calculate principal curvature.
Since principal curvature lies in a plane containing the surface-normal-vector, this is where it gets tricky.

My current approach is very brute force:

1. take a point and it's surrounding 8 points
2. calculate gradient at central point
3. rotate all points so that gradient is now in the xy-plane (moving all points off the grid!)
4. calculate intersection of triangles defined by these 8 points and a sample of planes normal to the xy-plane
5. find maximum

This is of course very cumbersome and slow. I'm furthermore not sure how many normal plains I should sample and weather there are special planes that give better estimates (e.g. planes containing the rotated points). Is there a more elegant way to do this?

This may be a duplicate: Numerical computation of surface curvature

• I found somewhat of a solution to this problem, though I'm not sure it's applicable: math.rug.nl/~veldman/Scripties/Lam-MasterTechWisk.pdf (chapter 4) This thesis uses value-fractions to calculate curvature on a purely 2D grid. Since I require this algorithm only for pretty graphics, I'll use it. I believe it isn't an actual solution to the problem, though! (Since the curvature is taken along x-y orientation instead of applying principal curvature calculations) Aug 6 '18 at 7:10

Curvature defined as $$\nabla \vec{e}_n$$ can be calculated with a suitably defined $$\vec n$$. Having a grid of Z-values, a vector normal to the surface can be found by using the gradient: $$\vec n = \binom{1}{\vec \nabla Z }$$ With $$\vec \nabla z$$ supplying $$x$$ and $$y$$ values of $$\vec n$$ and $$z=1$$. Thereby: $$\vec e _n = {\vec{n}\over{|\vec n|}} = {\vec{n}\over \sqrt{1+(\vec \nabla Z)^2}}$$