home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 5.ddi / SPLINES.DI$ / BSPLINE.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.6 KB  |  56 lines
  1. function pp=bspline(t,window)
  2. % BSPLINE    Generate the B-spline for a given knot sequence.
  3. %
  4. %                 [pp=] bspline(t[,window])
  5. %
  6. %  plots the B-spline with knot sequence  t  as well as the polynomial
  7. %  pieces of which it is composed, in the window specified.
  8. %  Also returns its pp form if requested.
  9.  
  10. % C. de Boor / latest change: May 22, 1989
  11. % C. de Boor / latest change: 11 May 1992 (correct misprint)
  12. % Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
  13.  
  14. %  Generate the spline description
  15.  
  16. k=length(t)-1;
  17. if (k>1)
  18.    adds=ones(1,k-1);
  19.    tt=[adds*t(1) t(:)' adds*t(k+1)];j=k+1;
  20.    a=[adds*0 1 adds*0]; 
  21.  
  22.    %  From this, generate the pp description
  23.  
  24.    inter=find( diff(tt)>0 );l=length(inter);
  25.    tx=ones(l,1)*[2-k:k-1]+inter'*ones(1,2*(k-1));tx(:)=tt(tx);
  26.    tx=tx-tt(inter)'*ones(1,2*(k-1));
  27.    b=ones(l,1)*[1-k:0]+inter'*ones(1,k);b(:)=a(b);
  28.    c=sprpp(tx,b);x=[tt(inter) tt(2*k)];
  29. else
  30.   l=1;x=t;c=1;
  31. end
  32.  
  33. if (nargout>0), pp=ppmak(x,c,1); end
  34.  
  35. %  Now generate a mesh ...
  36.  
  37. step=100;xx=x(1)+[-10:step+10]*(x(l+1)-x(1))/step;nstep=length(xx);
  38.  
  39. %  and generate and plot the polynomial pieces
  40.  
  41. if (nargin>1), subplot(window); else, clf; end
  42. xxx=[xx(1) xx(nstep)]; yyy=[-1 2];axis([xxx yyy]);hold on;
  43. for j=1:k+1; plot([t(j) t(j)],yyy);end
  44. bspl=zeros(1,step+21);find(xx>=x(1));jh=ans(1);
  45. jsmax=3;js=jsmax;
  46. for j=1:l;js=js+1;if (js>jsmax), js=1;end
  47.    jl=jh; find(xx>=x(j+1));jh=ans(1);
  48.    pval=polyval(c(j,:),xx-x(j));bspl(jl:jh)=pval(jl:jh);
  49.    if (js==1),     plot(xx,pval,':');
  50.    elseif (js==2), plot(xx,pval,'-.');
  51.    else,           plot(xx,pval,'--');
  52.  
  53.    end
  54. end
  55. plot(xx,bspl);hold off;
  56.