Curvature of a spline

[Deepak]

I used the Matlab 'spline' function to fit a piecewise cubic spline to my data. Now I want to find the local curvature of this spline at any point. What is the best way to do this?

Thanks,
Deepak

[Bruno Luong]

"Deepak" <deepak....@gmail.com> wrote in message <h9h7be$lh6$1...@fred.mathworks.com>...

> I used the Matlab 'spline' function to fit a piecewise cubic spline to my data. Now I want to find the local curvature of this spline at any point. What is the best way to do this?

K = s'' / (1 + s'^2)^(3/2)

Compute first and second derivatives s' and s" from the pp form.

Bruno

[Didi Cvet]

"Bruno Luong" <b.l...@fogale.findmycountry> wrote in message <h9htqu$qnk$1...@fred.mathworks.com>...

Hi
I opened this thread whit hope that I will find something for my problem but I am not assured. I have a formula for computing spline curvature as integral of ||s''||. And here is how I compute that:
>>Cprim=fnder(C);
>>Csek=fnder(Cprim);
>>ff1=fnval(Csek, param);
>>norma1=ff1(1,:).^2+ff1(2,:).^2;
>>curvature=trapz(norma1);
where C is my spline curve and param are parameters. I think that my algorithm is too sensitive in it's numeric calculations. Does anyone have better idea for calculation of spline curvature. I tried with the above formula (K = s'' / (1 + s'^2)^(3/2)) but I get same numeric instability.
Thanks
Didi

[Antonietta]

"Bruno Luong" <b.l...@fogale.findmycountry> wrote in message <h9htqu$qnk$1...@fred.mathworks.com>...

Dear Bruno,
Could you give further details on curvature of a spline calculation? I tried applying the equation K = s'' / (1 + s'^2)^(3/2) to my data but I got wired results.

F=spline(1:23,B); % B is a 3x23 matrix

% first derivative calculation
ppFirstDer = F;
A1=zeros(size(F.coefs));
A1(:,2)=3*F.coefs(:,1);
A1(:,3)=2*F.coefs(:,2);
A1(:,4)=1*F.coefs(:,3);
ppFirstDer.coefs=A1;

% second derivative calculation
ppSecondDer = F;
A=zeros(size(F.coefs));
A2(:,3)=6*F.coefs(:,2);
A2(:,4)=2*F.coefs(:,3);
ppSecondDer.coefs=A2;

kappa=F;
kappa.coefs= A2./(1 + A1.^2).^(3/2);

Did I made something wrong?

[Antonietta]

"Bruno Luong" <b.l...@fogale.findmycountry> wrote in message <h9htqu$qnk$1...@fred.mathworks.com>...

[Bruno Luong]

Here we go:

% Input:!: Generate the 10 points on the circle of radius 20
x = 1:10;
R = 20;
y = -sqrt(R^2-x.^2);

% The spline that interpolates (x,y)
s = spline(x,y);

% first derivative, pp form
s1 = s;
s1.order = s.order-1;
s1.coefs = bsxfun(@times, s.coefs(:,1:end-1), s1.order:-1:1);

% second derivative, pp form
s2 = s1;
s2.order = s1.order-1;
s2.coefs = bsxfun(@times, s1.coefs(:,1:end-1), s2.order:-1:1);

% Evaluate the curvature
% 100 points cover the same interval
x = linspace(x(1),x(end));
A1 = ppval(s1,x); % first derivative at x
A2 = ppval(s2,x); % second derivative at x
% curvature
kappa = A2./(1 + A1.^2).^(3/2);

% radius of curvature, should be close to 20
r = 1./kappa;

% Graphic check
plot(x, r);

% Bruno