function y = bsplineBasis(x, t, j, n)
% bsplineBasis: compute a bspline basis function.
%
%    y = bsplineBasis(x, t, j, n)
%
% Evaluates the j^th B-spline basis function of degree n for a given set of
% knots t at points x.
%
% Uses the Cox-de Boor recursion formula.
%
% Input Parameters:
%
%   x: Vector of points at which to evaluate the b-spline basis.
%
%   t: Vector of knots.  Normally the first and last will be repeated n
%   times.
%
%   j: Scalar specifying which basis function from the basis function
%   set.  0 <= j <= length(t) - n - 2.
%
%   n: Degree of basis functions.  Default n = 3 (cubics).
%
% Output Parameters:
%
%   y: Basis function evaluated at the points in x.
%
% Ian Mitchell for CPSC 303, 2/8/08.
  
  % Default parameter setting.
  if(nargin < 4)
    n = 3;
  end
  
  % There are m+1 knots.
  m = length(t) - 1;

  % Check validity of j.
  if((j < 0) || (j > m - n - 1))
    error([ 'Parameter j = ' num2str(j) ' not in range [ 0, ' ...
            num2str(m - n - 1) ' ]']);
  end
  
  % We actually construct the b-spline recursively.  Remember when
  % indexing into the knot vector that Matlab's indexing starts at 1.
  if(n == 0)
    % Base case.
    %if(j+2 == length(t))
    %  y = ((x >= t(j+1)) & (x <= t(j+2)));
    %else
      y = ((x >= t(j+1)) & (x < t(j+2)));
    %end
   
  else
    % If the two knots in the denominator are equal, we ignore the term.
    denom1 = t(j + 1 + n) - t(j + 1);
    if(denom1 == 0)
      y = zeros(size(x));
    else
      y = (x - t(j + 1)) .* bsplineBasis(x, t, j, n-1) / denom1;
    end
    
    denom2 = t(j + 2 + n) - t(j + 2);
    if(denom2 ~= 0)
      y = y + (t(j + 2 + n) - x) .* bsplineBasis(x, t, j+1, n-1) / denom2;
    end
  end