How to calculate curvature of isoline of image

I've got a discrete grayscale image (a.k.a. 2D rectangular array of floats). Its continuous representation (whether received via reconstruction with sinc, with a cubic kernel, with a triangle kernel...) is a continous 2D function, so it has isolines. I'm interested in calculating the magnitude of the derivative of the isoline at a given point.

Pic: Red line is part of an isoline, and green arrows are example derivatives of that isoline at 2 points.

Here's my approximate solution:

float getCurvature(vec2 p) {
if(grad = vec2(0, 0)) return 0;

float val = getBilinear(img, p);
return (val - (valLeft + valRight) * .5f);
}

It basically computes the second directional derivative in a direction perpendicular to the gradient. It works fine for my purposes, except that it's somewhat slow in the CPU implementation (and I have reasons not to port it to GPU).

Any faster (and/or more accurate) method?

Edit: I found a twice-faster method, still calculating the second derivative, but in a simpler way. It's equation (5) on this page. It works, but the result is more crude. I think I should switch to a higher-quality differentiation kernel than [+1, 0, -1].