Problem in curvature calculation using cubic spline interpolation of arc-length parameterized 2D curve

Hello all.

I have a 2D parametric curve C

, , t is a parameter which lies in [0,15].

I have numerically calculated the arc length (s) of the above curve by using discrete values of 't', say, t=0:(15/500):15, ie., using 501 points.

Thus, I got a set of points x(s) and y(s) and a set of arc-length (S),all containing 501 elements each.

I used cubic spline to fit the data set x(s) and y(s) using S.

Finally I have a x(s) and y(s) independantly piece-wise cubic polynomials in arc length.

The derivatives , , , and are now easily computed by differentiating the appropriate segments of the piece-wise cubic . The slopes and the curvatures of the parameterized curce (C) are then found out using the following equations.

I have compared the actual and fitted values for the curve. The plots are agreeing for function values and the first derivatives. But for curvature, the fitted values of the curvature shows unexpected oscillations at some points. My prime requirement is a smooth curvature function in arc length parameter. Please find the attached plots for comparison of values.

Kindly help me to rectify the issue in the curvature values. Note sure why this kind of variation happens in the curvature, eventhough the function and the slope values are smooth.

%% 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;
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_s*y_fitted_s_s-y_fitted_s*x_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