octave_packages/nurbs-1.3.6/basisfunder.m

   1 function dersv = basisfunder (ii, pl, uu, u_knotl, nders)
   2
   3 % BASISFUNDER:  B-Spline Basis function derivatives.
   4 %
   5 % Calling Sequence:
   6 %
   7 %   ders = basisfunder (ii, pl, uu, k, nd)
   8 %
   9 %    INPUT:
  10 %
  11 %      ii  - knot span index (see findspan)
  12 %      pl  - degree of curve
  13 %      uu  - parametric points
  14 %      k   - knot vector
  15 %      nd  - number of derivatives to compute
  16 %
  17 %    OUTPUT:
  18 %
  19 %      ders - ders(n, i, :) (i-1)-th derivative at n-th point
  20 %
  21 %    Adapted from Algorithm A2.3 from 'The NURBS BOOK' pg72.
  22 %
  23 %    Copyright (C) 2009,2011 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   dersv = zeros(numel(uu), nders+1, pl+1);
  39
  40   for jj = 1:numel(uu)
  41
  42     i = ii(jj)+1; %% convert to base-1 numbering of knot spans
  43     u = uu(jj);
  44
  45     ders = zeros(nders+1,pl+1);
  46     ndu = zeros(pl+1,pl+1);
  47     left = zeros(pl+1);
  48     right = zeros(pl+1);
  49     a = zeros(2,pl+1);
  50     ndu(1,1) = 1;
  51     for j = 1:pl
  52       left(j+1) = u - u_knotl(i+1-j);
  53       right(j+1) = u_knotl(i+j) - u;
  54       saved = 0;
  55       for r = 0:j-1
  56         ndu(j+1,r+1) = right(r+2) + left(j-r+1);
  57         temp = ndu(r+1,j)/ndu(j+1,r+1);
  58         ndu(r+1,j+1) = saved + right(r+2)*temp;
  59         saved = left(j-r+1)*temp;
  60       end
  61       ndu(j+1,j+1) = saved;
  62     end
  63     for j = 0:pl
  64       ders(1,j+1) = ndu(j+1,pl+1);
  65     end
  66     for r = 0:pl
  67       s1 = 0;
  68       s2 = 1;
  69       a(1,1) = 1;
  70       for k = 1:nders %compute kth derivative
  71         d = 0;
  72         rk = r-k;
  73         pk = pl-k;
  74         if (r >= k)
  75           a(s2+1,1) = a(s1+1,1)/ndu(pk+2,rk+1);
  76           d = a(s2+1,1)*ndu(rk+1,pk+1);
  77         end
  78         if (rk >= -1)
  79           j1 = 1;
  80         else
  81           j1 = -rk;
  82         end
  83         if ((r-1) <= pk)
  84           j2 = k-1;
  85         else
  86           j2 = pl-r;
  87         end
  88         for j = j1:j2
  89           a(s2+1,j+1) = (a(s1+1,j+1) - a(s1+1,j))/ndu(pk+2,rk+j+1);
  90           d = d + a(s2+1,j+1)*ndu(rk+j+1,pk+1);
  91         end
  92         if (r <= pk)
  93           a(s2+1,k+1) = -a(s1+1,k)/ndu(pk+2,r+1);
  94           d = d + a(s2+1,k+1)*ndu(r+1,pk+1);
  95         end
  96         ders(k+1,r+1) = d;
  97         j = s1;
  98         s1 = s2;
  99         s2 = j;
 100       end
 101     end
 102     r = pl;
 103     for k = 1:nders
 104       for j = 0:pl
 105         ders(k+1,j+1) = ders(k+1,j+1)*r;
 106       end
 107       r = r*(pl-k);
 108     end
 109
 110     dersv(jj, :, :) = ders;
 111
 112   end
 113
 114 end
 115
 116 %!test
 117 %! k    = [0 0 0 0 1 1 1 1];
 118 %! p    = 3;
 119 %! u    = rand (1);
 120 %! i    = findspan (numel(k)-p-2, p, u, k);
 121 %! ders = basisfunder (i, p, u, k, 1);
 122 %! sumders = sum (squeeze(ders), 2);
 123 %! assert (sumders(1), 1, 1e-15);
 124 %! assert (sumders(2:end), 0, 1e-15);
 125
 126 %!test
 127 %! k    = [0 0 0 0 1/3 2/3 1 1 1 1];
 128 %! p    = 3;
 129 %! u    = rand (1);
 130 %! i    = findspan (numel(k)-p-2, p, u, k);
 131 %! ders = basisfunder (i, p, u, k, 7);
 132 %! sumders = sum (squeeze(ders), 2);
 133 %! assert (sumders(1), 1, 1e-15);
 134 %! assert (sumders(2:end), zeros(rows(squeeze(ders))-1, 1), 1e-13);
 135
 136 %!test
 137 %! k    = [0 0 0 0 1/3 2/3 1 1 1 1];
 138 %! p    = 3;
 139 %! u    = rand (100, 1);
 140 %! i    = findspan (numel(k)-p-2, p, u, k);
 141 %! ders = basisfunder (i, p, u, k, 7);
 142 %! for ii=1:10
 143 %!   sumders = sum (squeeze(ders(ii,:,:)), 2);
 144 %!   assert (sumders(1), 1, 1e-15);
 145 %!   assert (sumders(2:end), zeros(rows(squeeze(ders(ii,:,:)))-1, 1), 1e-13);
 146 %! end
 147 %! assert (ders(:, (p+2):end, :), zeros(numel(u), 8-p-1, p+1), 1e-13)
 148 %! assert (all(all(ders(:, 1, :) <= 1)), true)
 149
 150
 151