Coefficents of piecewise polynomial in matlab

Hi I have a simple example of piecewise polynomial from MATLAB using the function 'cscvn'

cscvn( [1 0 -1    0 1;0 1 0   -1 0] )

ans : form: 'pp'

breaks: [0 1.1892 2.3784 3.5676 4.7568]

coefs: [8×4 double]

pieces: 4

order: 4

dim: 2

My question is : How do I read the coeffiecnt matrix as shown above. I can understand from the above coefs matrix that it has 4 columns corresponding to 4 coefficents of a 4th order polynomial. But why this matrix has 8 columns? I only have 4 pieces of polynomial, then why double the number of rows in coefs matrix?

Seocnd question is : How can I use the above coefs matrix to get the equation of piecewise poilynomial

Answer 1

[John D'Errico]

Do you understand the result will not be some simple polynomial like you are used to seeing?

plot([1 0 -1    0 1],[0 1 0   -1 0],'-o')

You want a curve that follows around this diamond shaped thing. So you CANNOT have a function of the form y(x). PERIOD. You cannot do so.

Instead, cscvn generates a curve as a function of a parameter, I'll call it t, where t is the arclength of the piecewise linear curve that connects those dots. So cscvn creates the pair of functions x(t), and y(t).

READ THE DOCS FOR CSCVN: cscvn

In there this scheme is explained fairly clearly.

Again, there is no simple equation of the form y(x). There cannot be, as if you look at the plot, the result cannot be a single valued function.

[Mukesh Singh Bisht]

Yes, I understand that cscvn will create the pair of functions x(t), and y(t). But I want to know in the coefs matricx, which rows correspond to coefficients for x(t) and which rows for y(t)?

[Torsten]

Does this help ?

c=cscvn( [1 0 -1    0 1;0 1 0   -1 0] );
c.coefs

ans = 8×4

0.2973 -1.0607 0.0000 1.0000 -0.2973 0.0000 1.2613 0 0.2973 0.0000 -1.2613 0 0.2973 -1.0607 -0.0000 1.0000 -0.2973 1.0607 -0.0000 -1.0000 0.2973 0 -1.2613 0 -0.2973 0.0000 1.2613 0 -0.2973 1.0607 0.0000 -1.0000

c.breaks

ans = 1×5

0 1.1892 2.3784 3.5676 4.7568

x = @(t)1-1.0607*t.^2+0.2973*t.^3;
y = @(t)1.2613*t-0.2973*t.^3;
t = linspace(0,1.1892,20);
hold on
plot(x(t),y(t))
x = @(t)-1.2613*(t-1.1892)+0.2973*(t-1.1892).^3;
y = @(t)1-1.0607*(t-1.1892).^2+0.2973*(t-1.1892).^3;
t = linspace(1.1892,2.3784,20);
plot(x(t),y(t))
x = @(t)-1+1.0607*(t-2.3784).^2-0.2973*(t-2.3784).^3;
y = @(t)-1.2613*(t-2.3784)+0.2973*(t-2.3784).^3;
t = linspace(2.3784,3.5676,20);
plot(x(t),y(t))
x = @(t)1.2613*(t-3.5676)-0.2973*(t-3.5676).^3;
y = @(t)-1+1.0607*(t-3.5676).^2-0.2973*(t-3.5676).^3;
t = linspace(3.5676,4.7568,20);
plot(x(t),y(t))
hold off

[Walter Roberson]

(The rows appear to alternate, one row of x, then one row of y, with the pairs repeated according to the number of pieces. So x coordinates would be (1:2:end,:) and y coordinates would be (2:2:end,:)

[Bruno Luong]

@John D'Errico "Instead, cscvn generates a curve as a function of a parameter, I'll call it t, where t is the arclength of the piecewise linear curve that connects those dots"

Actually it's the sum of the square-root of the linear arclength, this thing is called Eugene Lee centripetal parametrization according to cscvn doc

Answer 2

[Bruno Luong]

Similar answer as here

points = [1 0 -1    0 1;0 1 0   -1 0];

pp = cscvn(points)

pp = struct with fields:
      form: 'pp'
    breaks: [0 1.1892 2.3784 3.5676 4.7568]
     coefs: [8×4 double]
    pieces: 4
     order: 4
       dim: 2

m = pp.dim;
p = pp.pieces;
%D = diff(points,1,2);
%d = sqrt(sqrt(sum(D.^2,1)));
%t = cumsum([0 d]); % Eugene Lee's, equal to pp.breaks, see https://www.mathworks.com/help/curvefit/cscvn.html
t = pp.breaks;

ni = 200;
ti = linspace(t(1),t(end), ni);

loc = discretize(ti, t);
x = zeros(length(ti), m);
coefs = reshape(pp.coefs, m, p, []);
for k=1:p
    ink = loc == k;
    tik = ti(ink);
    dti = tik - t(k);
    for d=1:m
        P = coefs(d, k, :);
        x(ink,d) = polyval(P(:), dti);
    end
end

close all
hold on
if m == 2
    plot(points(1,:),points(2,:),'or'),
    plot(x(:,1), x(:,2),'b.')
elseif m == 3
    plot3(points(1,:),points(2,:),points(3,:),'or'),
    plot3(x(:,1), x(:,2), x(:,3), 'b.')
else
    warning('plotting not supported')
end