How to find curvature(k) of plane curve have having using it's position points (x & y) equally spaced with arc length s.
조회 수: 39 (최근 30일)
이전 댓글 표시
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
YOUNG MIN CHOI
2020년 7월 30일
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
채택된 답변
Bruno Luong
2019년 1월 3일
편집: Bruno Luong
2019년 1월 3일
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
Francesco Pignatelli
2022년 12월 23일
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
2022년 12월 24일
편집: Bruno Luong
2022년 12월 24일
@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.
추가 답변 (1개)
John D'Errico
2019년 1월 3일
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
YOUNG MIN CHOI
2020년 7월 30일
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
2022년 12월 23일
Hi all,
I have the same propled as yours. Did you fix it? Can I have a look into the code?
참고 항목
카테고리
Help Center 및 File Exchange에서 Spline Postprocessing에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!