octave_packages/nurbs-1.3.6/curvederivcpts.m

   1 function pk = curvederivcpts (n, p, U, P, d, r1, r2)
   2 % Compute control points of n-th derivatives of a B-spline curve.
   3 %
   4 % usage: pk = curvederivcpts (n, p, U, P, d)
   5 %        pk = curvederivcpts (n, p, U, P, d, r1, r2)
   6 %
   7 % If r1, r2 are not given, all the control points are computed.
   8 %
   9 %  INPUT:
  10 %         n+1 = number of control points
  11 %         p   = degree of the spline
  12 %         d   = maximum derivative order (d<=p)
  13 %         U   = knots
  14 %         P   = control points
  15 %         r1  = first control point to compute
  16 %         r2  = auxiliary index for the last control point to compute
  17 %  OUTPUT:
  18 %         pk(k,i) = i-th control point of (k-1)-th derivative, r1 <= i <= r2-k
  19 %
  20 % Adaptation of algorithm A3.3 from the NURBS book, pg98.
  21 %
  22 %    Copyright (C) 2009 Carlo de Falco
  23 %
  24 %    This program is free software: you can redistribute it and/or modify
  25 %    it under the terms of the GNU General Public License as published by
  26 %    the Free Software Foundation, either version 2 of the License, or
  27 %    (at your option) any later version.
  28
  29 %    This program is distributed in the hope that it will be useful,
  30 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  31 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  32 %    GNU General Public License for more details.
  33 %
  34 %    You should have received a copy of the GNU General Public License
  35 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  36
  37   if (nargin <= 5)
  38     r1 = 0;
  39     r2 = n;
  40   end
  41
  42   r = r2 - r1;
  43   for i=0:r
  44     pk(1, i+1) = P(r1+i+1);
  45   end
  46
  47   for k=1:d
  48     tmp = p - k + 1;
  49     for i=0:r-k
  50       pk (k+1, i+1) = tmp * (pk(k,i+2)-pk(k,i+1)) / ...
  51           (U(r1+i+p+2)-U(r1+i+k+1));
  52     end
  53   end
  54 end
  55
  56 %!test
  57 %! line = nrbmak([0.0 1.5; 0.0 3.0],[0.0 0.0 1.0 1.0]);
  58 %! pk   = curvederivcpts (line.number-1, line.order-1, line.knots,...
  59 %!                        line.coefs(1,:), 2);
  60 %! assert (pk, [0 3/2; 3/2 0], 100*eps);