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

9 visualizaciones (últimos 30 días)
SA-W
SA-W el 9 de Mzo. de 2024
Editada: Bruno Luong el 14 de Mzo. de 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 comentarios
Bruno Luong
Bruno Luong el 12 de Mzo. de 2024
@SA-W Because my answers do not address how fn2fm specifically works
SA-W
SA-W el 12 de Mzo. de 2024
You answered my question. If you want to repost it, I will accept.

Iniciar sesión para comentar.

Respuesta aceptada

Bruno Luong
Bruno Luong el 13 de Mzo. de 2024
Editada: Bruno Luong el 14 de Mzo. de 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 comentarios
SA-W
SA-W el 13 de Mzo. de 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 el 13 de Mzo. de 2024
Editada: Bruno Luong el 13 de Mzo. de 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.

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Spline Postprocessing en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by