Questions about SimmSpline

[Matthew Lee]

Dear OpenSim development team and OpenSim users,

Thank you all for your contributions to the biomechanics community! I have some questions about the class SimmSpline, I would appreciate it if anyone could give me some advice.


1. May I ask what is the meaning of SimmSpline() and the mathematical formula behind it? What are the boundary conditions for calculating this type of spline?
2. If the abscissa range of a group of points used to form a SimmSpline is [x1,x2], if an x<x1 is given, what does the value obtained by SimmSpline.calcValue(x) mean?What is the mathematical formula used for this calculation?

I would be very grateful if someone could give me some advice.

Best Wishes,
Matthew

[Ton van den Bogert]

Splines is one of the topics in my numerical methods class this semester, so I was curious about the SimmSpline.

It is a cubic spline, which means that it is a piecewise cubic polynomial; a different polynomial in each interval between two knots.

A spline is smooth, which means that the function, and the 1st and 2nd derivatives are required to be continuous at the knots. SimmSpline is an interpolating spline (as opposed to a smoothing spline), so it is also required to pass exactly through the knots. These requirements are not sufficient to fully specify the polynomial coefficients. Two additional boundary conditions are needed.

The code suggested that SimmSpline estimates the 3rd derivative at the first and last knot, from finite differences on the knot data, and then requires that the first and last polynomial have these 3rd derivative values. I ran some tests and that is indeed the case. I had not seen this type of boundary condition before, but it makes some sense.

More common are the following types of boundary conditions:

• natural: second derivatives at the first and last knot are required to be zero
• clamped: first derivatives at the first and last knot are required to be zero (or some specified non-zero value)
• not-a-knot: third derivative is required to be continuous at the second knot and the second last knot

Equations for the natural spline can be found here: https://en.wikipedia.org/wiki/Spline_interpolation.

The word "spline" comes from ship building where thin pieces of wood were used to draw smooth shapes. The natural spline is the mathematical representation of this. The endpoint second derivative is zero in a wooden spline because this minimizes the elastic energy due to bending. Here are some pictures: https://www.boatdesign.net/threads/spli ... 244/page-2

I ran some tests on a simple example to compare all four splines. The extrapolation is linear in all cases.

Here is the code:

Code: Select all

% define knots
x = [-1 0 1 2 3 4 5 6];
y = [0 0 0 1 1 0 0 0];

% construct the SimmSpline
import org.opensim.modeling.*;
sspline = SimmSpline();
for i = 1:numel(x)
    sspline.addPoint(x(i),y(i));
end

% evaluate the SimmSpline on 1000 points and extrapolate somewhat
xx = linspace(min(x)-1, max(x)+1, 1000)';
yy = zeros(numel(x),4);  % we will store 4 spline curves
for i = 1:numel(xx)
    yy(i,1) = sspline.calcValue(Vector(StdVectorDouble(xx(i))));
end

% use my own spline code, three different boundary conditions
p2 = makespline(x,y,'natural');
yy(:,2) = evalspline(p2,x,y,xx);
p3 = makespline(x,y,'clamped');
yy(:,3) = evalspline(p3,x,y,xx);
p4 = makespline(x,y,'not-a-knot');
yy(:,4) = evalspline(p4,x,y,xx);

% plot them all to compare
subplot(2,2,1)
plot(xx, yy, x,y,'o');
title('knots and splines')
legend('SimmSpline','natural','clamped','not-a-knot', 'data')

% use finite differences to get the derivatives of the splines)
% (we could do this analytically but this is OK)
yd = [NaN(1,4) ; diff(yy)./diff(xx)];
ydd = [diff(yd)./diff(xx) ; NaN(1,4)];
yddd = [diff(ydd)./diff(xx) ; NaN(1,4)];

% plot the derivatives
subplot(2,2,2)
plot(xx,yd);
title('1st derivative')
legend('SimmSpline','natural','clamped','not-a-knot')
subplot(2,2,3)
plot(xx,ydd);
title('2nd derivative')
subplot(2,2,4)
plot(xx,yddd);
title('3rd derivative')
ylim([-10 10])

% this code will confirm that SimmSpline uses 3rd derivative estimated from the data as a boundary condition
yd_data = [diff(y)./diff(x)  NaN];
ydd_data = [NaN  diff(yd_data)./diff(x)];
yddd_data = diff(ydd_data)./diff(x);
fprintf('3rd derivative estimated from first 4 data points: %6.2f\n', yddd_data(2))
fprintf('3rd derivative estimated from last  4 data points: %6.2f\n', yddd_data(end-1))

%=====================================================================
function p = makespline(x,y,bc)
% Makes an interpolating natural spline for data points x,y
%
% Inputs:
%	x		x values for the N knots
%	y		y values for the N knots
%   bc      boundary conditions: 'natural', 'clamped', 'not-a-knot'
%
% Outputs:
%	p		3*(N-1) spline coefficients, stored as [p11 p21 p31 p12 p22 p32 ...]
%           pij is the ith coefficient of the polynomial for the jth
%           time interval
%           The polynomial for time interval j (from x(j) to x(j+1) is:
%           y = y(j) + p1j*(x-x(j)) + p2j*(x-x(j))^2 + p3j*(x-x(j))^3

	N = numel(x);					% determine how many knots there are
	
	% make space to store the system of linear equations A*x = b
	A = zeros(3*(N-1),3*(N-1));		% matrix which will be used to solve p
	b = zeros(3*(N-1),1);			% right hand side
	
	% set up the 3*(N-1) linear equations
    row = 0;
	for i = 1:N-1
		col = 3*(i-1)+1;	% index of p1 for interval x(i) to x(i+1) 
        
		% equation 1: spline must go through x(i+1),y(i+1)
        row = row+1;
		A(row, col) = x(i+1)-x(i);
		A(row, col+1) = (x(i+1)-x(i))^2;
		A(row, col+2) = (x(i+1)-x(i))^3;
		b(row) = y(i+1)-y(i);
		
		% equations 2&3: 1st and 2nd derivative must be continous at x(i)
		if (i > 1)				% we can only do this at node 2 and higher
            row = row+1;
            % continuity of first derivative
			A(row, col-3) = 1;
			A(row, col-2) = 2*(x(i)-x(i-1));
			A(row, col-1) = 3*(x(i)-x(i-1))^2;
			A(row, col)   = -1;
			b(row) = 0;
            
            % continuity of second derivative
            row = row+1;
			A(row, col-2) = 2;
			A(row, col-1) = 6*(x(i)-x(i-1));
			A(row, col+1) = -2;
			b(row) = 0;
        end
        
        % for first or last interval, use boundary conditions specified in
        % the input bc
        if strcmp(bc, 'natural')
		
            % at first interval, second derivative must be zero at x(i)
            if (i == 1)				
                row = row+1;
                A(row, col+1) = 2;		% this means p2 must be zero
                b(row) = 0;
            end

            % at last interval, second derivative must be zero at x(i+1)
            if (i == (N-1))				
                row = row+1;
                A(row, col+1) = 2;		% this means p2+6*p3*(x(i+1)-x(i)) must be zero
                A(row, col+2) = 6*(x(i+1)-x(i));
                b(row) = 0;
            end
            
        elseif strcmp(bc, 'clamped')
            
            % at first interval, first derivative must be zero at x(i)
            if (i == 1)				
                row = row+1;
                % this row says p1 must be zero for this interval
                A(row, col) = 1;		
                b(row) = 0;
            end

            % at last interval, first derivative must be zero at x(i+1)
            if (i == (N-1))				
                row = row+1;
                % this row says p1 + 2*p2*((x(i+1)-x(i))) +
                % 3*p3*(x(i+1)-x(i))^2 must be zero
                A(row, col)   = 1;		
                A(row, col+1) = 2*(x(i+1)-x(i));		
                A(row, col+2) = 3*(x(i+1)-x(i))^2;
                b(row) = 0;
            end
           
        elseif strcmp(bc, 'not-a-knot')
            
            % third derivative must be continuous at x(2)
            if (i == 1)				
                row = row+1;
                % this row says 6*p3 = 6*p3(from next interval)		
                A(row, col+2) = 1;		
                A(row, col+5) = -1;
                b(row) = 0;
            end

            % third derivative must be continuous at x(N-1)
            if (i == (N-1))				
                row = row+1;
                % this row says 6*p3 = 6*p3(from next interval)		
                A(row, col-1) = 1;		
                A(row, col+2) = -1;
                b(row) = 0;
            end
           
        end
	end
	
	% solve the polynomial coefficients
	p = A\b;		

end
%=======================================================================
function j = locate(x,xq)
    % find the integer j, such that x(j) <= xq < x(j+1)
    %
    % Inputs:
    %   x............Array of knots (must be increasing)
    %   xq...........Query point
    %
    % Output:
    %   j............Index of knot interval where xq is found

    % check that x values are increasing
    if min(diff(x)) < 0
        error('locate: knot array x is not increasing');
    end
    
    % check to see if xq is within the range defined by x
    N = numel(x);
    if xq < x(1)
        j = 0;
    elseif xq > x(N)
       j = N;
    else    
        % find the two surrounding knots by bisection
        ilow = 1;
        ihigh = N;
        while (ihigh - ilow) > 1
            imid = round((ihigh+ilow)/2);
            if xq < x(imid)
                ihigh = imid;
            else
                ilow = imid;
            end
        end
        % we now know that xq is between x(ilow) and x(ilow+1)
        j = ilow;
    end
end
%=========================================================================
function yq = evalspline(p,x,y,xq)
	% Evaluates a cubic spline at a series of query points xq
	%
	% Inputs:
	%	p		3*(N-1) spline coefficients, stored as [p11 p21 p31 p12 p22 p32 ...]
    %           pij is the ith coefficient of the polynomial for the jth
    %           time interval
    %           The polynomial for time interval j (from x(j) to x(j+1) is:
    %           y = y(j) + p1j*(x-x(j)) + p2j*(x-x(j))^2 + p3j*(x-x(j))^3
	%   x		the N x-values at the knots (must be increasing)
	%   y		the N y-values at the knots
	%	xq		a list of query points where the spline must be evaluated
	%
	% Outputs:
	%	yq		spline function values, evaluated at all elements of xq
	
    N = numel(x);       % the number of knots
	Nq = numel(xq);		% determine how many elements in xq
	yq = zeros(size(xq));	% create space, same size as xq, to store yq
	
	% go through all elements of xq, and compute yq for each of them
	for iq = 1:Nq
		% find the knot x(j) that is just below xq(iq)
		j = locate(x,xq(iq));

        % if j is outside of the intervals 1..N-1, do linear extrapolation
        if j==0
            slope = p(1);
		    yq(iq) = y(1) + slope*(xq(iq)-x(1));
        elseif j==N
            slope = p(3*(N-2)+1) + 2*p(3*(N-2)+2)*(x(N)-x(N-1)) + 3*p(3*(N-2)+3)*(x(N)-x(N-1))^2;
		    yq(iq) = y(N) + slope*(xq(iq)-x(N));
        else
            % use the cubic polynomial for interval j
		    % from p, extract the 3 polynomial coefficients for the 
		    % interval x(i) to x(i+1)
		    p1 = p(3*(j-1) + 1);
		    p2 = p(3*(j-1) + 2);
		    p3 = p(3*(j-1) + 3);
		    
		    % compute the polynomial value at xq(iq)
		    t = xq(iq) - x(j);
		    yq(iq) = y(j) + p1*t + p2*t^2 + p3*t^3;
        end
	end

end

And the results:

untitled.png (44.11 KiB) Viewed 887 times

The last plot shows that the SimmSpline has 3rd derivative values of +1 and -1 at the ends, which are the same as the result of a finite difference approximation on the knot data.

Conclusions:

• if you want to extrapolate, it's probably best to use a natural spline
• if you want to use the spline in the first or last knot interval, natural spline is safest also
• elsewhere, it does not matter much which spline is used

In OpenSim, SimmSpline is typically used to define a relationship between two variables. In OpenSim models, you often see extra points added before or after the region of interest. This is to reduce the influence of the boundary conditions, as the examples show.

[Matthew Lee]

Dear Professor Ton,

Thank you very much for your answer! It is a great honor for me to receive your advice!Your answer is wonderful and detailed, which is much more exciting than the textbooks I have read! From studying your answer, I believe I have a pretty clear understanding of cubic splines, thank you!

By the way, for the OpenSim developers, I seem to have found a small bug about the extrapolation slope of the SimmSpline at the right endpoint (Usually this probably won't affect anything). On line 538 of SimmSpline.cpp,https://github.com/opensim-org/opensim- ... e.cpp#L538

The linear extrapolation slope should be

Code: Select all

return _y[n-1] + (aX - _x[n-1])*(_b[n-1]+ 2*_c[n-1] *(aX - _x[n-1])+3*_d[n-1]*(aX - _x[n-1])* (aX - _x[n-1]));.

Best Wishes,
Matthew

[Ton van den Bogert]

You're welcome!

I am not sure it's a bug. The extrapolation always looks correct in my tests with SimmSpline.

I have not thoroughly analyzed the code, but keep in mind that an array in C++ starts with index 0. With n knots, there are n-1 knot intervals, indexed from 0 to n-2. The polynomial with index (n-1) is the extrapolation, not the final knot interval. The first derivative at the final knot is then b[n-1], the coefficient of the linear term in that polynomial.

[Matthew Lee]

Dear Professor Ton,

Thank you very much for your reply! I think you are right, the code that I mentioned before is actually correct, glad that with your help I figured out this confusion!

I still have a few more questions and still need your help. After running inverse kinematics, what we actually get is joint angle data at several discrete time points (assuming these points belong to set A). If I want to obtain joint angle data at a time point that is not belong to set A, interpolation is required to obtain it. May I ask which type of spline curve is used by default in OpenSim to interpolate the results of inverse kinematics?

Since the mot file output by inverse kinematics only contains coordinate (joint angle) information,if I want to further obtain joint angular velocity and joint angular acceleration information, is there an easy way in OpenSim to do this directly?

I really appreciate your professional answers and help, thanks again for your patience and time!

Best Wishes,
Matthew