I have a 2D curve in the $xy$-plane, which was arc length parameterized numerically, and fitted by cubic splines for both $x$ and $y$. If one of the segments of the cubic spline is:
\begin{align} x&=a_1s^3 + a_2s^2 + a_3s + a_4 \\ y&=b_1s^3 + b_2s^2 + b_3s + b_4, \end{align}
where $a_1$, $a_2$, $a_3$, $a_4$, $b_1$, $b_2$, $b_3$, and $b_4$ are constants and $s$ is the arc length parameter.
How can I find the curvature $k(s)$ of the curve using the above equations?
Update in the question:
Please find the attached image for the details of my problem.
A MATLAB code for the problem is also given below:
%% Original curve
Lt=15; %length of the parameter 't'
N=500; %Number of points on the curve
h=Lt/N; %Step size along t
t = 0:h:Lt;
% Function definition
x= log(2+t);
y= log(1+t);
%% Actual derivative, and curvature
slope=gradient(y)./gradient(x); % derivative
second_derivative=gradient(slope)./gradient(x); % second derivative
curv=second_derivative./(1+(slope).^2).^(3/2); % curvature at any point
%% Arc length parameterization
x_t = gradient(x);
y_t = gradient(y);
s = cumtrapz( sqrt(x_t.^2 + y_t.^2 ) ); % Arc length
X = s.';
V = [ t.', x.', y.' ];
L = s(end); % Total length of the arc
s0 = s(1);
Ni = length(s);
Xq = linspace(s0,L,Ni); % equally spaced indices
Vq = interp1(X,V,Xq);
xs = Vq(:,2); % arc length parameterized x
ys = Vq(:,3); %arc length parameterized y
%% Cubic spline interpolation
pp1 = csaps(X, xs);
Cx=pp1.coefs; % getting the coefficients of the piece-wise polynomials
pp2 = csaps(X, ys);
Cy=pp2.coefs;
%% Comparing the actual function with the fitted ones
fitted_x=[];fitted_y=[];fitted_slope=[];fitted_curvature=[];
for i=1:N
out=fitted_values(Cx,Cy,X,i);
fitted_x=cat(1,fitted_x,out(1));
fitted_y=cat(1,fitted_y,out(2));
fitted_slope=cat(1,fitted_slope,out(3));
fitted_curvature=cat(1,fitted_curvature,out(4));
end
figure;
plot(x,y,'b'); % Actual function
hold on;
plot(fitted_x,fitted_y,'r-'); % Fitted function
legend('Actual function','Fitted function');
figure;
plot(x,slope,'b'); % Actual slope
hold on;
plot(fitted_x,fitted_slope,'r-'); % Fitted slope
legend('Actual slope','Fitted slope');
figure;
plot(x,curv,'b'); % Actual curvature
hold on;
plot(fitted_x,fitted_curvature,'r-'); % Fitted slope
legend('Actual curvature','Fitted curvature');
function fval=fitted_values(Cx,Cy,X,i)
s=X(i); % arc length parameter
if i==1
s1=s;
else
s1=X(i-1);
end
spline_segment_number=i-1; % At which cubic spline segment the evaluation is to be performed
if i==1
spline_segment_number=1;
end
coeff_x=Cx(spline_segment_number,:);
coeff_y=Cy(spline_segment_number,:);
% Fitted function values
x_fitted=coeff_x(1)(s-s1)^3+coeff_x(2)(s-s1)^2+coeff_x(3)(s-s1)+coeff_x(4);
y_fitted=coeff_y(1)(s-s1)^3+coeff_y(2)(s-s1)^2+coeff_y(3)(s-s1)+coeff_y(4);
% Fitted slope values
x_fitted_s=3*coeff_x(1)*(s-s1)^2+2*coeff_x(2)*(s-s1)^1+coeff_x(3);
y_fitted_s=3*coeff_y(1)*(s-s1)^2+2*coeff_y(2)*(s-s1)^1+coeff_y(3);
slope_fitted=y_fitted_s/x_fitted_s;
% Fitted curvature values
x_fitted_s_s=6*coeff_x(1)*(s-s1)^1+2*coeff_x(2);
y_fitted_s_s=6*coeff_y(1)*(s-s1)^1+2*coeff_y(2);
curvature_fitted=(x_fitted_sy_fitted_s_s-y_fitted_sx_fitted_s_s)/(x_fitted_s^2+y_fitted_s^2)^(3/2); % using the curvature formula
fval=[x_fitted;y_fitted;slope_fitted;curvature_fitted];
end
First we differentiate $(x,y)$ with respect to arc $s$ for input to curvature formula.
$$x(s)= a_1 s^3 + a_2 s^2 + a_3 s + a_4 ;\; $$ $$ x'(s)= 3 a_1 s^2 + 2 a_2 s + a_3;\; x''(s)=6 a_1 s + 2 a_2 ;\; $$
similarly find $ y'(s),y''(s) $. Compute curvature with formula:
$$ k_g(s) =\dfrac{x'y''-y'x''}{(x^{'2}+ y^{'2})^\frac32} = \dfrac{d\phi}{ds}$$
$$ \phi= \int k_g\;ds\;$$
where $\phi $ is slope obtained by integrating back numerically. Next $(x,y)$ are integrated for, using any CAS:
$$x(s)= \int \cos \phi\; ds, \quad y(s)= \int \sin \phi \; ds. $$