How to find curvature(k) of plane curve have having using it's position points (x & y) equally spaced with arc length s.

Dear Seniors,

I'm researching on a topic and from many days, i'm so confused, i'm learning to get strong programmming skills, but rightnow i don't know how to get out of this. I used ''diff'' and ''syms'' to solve, but it is not effective.

I have a cubic spline curve having arrays of psition x and y equally parsmertize with arclength s. I want to get derivative of x points with respect to s, derivative of y point w.r.t s. and then put in formula to find curvature (k), but i can't get dy/ds & dx/ds.

points(x&y)=[328.1228 201.8773; 326.2455 203.7545; 324.3683 205.6317; 322.4910 207.5090; 320.6138 209.3826; 318.7365 211.22634; 316.8593 213.1407; 314.9820 215.0180; 313.1048 216.8952; 311.2276 218.7725; 309.3504 220.6489; 307.4731 222.5269; 305.5958 224.4041; 303.7189 226.2818; 301.8424 228.1598]; each position is equally spaced difference of s=2.6548.

for continues acrlength w.r.t position(x&y) along curve i used

acrlen=sqrt((diff(points(:,1))).^2+(diff(points(:,2))).^2);
s=[0 cumsum(arclen')];

[YOUNG MIN CHOI]

Hi, Xinyue Lan

I have a same problem with you.

"cubic spline curve arc-length paramertization"

If you have solved this problem, Can you provide MATLAB code?

[e-mail]

Answer 1

[Bruno Luong]

Big Curvature -> red

Small Curvature -> blue

% Dummy test data
t = 1:7;
x = rand(size(t));
y = rand(size(t));

fx = spline(t,x);
fy = spline(t,y);
d1fx = differentiate(fx);
d1fy = differentiate(fy);
d2fx = differentiate(d1fx);
d2fy = differentiate(d1fy);

ti = linspace(min(t),max(t),1000);
xi = ppval(fx,ti);
yi = ppval(fy,ti);
d1xi = ppval(d1fx,ti);
d1yi = ppval(d1fy,ti);
d2xi = ppval(d2fx,ti);
d2yi = ppval(d2fy,ti);
% EDIT 30/May/2024 BUG fix, signale by @Robert Owens; thanks
K = abs(d1xi.*d2yi - d2xi.*d1yi) ./ sqrt(d1xi.^2+d1yi.^2).^3; 

figure(1)
clf
plot(x,y,'ow');
hold on
% FEX https://fr.mathworks.com/matlabcentral/fileexchange/60402-plotcolor
plot_color(xi,yi,log(K),jet,[],'Linewidth',2);
set(gca,'Color','k');
axis equal

function ppdf = differentiate(ppf)
% Spline Derivative
ppdf = ppf;
ppdf.order=ppdf.order-1;
ppdf.coefs=ppdf.coefs(:,1:end-1).*(ppdf.order:-1:1);
end

[Francesco Pignatelli]

Hi @Bruno Luong, thanks for this code, it is helping me a lot. However I have a small doubt: if I have set of data (x,y) how is it possible to perform a cubic spline curve arc-length paramertization before differentianting?

[Bruno Luong]

@Francesco Pignatelli I'm not sure if your issue is related to the origonal question. It seem you should take a look at John submission interparc and see it it fits your need.

[Robert Owens]

Hi @Bruno Luong, please note a minor error in the formula for the curvature.

Answer 2

[John D'Errico]

Why does it seem you are trying to do this the hard way? Seriously, you are working WAY too hard on this.

plot(x,y,'o-'),grid on

Your data is so close to being a straight line, that I have a hard time seeing any deviation from a line.

So why in the name of god and little green apples did you want to make a parametric function out of this? Seriously.

spl = spline(x,y);

Now, I'll assume that you wish to plot the curvature of this relation, as opposed to the radius of curvature.

https://en.wikipedia.org/wiki/Curvature

There we see that the signed curvature is:

K = f' '(x)/(1 + (f'(x))^2)^(3/2)

If you don't care about the sign, then take the absolute value. Regardless, it is easy enough to write now.

fp = fnder(spl);
fpp = fnder(spl,2);
K = @(x) fnval(fpp,x)./(1 + fnval(fp,x).^2).^(3/2);
ezplot(K,[min(x),max(x)])
axis([min(x),max(x),-1e-2,1.5e-2])
grid on

You can't expect too much more than a piecewise linear function from a simple simple here. I suppose I could have done the original fit using a higher order spline. This would work for example:

spl = spapi(5,x,y)

fnder is part of the curve fitting toolbox, but is not that difficult to implement otherwise.

[Xinyue Lan]

Dear Sir

thank you so much for your reply, but i shared just a little part of my curve. I'm doing path planning functionality of mobile robot,theresofe i need to find predifined waypoints data, to guide robot(including curvature and acrlength).

This is my curve from x & y points obtain from cubic spline curve arclength paramertization .

[Xinyue Lan]

I'm talking about this formula , here t=s

[YOUNG MIN CHOI]

Hi, Xinyue Lan

I have a same problem with you.

"cubic spline curve arc-length paramertization"

If you have solved this problem, Can you provide MATLAB code?

[e-mail]

[Francesco Pignatelli]

Hi all,

I have the same propled as yours. Did you fix it? Can I have a look into the code?