home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 5.ddi / SPLINES.DI$ / SPVAL.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.6 KB  |  62 lines
  1. function v = spval(sp,x)
  2. % SPVAL    Evaluate spline.
  3. %
  4. %          v = spval(sp,x)
  5. %
  6. %  returns the vector  sp(x) .
  7.  
  8. % C. de Boor / latest change: May 23, 1989
  9. % Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
  10.  
  11.  
  12. %  if necessary, sort  x 
  13. xs=x(:)';lx=length(xs);tosort=0;
  14. if (~isempty(find(diff(xs)<=0))),
  15.    tosort=1;[xs,indx]=sort(xs);
  16. end
  17.  
  18. %  take apart  spline
  19. [t,a,n,k,d]=spbrk(sp);
  20.  
  21. %  and augment the knot sequence so that first and last knot each have
  22. %  multiplicity  k .
  23.  
  24. index=find(diff(t)>0);addl=k-index(1);addr=index(length(index))-n;
  25. if (addl>0|addr>0),
  26.    t=[t(1)*ones(1,addl) t(:)' t(n+k)*ones(1,addr)];
  27.    a=[zeros(d,addl) a zeros(d,addr)];
  28.    n=n+addl+addr;
  29. end
  30.  
  31. % for each data point, compute its knot interval
  32. index=max([sorted(t(1:n),xs);k*ones(1,lx)]);
  33.  
  34. % set up and carry out in lockstep the first spline evaluation algorithm:
  35.  
  36. b=a(:,index);
  37. if (k>1)
  38.    ones(d,1)*index;dindex=ans(:);
  39.    tx=ones(d*lx,1)*[2-k:k-1]+dindex*ones(1,2*(k-1));tx(:)=t(tx);
  40.    ones(d,1)*xs;tx=tx-ans(:)*ones(1,2*(k-1));
  41.    a=a(:);d*ones(d,1)*index+[1-d:0]'*ones(1,lx);dindex(:)=ans(:);
  42.    b=d*ones(d*lx,1)*[1-k:0]+dindex*ones(1,k);b(:)=a(b);
  43.  
  44.    % the following loop is taken from  sprpp .
  45.  
  46.    for r=1:k-1;
  47.       for i=1:k-r;
  48.          b(:,i) =(tx(:,i+k-1).*b(:,i)-tx(:,i+r-1).*b(:,i+1))./...
  49.                   (tx(:,i+k-1)-tx(:,i+r-1));
  50.       end
  51.    end
  52.  
  53. end
  54.  
  55. v=zeros(d,lx);v(:)=b(:,1);
  56. if tosort>0,[junk,indx]=sort(indx); v=v(:,indx);end
  57.  
  58. % finally, zero out all values for points outside the knot interval
  59. index=find(x<t(1)|x>t(n+1));
  60. if (~isempty(index)),
  61.    v(:,index)=zeros(d,length(index)); end
  62.