8
$\begingroup$

I've got a 2D cubic spline (Bézier) and I have the polygon-line that's a discretization of that spline.

Is there an efficient and simple to implement a way to calculate the maximum curvature of the spline either by using its smooth representation or the discrete poly-line? It doesn't have to be exact since it's for a game, an error around 10% would be totally acceptable.

Different wording: I would like to know: if I'd drive a car driving at constant speed along the spline, what would be the maximum degree I had to turn the steering wheel to stay on it?

Improve this question
$\endgroup$
16
  • 1
    $\begingroup$ Taking derivative of the spline that you know did not helped? Two statements are different in meaning. Do you have spline type? Like Bézier? $\endgroup$
    – Evil
    Commented Jan 13, 2016 at 20:48
  • $\begingroup$ Sure, It's a cubic Bezier spline. AFAIK there is no analytical way to determin the max curvature of a cubic spline. $\endgroup$
    – user3049681
    Commented Jan 13, 2016 at 21:00
  • $\begingroup$ Cubic spline means it has cube coefficient. Type Bézier, ok. This is 2D case? 3D case? (For car its rather 2D). First step is to clear what you have, and what you really need. $\endgroup$
    – Evil
    Commented Jan 13, 2016 at 21:02
  • $\begingroup$ Thanks, it's the 2D case, I've got the cubic bezier spline and I've got a discretized Version in form of a polyline. From those two Things I'd like to somehow roughly determin the Maximum curvature. $\endgroup$
    – user3049681
    Commented Jan 13, 2016 at 21:08
  • 1
    $\begingroup$ 1. How do you want to define "curvature"? 2. Note that your second question (about the car) is not well-defined, since you haven't told us anything about how the speed of the car varies over time. Did you intend to specify that the car travels at constant velocity? If so, that question is equivalent to asking at what point maximizes the magnitude of the second derivative $[x''(t)]^2 + [y''(t)]^2$, assuming you choose a parameterization so that the velocity is constant ($[x'(t)]^2 + [y'(t)]^2 = 1$). $\endgroup$
    – D.W. ♦
    Commented Jan 13, 2016 at 21:38

2 Answers 2

Reset to default
2
$\begingroup$

Step 1: Express the points on the spline parametrically, so the spline is the set of points of the form $(x(t),y(t))$, where $t$ is a parameter. Here $x(t)$ represents the $x$-coordinate (as a function of the parameter $t$) and $y(t)$ represents the $y$-coordinate. Since this is a cubic spline, you can find functions $x(t),y(t)$ that are cubic polynomials with known coefficients that provide this parametric expression.

Step 2: Use the formula on Wikipedia for computing the curvature, given a parametric representation of the curve. This gives you a formula for the curvature as a function of $t$, namely, $\kappa(t)$. Note that since $x(t)$ and $y(t)$ are cubic polynomials, you can explicitly compute their first and second derivatives, so you can analytically compute an explicit expression for $\kappa(t)$, i.e., for the curvature as a function of $t$.

Step 3: Find the value of $t$ that maximizes $\kappa(t)$. Note that we are now dealing with a function $\kappa : \mathbb{R} \to \mathbb{R}$, i.e., we are in the one-dimensional case. Thus, we can find the maximum numerically using any of a number of methods: gradient descent, Newton's method, or a number of other methods.

Alternatively, you can analytically compute the derivative of $\kappa(t)$, and then solve the equation $\kappa'(t)=0$ for $t$. This might allow an analytic solution that identifies a list of candidate maxima of $\kappa(t)$. Make sure to also check the endpoints (the lower and upper bounds for $t$). Evaluate $\kappa$ at each candidate and pick the one that makes the value of $\kappa(t)$ as large as possible.

Improve this answer
$\endgroup$
2
$\begingroup$

I followed D.W.'s approach and found a solution as follows:

In order to use an online derivative calculator, I expressed the cubic spline with points P0, P1, P2 and P3 using variables a through h:

P0 = (a, e)
P1 = (b, f)
P2 = (c, g)
P3 = (d, h)

The x-y-coordinates and their derivatives are as follows:

x(t)   = a*(1-t)^3+3*b*(1-t)^2*t+3*c*(1-t)*t^2+d*t^3
x'(t)  = 3*((d-3*c+3*b-a)*t^2+2*(c-2*b+a)*t+(b-a))
x''(t) = 6*((d-3*c+3*b-a)*t+(c-2*b+a))

y(t)   = e*(1-t)^3+3*f*(1-t)^2*t+3*g*(1-t)*t^2+h*t^3
y'(t)  = 3*((h-3*g+3*f-e)*t^2+2*(g-2*f+e)*t+(f-e))
y''(t) = 6*((h-3*g+3*f-e)*t+(g-2*f+e))

In order to keep the equations in a more compact form, I introduced the following substitutions:

m = d-3*c+3*b-a
n = c-2*b+a
o = b-a
p = h-3*g+3*f-e
q = g-2*f+e
r = f-e

This leads to the following simplified derivatives:

x'(t)  = 3*(m*t^2+2*n*t+o)
x''(t) = 6*(m*t+n)
y'(t)  = 3*(p*t^2+2*q*t+r)
y''(t) = 6*(p*t+q)

Using the mentioned curvature equation and the online derivative calculator I obtained the following expression with its derivative:

k(t)   = (x' * y'' - y' * x'') / (x'^2 + y'^2)^(3/2)
       = ((3*(m*t^2+2*n*t+o)) * (6*(p*t+q)) - (3*(p*t^2+2*q*t+r)) * (6*(m*t+n))) / ((3*(m*t^2+2*n*t+o))^2 + (3*(p*t^2+2*q*t+r))^2)^(3/2)
k'(t)  = (-18*m*(p*t^2+2*q*t+r)+18*p*(m*t^2+2*n*t+o)-18*(m*t+n)*(2*p*t+2*q)+18*(2*m*t+2*n)*(p*t+q))/(9*(p*t^2+2*q*t+r)^2+9*(m*t^2+2*n*t+o)^2)^(3/2)-(3*(18*(p*t+q)*(m*t^2+2*n*t+o)-18*(m*t+n)*(p*t^2+2*q*t+r))*(18*(2*p*t+2*q)*(p*t^2+2*q*t+r)+18*(2*m*t+2*n)*(m*t^2+2*n*t+o)))/(2*(9*(p*t^2+2*q*t+r)^2+9*(m*t^2+2*n*t+o)^2)^(5/2))

Unfortunately, the online tool is not able to yield roots for this expression. But if you look closely (preferably at some rendered representation), you'll find a difference with almost identical denominators (...)^(3/2) and 2*(...)^(5/2). Multiplying the left-hand side with 2*(...) yields a common denominator and the following nominator, which should be set to zero:

2*(-18*m*(p*t^2+2*q*t+r)+18*p*(m*t^2+2*n*t+o)-18*(m*t+n)*(2*p*t+2*q)+18*(2*m*t+2*n)*(p*t+q))*(9*(p*t^2+2*q*t+r)^2+9*(m*t^2+2*n*t+o)^2) - (3*(18*(p*t+q)*(m*t^2+2*n*t+o)-18*(m*t+n)*(p*t^2+2*q*t+r))*(18*(2*p*t+2*q)*(p*t^2+2*q*t+r)+18*(2*m*t+2*n)*(m*t^2+2*n*t+o)))

This tool can simplify this expression and yields the following polynomial in text form:

  ((1296*m*p^2+1296*m^3)*q-1296*n*p^3-1296*m^2*n*p)*t^5
+ ((1620*m*p^2+1620*m^3)*r+3240*m*p*q^2+(3240*m^2*n-3240*n*p^2)*q-1620*o*p^3+((-1620*m^2*o)-3240*m*n^2)*p)*t^4
+ ((5184*m*p*q+1296*n*p^2+6480*m^2*n)*r+1296*m*q^3-1296*n*p*q^2+((-6480*o*p^2)-1296*m^2*o+1296*m*n^2)*q+((-5184*m*n*o)-1296*n^3)*p)*t^3
+ (1296*m*p*r^2+(1944*m*q^2+6480*n*p*q-1296*o*p^2+1296*m^2*o+8424*m*n^2)*r-8424*o*p*q^2-6480*m*n*o*q+((-1296*m*o^2)-1944*n^2*o)*p)*t^2
+ (2592*n*p*r^2+(3888*n*q^2-2592*o*p*q+2592*m*n*o+3888*n^3)*r-3888*o*q^3+((-2592*m*o^2)-3888*n^2*o)*q)*t
+ (-324*m*r^3+(1944*n*q+324*o*p)*r^2+((-1944*o*q^2)-324*m*o^2+1944*n^2*o)*r-1944*n*o^2*q+324*o^3*p)

Now you can choose any solver for 5th-degree polynomials. In my case, I could easily use the roots function from Python's NumPy:

poly = [
    (1296*m*p**2+1296*m**3)*q-1296*n*p**3-1296*m**2*n*p,
    (1620*m*p**2+1620*m**3)*r+3240*m*p*q**2+(3240*m**2*n-3240*n*p**2)*q-1620*o*p**3+((-1620*m**2*o)-3240*m*n**2)*p,
    (5184*m*p*q+1296*n*p**2+6480*m**2*n)*r+1296*m*q**3-1296*n*p*q**2+((-6480*o*p**2)-1296*m**2*o+1296*m*n**2)*q+((-5184*m*n*o)-1296*n**3)*p,
    1296*m*p*r**2+(1944*m*q**2+6480*n*p*q-1296*o*p**2+1296*m**2*o+8424*m*n**2)*r-8424*o*p*q**2-6480*m*n*o*q+((-1296*m*o**2)-1944*n**2*o)*p,
    2592*n*p*r**2+(3888*n*q**2-2592*o*p*q+2592*m*n*o+3888*n**3)*r-3888*o*q**3+((-2592*m*o**2)-3888*n**2*o)*q,
    -324*m*r**3+(1944*n*q+324*o*p)*r**2+((-1944*o*q**2)-324*m*o**2+1944*n**2*o)*r-1944*n*o**2*q+324*o**3*p,
]
roots = np.roots(poly)

The resulting list contains all values t, where the curvature k(t) is at a local minimum or maximum. There could, however, be imaginary solutions that should be ignored.

Example:

example

Regarding D.W.'s hint about endpoints: I'm not sure if the curvature could be extrem at these points. But if in doubt, make sure to check the endpoints explicitly.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.