# Cubic Bézier radius of curvature calculation

I modelled a racetrack using cubic Bézier curves. I understand that the representation is parametric: x and y vary with t. I have successfully created the cartesian coordinates thanks to this: https://blogs.mathworks.com/graphics/2014/10/13/bezier-curves/

Now I'm trying to calculate the radius of the curvature at each of the x points I have created. I would like to use this expression to calculate the radii (I think this is the easiest way to do this, if there is a simpler way, please feel free to let me know).

Here is how I plan on doing this:

1/ Use the y(t) function for every curve (p24 of this PDF)

2/ Evaluate it at every t

I should have the curvature for every x,y point, right? I'm worried about the 'scaling': t varies between 0 and 1 but x varies between way higher values (for example between 300 and 555). Isn't this going to make the calculation of the radius false?

Thank you very much in advance for your precious help,

Regards, LD

## 1 Answer

That formula is great when $$y$$ can be written as a function of $$x$$, but that doesn't work in cases where the assumption fails...for things like, say, an oval racetrack, where there are two $$y$$-values for every $$x$$-value, or for any more generic curve that might contain a vertical line segment (e.g., a Bezier!).

Instead, you want to use the parametric curvature formula: $$k(t) = \frac{x''(t) y'(t) - y''(t)x'(t)} {(x'(t)^2 + y'(t)^2)^\frac{3}{2}}$$

(This formula assumes that the curve goes counterclockwise; if your curve goes clockwise, then you need to negate it.)

The radius of curvature is then $$1/k(t)$$.

The only problem arises when the denominator in $$k$$ is zero, which happens only when both the $$x$$ and $$y$$ derivatives are zero, i.e., when the curve is not "regular"; in this case, you get a $$0/0$$ form, and the formula's no good to you. In your project, it seems to me that this is unlikely to happen.

Of course, you're inverting $$k$$, so you'll have a problem if $$x''(t) y'(t) - y''(t)x'(t) = 0$$ at some time $$t$$; again, I expect this not to be a problem unless your racetrack has some long "straightaway," where the radius of curvature is infinite.

Here's a bit of sample code to show you what's going on.

function sillycurve()
% Here's an example of the code for you. It actually has the "denominator = 0"
% problem, which I've addressed (crudely) here by adding .001 to the
% denominator, just so you can see the problem nicely.

clear;
clf;
N = 300; % use 300 points
t = linspace(0, 2*pi, N);
dt = 2*pi/N;
% Substitute any formula you like for x and y below.
x = sin(t);
y = cos(2*t) - cos(t);
figure(1);
% Show the curve, growing...
for i=[1:10:N]
plot(x(1:i), y(1:i));
set(gca, 'XLim', [min(x), max(x)], 'YLim', [min(y), max(y)])
figure(gcf);
pause(0.1);
end

figure(gcf);
%
% Compute $$x'(t)$$ and $$x''(t)$$, etc.
xp =  (circshift(x,  1) - x )/dt;
xpp = (circshift(xp, 1) - xp)/dt;
yp =  (circshift(y,  1) - y )/dt;
ypp = (circshift(yp, 1) - yp)/dt;

% Numerator and denominator for curvature formula
kn = xpp .* yp - ypp .* xp;
kd = (xp.^2 + yp.^2).^1.5 + .001; % the .001 is added to avoid divide-by-zero problems

% A plot of the denominator, so that you can see where it goes to zero:
figure(2);
plot(t, kd);
figure(gcf);

% And finally, a plot of the curvature
k = kn ./ kd;
R = 1 ./ k;

figure(3);
plot(t, k);
figure(gcf);

• Thank you very much for your help. Unfortunately, I already know that some curves are pure straights... I'm starting to think that using cartesian coordinates might not be the best of ideas. I saw a few people use curvilinear coordinates and obtain the curvature from those. Would it be easier? I have never used said coordinates... – LTS_student Jun 28 '19 at 10:43
• Well, if you've got straight line segments, then the radius of curvature, along those segments, is infinity. Changing coordinate systems isn't going to change that, as curvature is an implicit notion, independent of coordinates! – John Hughes Jun 28 '19 at 10:46
• If I get infinite values, that's okay. If the code crashes, that's a problem. If x'(t) y''(t) - y'(t)x''(t)=0, then the radius is infinite, right? If that condition is true, then I could skip the rest of the calculation and directly say that the result (the radius) is equal to inf. But will the radius value be correct the rest of the time? This method requires the formula to be evaluated at t, not at the real x point. Will that cause an issue? – LTS_student Jun 28 '19 at 11:03
• That's almost right. If $x'(t) = y'(t) = 0$ as well, then you've got $0/0$, which is an indeterminate form. I very much doubt this will happen. For "the rest of the time," the formula is indeed correct. You're correct that you need the $t$-value in addition to $x$ and $y$. Since you used $t$ to compute $x$ and $y$, you could just store it in an additional field, i.e., you could compute $(x, y, t)$ points, but plot just the $x$- and $y$-coordinates. Then you'd have $t$ around when you need it to compute the curvature. (Or you could compute, at each point, $(x, y, t, k, R)$ and be done!) – John Hughes Jun 28 '19 at 11:11
• Additionally, I could use a polynomial curve fitting function to approximate the cartesian polynomial. But what should be the order of such a polynomial? In other words, if we have a cubic Bézier (order 4), then could the output in cartesian coordinates be also approximated by a 4th order polynomial? – LTS_student Jun 28 '19 at 11:14