How to calculate curvature of isoline of image

Question

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) {
    vec2 grad = getBilinear(gradients, p);
    if(grad = vec2(0, 0)) return 0;
    else grad = normalize(grad);

    vec2 gradP = perpendicularLeft(grad);

    float val = getBilinear(img, p);
    float valLeft = getBilinear(img, p+gradP);
    float valRight = getBilinear(img, p-gradP);
    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].

Answer 1

Gradient magnitude and angle can be found using a Sobel operator. It's more sophisticated than the [1, 0, -1]. This is commonly used either directly for edge detection, or as part of a more complex algorithm.

Or, if I've misunderstood, and you want only the curvature of the 2D isoline, this answer on Math Overflow might be of use.