Computing the Radius of Curvature of a Bezier Curve Given Control Points

Question

I'm trying to write a genetic algorithm with curvature as one of the optimization parameters. I want to calculate curvature based on the control points of a bezier curve. I have a minimum radius of curvature that I want to be optimized for. I've been referencing this paper: https://arxiv.org/pdf/1503.01524.pdf

In the paper there is a function that takes the side lengths of a triangle to get the implied radius of curvature, which I have implemented. Here is my current code:

// Computes the curvature implied by 3 control points of a bezier curve
float curvature(float4 p0, float4 p1, float4 p2) {

    // Get the triangle side lengths
    float a = distance(p0, p1);
    float b = distance(p1, p2);
    float c = distance(p2, p0);

    // Do the curvature calculation
    float num = a * b * c;
    float denom = (a + b + c) * (b + c - a) * (a - b + c) * (a + b - c);

    return num / sqrt(denom);

}

The results of this function seem to be incorrect. I run this function for every point in the path, save the last two and then get the minimum radius out of all of them. When I graph the path, there seems to be a major discrepancy between this function's calculation and what I can see visually. Whats the right method to do this?

EDIT: I was looking to calculate the radius of curvature between three control points, not at a given point in the curve, apologies if this was unclear.

Answer 1

The radius of curvature R(t) is equal to 1/κ(t), where κ(t) is the curvature of the curve at point t, which for a parametric planar curve is:

            x'y" - y'x"
κ(t)  = --------------------
         (x'² + y'²)^(3/2)

(Where the ^(3/2) is really ^3/2, but you can't use html for superscript formatting in code blocks, and where the (t) part has been left off the functions for x and y because that makes things unnecessarily hard to read)

As such, for a Quadratic Bezier curve with controls P₁, P₂, and P₃, the first and second derivatives use the following control points:

B(t)': P₁' = 2(P₂ - P₁), and P₂' = 2(P₃ - P₂)
B(t)": P₁" = (P'₂ - P'₁)

Evaluating these for x and y are literally just "using the x or y coordinates", so:

x' = Px₁'(t-1) + Px₂'(t)
y' = Py₁'(t-1) + Py₂'(t)
x" = Px₁"
y" = Py₁"

Noting that x" and y" are just constants. We plug these values into the function for κ(t), provided that the denominator isn't zero, of course (which indicates a line segment, which has no radius of curvature) and then we know what R(t) is because that's just the inverse value.

But we don't actually need to do any of that, because we can make a simple observation: treating the circle as a parametric curve as well, we can find its radius pretty much immediately based on the fact that the vector length of the derivative of the curve, and the derivative of the circle at the same point, are equal. Additionally, we know that the length of the derivative is the same anywhere on a circle (because it's a circle, which by definition is radially symmetric about the origin). As such, we can solve the following:

B(t)  = P1(1-t)^2 + 2*P2(1-t)t + P3t^2
B(t)' = 2*(P2-P1)(1-t) + 2*(P3-P2)t

C(s)  = { r*sin(s), r*cos(s) }independent of 't')
C(s)' = { r*cos(s), -r*sin(s) }

d = |B(t)'| = |C(s = any value, so let's pick 0)'|
d = C(0)' = | (r,0) | = r

And we're done. The radius of curvature at point B(t) is equal to the vector length of B(t)', which is much easier to implement, and much faster to run.

Finally, if what you actually want is to just approximate sections of curve with circular arcs, then https://pomax.github.io/bezierinfo/#arcapproximation should cover the "how to" for that particular use-case.