octave_packages/nurbs-1.3.6/nrbbasisfun.m

   1 function [B, id] = nrbbasisfun (points, nrb)
   2
   3 % NRBBASISFUN: Basis functions for NURBS
   4 %
   5 % Calling Sequence:
   6 %
   7 %    B     = nrbbasisfun (u, crv)
   8 %    B     = nrbbasisfun ({u, v}, srf)
   9 %   [B, N] = nrbbasisfun ({u, v}, srf)
  10 %   [B, N] = nrbbasisfun (p, srf)
  11 %
  12 %    INPUT:
  13 %
  14 %      u or p(1,:,:)  - parametric points along u direction
  15 %      v or p(2,:,:)  - parametric points along v direction
  16 %      crv - NURBS curve
  17 %      srf - NURBS surface
  18 %
  19 %    OUTPUT:
  20 %
  21 %      B - Value of the basis functions at the points
  22 %          size(B)=[numel(u),(p+1)] for curves
  23 %          or [numel(u)*numel(v), (p+1)*(q+1)] for surfaces
  24 %
  25 %      N - Indices of the basis functions that are nonvanishing at each
  26 %          point. size(N) == size(B)
  27 %
  28 %
  29 %    Copyright (C) 2009 Carlo de Falco
  30 %
  31 %    This program is free software: you can redistribute it and/or modify
  32 %    it under the terms of the GNU General Public License as published by
  33 %    the Free Software Foundation, either version 2 of the License, or
  34 %    (at your option) any later version.
  35
  36 %    This program is distributed in the hope that it will be useful,
  37 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  38 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  39 %    GNU General Public License for more details.
  40 %
  41 %    You should have received a copy of the GNU General Public License
  42 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  43
  44   if (   (nargin<2) ...
  45       || (nargout>2) ...
  46       || (~isstruct(nrb)) ...
  47       || (iscell(points) && ~iscell(nrb.knots)) ...
  48       || (~iscell(points) && iscell(nrb.knots) && (size(points,1)~=2)) ...
  49       || (~iscell(nrb.knots) && (nargout>1)) ...
  50       )
  51     error('Incorrect input arguments in nrbbasisfun');
  52   end
  53
  54   if (~iscell(nrb.knots))         %% NURBS curve
  55
  56     [B, id] = nrb_crv_basisfun__ (points, nrb);
  57
  58   elseif size(nrb.knots,2) == 2 %% NURBS surface
  59     if (iscell(points))
  60       [v, u] = meshgrid(points{2}, points{1});
  61       p = [u(:), v(:)]';
  62     else
  63       p = points;
  64     end
  65
  66     [B, id] = nrb_srf_basisfun__ (p, nrb);
  67
  68   else                            %% NURBS volume
  69     error('The function nrbbasisfun is not yet ready for volumes')
  70   end
  71 end
  72
  73 %!demo
  74 %! U = [0 0 0 0 1 1 1 1];
  75 %! x = [0 1/3 2/3 1] ;
  76 %! y = [0 0 0 0];
  77 %! w = [1 1 1 1];
  78 %! nrb = nrbmak ([x;y;y;w], U);
  79 %! u = linspace(0, 1, 30);
  80 %! B = nrbbasisfun (u, nrb);
  81 %! xplot = sum(bsxfun(@(x,y) x.*y, B, x),2);
  82 %! plot(xplot, B)
  83 %! title('Cubic Bernstein polynomials')
  84 %! hold off
  85
  86 %!test
  87 %! U = [0 0 0 0 1 1 1 1];
  88 %! x = [0 1/3 2/3 1] ;
  89 %! y = [0 0 0 0];
  90 %! w = rand(1,4);
  91 %! nrb = nrbmak ([x;y;y;w], U);
  92 %! u = linspace(0, 1, 30);
  93 %! B = nrbbasisfun (u, nrb);
  94 %! xplot = sum(bsxfun(@(x,y) x.*y, B, x),2);
  95 %!
  96 %! yy = y; yy(1) = 1;
  97 %! nrb2 = nrbmak ([x.*w;yy;y;w], U);
  98 %! aux = nrbeval(nrb2,u);
  99 %! %figure, plot(xplot, B(:,1), aux(1,:).', w(1)*aux(2,:).')
 100 %! assert(B(:,1), w(1)*aux(2,:).', 1e-6)
 101 %!
 102 %! yy = y; yy(2) = 1;
 103 %! nrb2 = nrbmak ([x.*w;yy;y;w], U);
 104 %! aux = nrbeval(nrb2, u);
 105 %! %figure, plot(xplot, B(:,2), aux(1,:).', w(2)*aux(2,:).')
 106 %! assert(B(:,2), w(2)*aux(2,:).', 1e-6)
 107 %!
 108 %! yy = y; yy(3) = 1;
 109 %! nrb2 = nrbmak ([x.*w;yy;y;w], U);
 110 %! aux = nrbeval(nrb2,u);
 111 %! %figure, plot(xplot, B(:,3), aux(1,:).', w(3)*aux(2,:).')
 112 %! assert(B(:,3), w(3)*aux(2,:).', 1e-6)
 113 %!
 114 %! yy = y; yy(4) = 1;
 115 %! nrb2 = nrbmak ([x.*w;yy;y;w], U);
 116 %! aux = nrbeval(nrb2,u);
 117 %! %figure, plot(xplot, B(:,4), aux(1,:).', w(4)*aux(2,:).')
 118 %! assert(B(:,4), w(4)*aux(2,:).', 1e-6)
 119
 120 %!test
 121 %! p = 2;   q = 3;   m = 4; n = 5;
 122 %! Lx  = 1; Ly  = 1;
 123 %! nrb = nrb4surf   ([0 0], [1 0], [0 1], [1 1]);
 124 %! nrb = nrbdegelev (nrb, [p-1, q-1]);
 125 %! aux1 = linspace(0,1,m); aux2 = linspace(0,1,n);
 126 %! nrb = nrbkntins  (nrb, {aux1(2:end-1), aux2(2:end-1)});
 127 %! u = rand (1, 30); v = rand (1, 10);
 128 %! u = u - min (u); u = u / max (u);
 129 %! v = v - min (v); v = v / max (v);
 130 %! [B, N] = nrbbasisfun ({u, v}, nrb);
 131 %! assert (sum(B, 2), ones(300, 1), 1e-6)
 132 %! assert (all (all (B<=1)), true)
 133 %! assert (all (all (B>=0)), true)
 134 %! assert (all (all (N>0)), true)
 135 %! assert (all (all (N <= prod (nrb.number))), true)
 136 %! assert (max (max (N)),prod (nrb.number))
 137 %! assert (min (min (N)),1)