%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% x = evaluate_spline_prime(coefs,s)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This function finds the second derivative at each point s of the
% splines specified by coefs.  
%
% Input:
% coefs is a (L + 3) x nsplines matrix, defining the 1-D B-spline.
% s is a ns x 1 array, containing the points at which to evaluate
% the derivative.
%
% Global variables required:
% BS_sigma and G_sigma: precomputed matrices needed for *any* cubic
% B-spline with L curve pieces.  These are precomputed to speed up
% computation.
%
% Output:
% x is a ns x nsplines matrix, containing the 2nd derivative of each of
% the splines evaluated at each point.  
%
% No global variables set.
%
% This function goes through each curve piece of the splines, finds
% points that lie on the current curve piece, and calls
% evaluate_spline_double_prime_sigma to evaluate the curve 2nd
% derivative at the current set of points.  
%
% Functions called by this function:
% evaluate_spline_prime_sigma
function x = evaluate_spline_double_prime(coefs,s)

% number of curve pieces
L = size(coefs,2) - 3;
  
% s is cyclic between 3 and L + 3
s = mod(s-3,L)+3;

x = nan*ones(size(s));

% Go through each of the curve pieces
for sigma = 3:L+2,
  % Find points that are on the current curve piece
  [sinds,Ninds] = find((s >= sigma) & (s < sigma+1));
  if ~isempty(Ninds),
    inds = sub2ind([ns,N],sinds,Ninds);
    % Evaluate the tangent at the points on the current curve piece
    x(inds) = evaluate_spline_double_prime_sigma(coefs(Ninds,:),s(inds));
  end;
end;