Mercurial > forge
view extra/nurbs/inst/basisfunder.m @ 12662:2a027badd794 octave-forge
Added a message id in the warning
author | rafavzqz |
---|---|
date | Tue, 07 Jul 2015 13:52:48 +0000 |
parents | 37d08939bb7b |
children |
line wrap: on
line source
function dersv = basisfunder (ii, pl, uu, u_knotl, nders) % BASISFUNDER: B-Spline Basis function derivatives. % % Calling Sequence: % % ders = basisfunder (ii, pl, uu, k, nd) % % INPUT: % % ii - knot span index (see findspan) % pl - degree of curve % uu - parametric points % k - knot vector % nd - number of derivatives to compute % % OUTPUT: % % ders - ders(n, i, :) (i-1)-th derivative at n-th point % % Adapted from Algorithm A2.3 from 'The NURBS BOOK' pg72. % % See also: % % numbasisfun, basisfun, findspan % % Copyright (C) 2009,2011 Rafael Vazquez % % This program is free software: you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by % the Free Software Foundation, either version 3 of the License, or % (at your option) any later version. % This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with this program. If not, see <http://www.gnu.org/licenses/>. dersv = zeros(numel(uu), nders+1, pl+1); for jj = 1:numel(uu) i = ii(jj)+1; %% convert to base-1 numbering of knot spans u = uu(jj); ders = zeros(nders+1,pl+1); ndu = zeros(pl+1,pl+1); left = zeros(pl+1); right = zeros(pl+1); a = zeros(2,pl+1); ndu(1,1) = 1; for j = 1:pl left(j+1) = u - u_knotl(i+1-j); right(j+1) = u_knotl(i+j) - u; saved = 0; for r = 0:j-1 ndu(j+1,r+1) = right(r+2) + left(j-r+1); temp = ndu(r+1,j)/ndu(j+1,r+1); ndu(r+1,j+1) = saved + right(r+2)*temp; saved = left(j-r+1)*temp; end ndu(j+1,j+1) = saved; end for j = 0:pl ders(1,j+1) = ndu(j+1,pl+1); end for r = 0:pl s1 = 0; s2 = 1; a(1,1) = 1; for k = 1:nders %compute kth derivative d = 0; rk = r-k; pk = pl-k; if (r >= k) a(s2+1,1) = a(s1+1,1)/ndu(pk+2,rk+1); d = a(s2+1,1)*ndu(rk+1,pk+1); end if (rk >= -1) j1 = 1; else j1 = -rk; end if ((r-1) <= pk) j2 = k-1; else j2 = pl-r; end for j = j1:j2 a(s2+1,j+1) = (a(s1+1,j+1) - a(s1+1,j))/ndu(pk+2,rk+j+1); d = d + a(s2+1,j+1)*ndu(rk+j+1,pk+1); end if (r <= pk) a(s2+1,k+1) = -a(s1+1,k)/ndu(pk+2,r+1); d = d + a(s2+1,k+1)*ndu(r+1,pk+1); end ders(k+1,r+1) = d; j = s1; s1 = s2; s2 = j; end end r = pl; for k = 1:nders for j = 0:pl ders(k+1,j+1) = ders(k+1,j+1)*r; end r = r*(pl-k); end dersv(jj, :, :) = ders; end end %!test %! k = [0 0 0 0 1 1 1 1]; %! p = 3; %! u = rand (1); %! i = findspan (numel(k)-p-2, p, u, k); %! ders = basisfunder (i, p, u, k, 1); %! sumders = sum (squeeze(ders), 2); %! assert (sumders(1), 1, 1e-15); %! assert (sumders(2:end), 0, 1e-15); %!test %! k = [0 0 0 0 1/3 2/3 1 1 1 1]; %! p = 3; %! u = rand (1); %! i = findspan (numel(k)-p-2, p, u, k); %! ders = basisfunder (i, p, u, k, 7); %! sumders = sum (squeeze(ders), 2); %! assert (sumders(1), 1, 1e-15); %! assert (sumders(2:end), zeros(rows(squeeze(ders))-1, 1), 1e-13); %!test %! k = [0 0 0 0 1/3 2/3 1 1 1 1]; %! p = 3; %! u = rand (100, 1); %! i = findspan (numel(k)-p-2, p, u, k); %! ders = basisfunder (i, p, u, k, 7); %! for ii=1:10 %! sumders = sum (squeeze(ders(ii,:,:)), 2); %! assert (sumders(1), 1, 1e-15); %! assert (sumders(2:end), zeros(rows(squeeze(ders(ii,:,:)))-1, 1), 1e-13); %! end %! assert (ders(:, (p+2):end, :), zeros(numel(u), 8-p-1, p+1), 1e-13) %! assert (all(all(ders(:, 1, :) <= 1)), true)