Calculating arc length of piece wise cubic polynomial

%
%
% Spline Generator
%   Takes an excel file with a list of inflection points for a specific 
%   course. The inflection points model a spline that follows the racing 
%   line from start to finish for autocross or a complete lap. The racecar
%   will move along the spline, representing its path along the course. 
%

clc;
clearvars;

%Imports (x,y) coordinants for inflection points from excel file
InflectionPoints = readmatrix('Auto2017InflectionPoints.xlsx');

%Time vector, interval, and inflection points
tpts = [0 600];
tvec = 0:0.1:600;
cp = transpose(InflectionPoints);

%Creates a peicewise polynomial from spline inflection points
[Position, qd, qdd, path] = bsplinepolytraj(cp,tpts,tvec);
 
 j = 1;
 length = 0;
 for i = 2:2:118
    %Flag for checking arc length
     checkedal = false;
    
    %Setup parametric eqs for current section 
    a = path.coefs(i-1,1);
    b = path.coefs(i-1,2);
    c = path.coefs(i-1,3);
    d = path.coefs(i-1,4);
    
    e = path.coefs(i,1);
    f = path.coefs(i,2);
    g = path.coefs(i,3);
    h = path.coefs(i,4);
    
    ent = path.breaks(1,j);
    ext = path.breaks(1,j+1);
    
    t = ent:0.01:ext;
    
    x = a*(t-ent).^3 + b*(t-ent).^2 + c*(t-ent) + d;
    y = e*(t-ent).^3 + f*(t-ent).^2 + g*(t-ent) + h;
    
    %Overlay section onto plot of sections
    hold on;
    plot(x,y);
   
   %Calculate Curvature and Arc Length of the section
   u=1;
   for t = ent:0.1:ext
   
    xd = 3*a*(t-ent)^2 + 2*b*(t-ent)+ c;
    xdd = 6*a*(t-ent) + 2*b;
    
    yd = 3*e*(t-ent)^2 + 2*f*(t-ent)+ g;
    ydd = 6*e*(t-ent) + 2*f;
    
    K = (abs(xd*ydd-yd*xdd))/((xd^2+yd^2)^(3/2));
    R = 1/K;
    
    radiusdata(u,2*j-1) = R;
    
    if checkedal == false
    length = length + arclength(a,b,c,e,f,g,ent,ext);
    radiusdata(u,2*j) = length;
    checkedal = true;
    end
    u=u+1;
   end 
    j=j+1;
 end
 

 function [arcl] = arclength(ax,bx,cx,ay,by,cy,a,b)
%arclength Calculates the arclength for a parametric cubic polynomial
%   INPUT:
%       Coefficients of parametric equations()
%       

f = @(t) sqrt((3*ax*t.^2 + 2*bx*t+ cx).^2 + (3*ay*t.^2 + 2*by*t+ cy).^2);
arcl = integral(f,a,b);

end

I am trying to calculate the arc length of the peices to ultimately get a matrix of, the distance along the curve and the radius at that distance, at specified intervals. I know the entire curve should have an arc length of roughly 32,000, however I am off by orders of magnitude. The arc length of each succesive peice is no longer than the previous one, but the peices should be within an order of magnitude in length.

I have tried looking at the arc length of the first t interval of each section and generally the arc length would become bigger the further down the curve the section is.

Can anyone spot what I'm doing wrong? I've been at this for hours and am out of ideas.

Thanks in advance!

[David Goodmanson]

Hi Robert,

It appears that the problem is basically that you are subdividing the interval from ent to ext with

t = ent:0.01:ext

and later with

t = ent:0.1:ext

so as to integrate the subintervals in some manner, but your integral always goes over the entire interval, ent to ext.

[Robert Carter]

Thanks I didn't catch the t = ent:0.01:ex so I changed that to t = ent:0.1:ex like the others, I also changed the intervals of integration:

%     
%     if checkedal == false
    length = length + arclength(a,b,c,e,f,g,t,t+0.1);
    radiusdata(u,2*j) = length;
%     checkedal = true;
%     end
    u=u+1;
   end 
    j=j+1;
 end

I think this is what your comment was suggesting. I am still way off unfortunately, I'm getting a total arc length of about 32,000,000 when it should be around 32,000. I'm not sure if it's a coincidence or if i'm just of by a factor of ten. I looked at the length increase between peices and the first peice has a length of 671 which is what it should be, after that the total length almost doubles from one peice to the next.

[Robert Carter]

Maybe because t is in intervals of 0.1 instead of 1 I am off by a factor of 10, and because it's cubic it's off by 10^3. I changed it to this:

length = length + arclength(a,b,c,e,f,g,t,t+0.1);
radiusdata(u,2*j) = length/1000;

I get a total length around what it should be, but now the first section is way off so it must be off by some amount relative to where the section is on the curve or something.

Answer 1

[darova]

I think the problem is in your arclength() function

   for t = ent:0.1:ext
    xd = 3*a*(t-ent)^2 + 2*b*(t-ent)+ c;                 % t = 0:ext-ent
    %%
    length = length + arclength(a,b,c,e,f,g,ent,ext);    % t = ent:ext
    %%
   end 

[Robert Carter]

Thank you! I changed the endpoints for arclength() and it works!

[darova]

YOu can also use polyval() and diff(polyval()) for your purposes