steps to convert spline from B-form to pp-form (fn2fm)

11 vues (au cours des 30 derniers jours)
SA-W
SA-W le 9 Mar 2024
Modifié(e) : Bruno Luong le 14 Mar 2024
x = [3.0,4.5,6.0,7.5,9.0,12.0,15.0];
y = [0 0.0343653 0.0694232 0.105143 0.141178 0.246013 0.630537];
f_bm = spapi(5, x, y);
f_pp = fn2fm(f_bm, 'pp');
Based on f_bm.knots, it is easy to compute f_pp.breaks. But how is f_pp.coefs calculated? Is there a linear system solved or is there an analytical equation for computing the coefficients based on the control points + knots + degree?
Thank you!
  12 commentaires
Bruno Luong
Bruno Luong le 12 Mar 2024
@SA-W Because my answers do not address how fn2fm specifically works
SA-W
SA-W le 12 Mar 2024
You answered my question. If you want to repost it, I will accept.

Connectez-vous pour commenter.

Réponse acceptée

Bruno Luong
Bruno Luong le 13 Mar 2024
Modifié(e) : Bruno Luong le 14 Mar 2024
The pp-form stores the polynomial of each each subiterval, the variable is x := (t-ti) where ti is the left knot of interval #i. https://www.mathworks.com/help/curvefit/the-ppform.html
there is a function private/Bspline2pp.m that just does the conversion to pp
It evalutes the spline in k points inside each subitervals then invert the Vandermond matrix to compute a row of pp.coefs. The last stage can be done with polyfit for coding convenience.
You can also compute from Taylor expansion of the polynomial wrt th the left knot as Torsen has suggested.
It requires to compute 0 to (k-1) order derivatives
Code to avoid using fn2fm (why?!) if this idea
WARNING this only works for non-duplicated knot vectors
x = [3.0,4.5,6.0,7.5,9.0,12.0,15.0];
y = [0 0.0343653 0.0694232 0.105143 0.141178 0.246013 0.630537];
f_bm = spapi(5, x, y);
f_pp = fn2fm(f_bm, 'pp'); % just used for checking correctness
k = f_bm.order;
breaks = f_bm.knots(k:end-k+1);
pieces = length(breaks)-1;
coefs = zeros(pieces,k);
Sd = f_bm;
xi = breaks(1:pieces);
for j=0:k-1
coefs(:,k-j) = fnval(Sd, xi)./factorial(j);
Sd = fnder(Sd);
end
pp = struct('form', 'pp', ...
'breaks', breaks, ...
'coefs', coefs,...
'piecs', pieces, ...
'order', k, ...
'dim', 1);
% Check the result
format long
f_pp.coefs
ans = 3×5
-0.000010839458734 0.000095398958747 -0.000104662728187 0.022889129608328 -0.000000000000000 0.000030128251860 -0.000067192922266 0.000053996227019 0.023842354392321 0.087249675574952 0.000176086472898 0.000158768966683 0.000311553851941 0.024130562157454 0.132073372169513
pp.coefs
ans = 3×5
-0.000010839458734 0.000095398958747 -0.000104662728187 0.022889129608328 -0.000000000000000 0.000030128251860 -0.000067192922266 0.000053996227019 0.023842354392321 0.087249675574952 0.000176086472898 0.000158768966683 0.000311553851941 0.024130562157453 0.132073372169513
% Note different last digit in coefs(3,4)
  3 commentaires
SA-W
SA-W le 13 Mar 2024
Thanks a lot.
Actually, my motivtion to switch to pp-form is that the evaluation of the spline and its derivatives is faster than in B-form because the B-spline basis functions are defined recursively. This is of course not very efficient in a computer program.
But I found De Boor's algorithm, see e.g wiki, which only needs a for loop to evaluate the spline; no recursion anymore. The algorithm showed at the bottom can be easily adjusted to also compute derivatives up to nth order.
Bruno Luong
Bruno Luong le 13 Mar 2024
Modifié(e) : Bruno Luong le 13 Mar 2024
I'm not sure, the De Boor algorithm and B-spline basis function evaluation are usually both based on Cox–de Boor recursion formula. It is just implement with different workfow. We never need to compute the explitly the basis functions if spline values are required to be calculated.
I suppose (not sure) when you call fnval with B-spline form input, it would use De Boor algorithm.
In my BSFK toolbox it is implement in function private/Berntein.m with vectorasation, if you want to dig in and see how it is implemented in MATLAB.
To me the most efficient evaluation is with pp-form. This is the whole point why this form is invented.

Connectez-vous pour commenter.

Plus de réponses (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by