function varargout=intersub(x,xi,y)
% c=INTERSUB(x,xi)
% [yi,P]=INTERSUB(x,xi,y)
%
% Interpolation subdivision for (among other things)
% boundary implementation of fast wavelet transform
% by the lifting scheme.
%
% Can (only) use this to find prediction coefficients of
% increasingly higher-order INTERPOLATING wavelet
% transforms with Deslauriers-Dubuc SCALING FUNCTIONS.
%
% USE I: Give locations only
% For a set of points 'x' and a location 'xi'
% finds the length(x) filter coefficients
% necessary for the lifting wavelet transform
% with polynomial cancellation of order length(x)-1.
%
% USE II: Give locations and values
% Returns the interpolated values and the coefficients of
% the polynomial. Like POLYFIT/POLYVAL.
%
% Then, later on, if you want to compute yi given xi
% from x and y, all you do is c(:)'*y(:) and you've got it.
% This works even outside the domain of x, unlike INTERP1.
%
% If the samples spacings are integer, and the relative location
% of xi with respect to its neighbours is unchanged, you only need
% to calculate these coefficients once and the "stencil" can be
% used subsequently in all levels of successive refinement.
%
% Related terms: refinement, pyramid scheme, Neville's
% algorithm, Deslauriers-Lubuc scheme.
%
% Example I: (from Fernandez et al. (1996) Proc. SPIE 2825)
%
% x=[1 3 5 7];  % Figure 4 and Table 2
% intersub(x,8) % (exterior point)
% intersub(x,6) % 
% intersub(x,4) % (Case a; interior point)
% intersub(x,2) % (Case b; boundary point)
% intersub(x,0) % (exterior point)
%
% Example II:
%
% x=[1 3];
% intersub(x,2) % (Case a; regularly sampled)
% x=[1.4 2.8];
% intersub(x,2) % (Case a; regularly sampled)
%
% Example III:
%
%% The so-called 4-point scheme "stencil" is given by
% c=intersub([1 3 5 7],4)
%
% Example IV:
%
% x=sort(rand(4,1)*10);
% y=polyval([4 -3 1 0],x);
% xi=linspace(min(x),max(x),100);
% [yi,P]=intersub(x,xi,y);
% plot(x,y,'s-'); hold on
% plot(xi,yi,'b-')
% xii=mean(x);
% c=intersub(x,xii);
% hold on
% plot(xii,c(:)'*y(:),'*')
% hold off
%
% Example V:
%
% Note that by definition, the interpolated curve goes through the sample
% point: intersub([1 2 3 4],2) is what you expect in that case.
%
% Last modified by fjsimons-at-alum.mit.edu, 2/27/2003


if nargin==2
  % Only work with locations and calculate interpolating coefficients
  % Note that this is NOT Neville's algorithm... that'd work too, I guess.
  c=xi.^[length(x)-1:-1:0]*inv(vander(x(:)'))*eye(length(x));
  varargout{1}=c;
elseif nargin==3
  % This is really nothing but POLYFIT and POLYVAL
  % Work with locations and values, and calculate interpolated
  % values and polynomial coefficients. Except POLYFIT inverts the
  % Vandermonde matrix by least squares, in case it is tough, I guess.
  P=inv(vander(x(:)'))*y(:);
  varargout{2}=P(:)';
  if length(xi(:))==1
    yi=P(:)'*(xi.^[length(x)-1:-1:0])';
  else
    yi=polyval(P,xi);
  end
  varargout{1}=yi(:);
end
