How to find curvature(k) of plane curve have having using it's position points (x & y) equally spaced with arc length s.
40 views (last 30 days)
Show older comments
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')];
1 Comment
YOUNG MIN CHOI
on 30 Jul 2020
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?
coym94@naver.com
Accepted Answer
Bruno Luong
on 3 Jan 2019
Edited: Bruno Luong
on 3 Jan 2019
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);
K = abs(d1xi.*d2yi - d2xi.*d2yi) ./ 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
2 Comments
Francesco Pignatelli
on 23 Dec 2022
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
on 24 Dec 2022
Edited: Bruno Luong
on 24 Dec 2022
@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.
More Answers (1)
John D'Errico
on 3 Jan 2019
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.
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.
4 Comments
YOUNG MIN CHOI
on 30 Jul 2020
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?
coym94@naver.com
Francesco Pignatelli
on 23 Dec 2022
Hi all,
I have the same propled as yours. Did you fix it? Can I have a look into the code?
See Also
Categories
Find more on Spline Postprocessing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!