octave_packages/nurbs-1.3.6/nrbcrvderiveval.m

   1 function ck = nrbcrvderiveval (crv, u, d)
   2
   3 %
   4 % NRBCRVDERIVEVAL: Evaluate n-th order derivatives of a NURBS curve.
   5 %
   6 % usage: skl = nrbcrvderiveval (crv, u, d)
   7 %
   8 %   INPUT:
   9 %
  10 %   crv : NURBS curve structure, see nrbmak
  11 %
  12 %   u   : parametric coordinate of the points where we compute the derivatives
  13 %
  14 %   d   : number of partial derivatives to compute
  15 %
  16 %
  17 %   OUTPUT:
  18 %
  19 %   ck (i, j, l) = i-th component derived j-1 times at the l-th point.
  20 %
  21 % Adaptation of algorithm A4.2 from the NURBS book, pg127
  22 %
  23 %    Copyright (C) 2010 Carlo de Falco, Rafael Vazquez
  24 %
  25 %    This program is free software: you can redistribute it and/or modify
  26 %    it under the terms of the GNU General Public License as published by
  27 %    the Free Software Foundation, either version 2 of the License, or
  28 %    (at your option) any later version.
  29
  30 %    This program is distributed in the hope that it will be useful,
  31 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  32 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  33 %    GNU General Public License for more details.
  34 %
  35 %    You should have received a copy of the GNU General Public License
  36 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  37
  38  ck = zeros (3, d+1, numel(u));
  39
  40  for iu = 1:numel(u);
  41    wders = squeeze (curvederiveval (crv.number-1, crv.order-1,  ...
  42                                      crv.knots, squeeze (crv.coefs(4, :)), u(iu), d));
  43
  44    for idim = 1:3
  45      Aders = squeeze (curvederiveval (crv.number-1, crv.order-1,  ...
  46                         crv.knots, squeeze (crv.coefs(idim, :)), u(iu), d));
  47      for k=0:d
  48        v = Aders(k+1);
  49        for i=1:k
  50              v = v - nchoosek(k,i)*wders(i+1)*ck(idim, k-i+1, iu);
  51        end
  52            ck(idim, k+1, iu) = v/wders(1);
  53      end
  54    end
  55  end
  56 end
  57
  58 %!test
  59 %! knots = [0 0 0 1 1 1];
  60 %! coefs(:,1) = [0; 0; 0; 1];
  61 %! coefs(:,2) = [1; 0; 1; 1];
  62 %! coefs(:,3) = [1; 1; 1; 2];
  63 %! crv = nrbmak (coefs, knots);
  64 %! u = linspace (0, 1, 10);
  65 %! ck = nrbcrvderiveval (crv, u, 2);
  66 %! w  = @(x) 1 + x.^2;
  67 %! dw = @(x) 2*x;
  68 %! F1 = @(x) (2*x - x.^2)./w(x);
  69 %! F2 = @(x) x.^2./w(x);
  70 %! F3 = @(x) (2*x - x.^2)./w(x);
  71 %! dF1 = @(x) (2 - 2*x)./w(x) - 2*(2*x - x.^2).*x./w(x).^2;
  72 %! dF2 = @(x) 2*x./w(x) - 2*x.^3./w(x).^2;
  73 %! dF3 = @(x) (2 - 2*x)./w(x) - 2*(2*x - x.^2).*x./w(x).^2;
  74 %! d2F1 = @(x) -2./w(x) - 2*x.*(2-2*x)./w(x).^2 - (8*x-6*x.^2)./w(x).^2 + 8*x.^2.*(2*x-x.^2)./w(x).^3;
  75 %! d2F2 = @(x) 2./w(x) - 4*x.^2./w(x).^2 - 6*x.^2./w(x).^2 + 8*x.^4./w(x).^3;
  76 %! d2F3 = @(x) -2./w(x) - 2*x.*(2-2*x)./w(x).^2 - (8*x-6*x.^2)./w(x).^2 + 8*x.^2.*(2*x-x.^2)./w(x).^3;
  77 %! assert ([F1(u); F2(u); F3(u)], squeeze(ck(:, 1, :)), 1e2*eps);
  78 %! assert ([dF1(u); dF2(u); dF3(u)], squeeze(ck(:, 2, :)), 1e2*eps);
  79 %! assert ([d2F1(u); d2F2(u); d2F3(u)], squeeze(ck(:, 3, :)), 1e2*eps);