%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% x = evaluate_spline_double_prime_sigma(coefs,s)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This function finds the 2nd derivative at each point s of the
% spline curves 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.
%
function x = evaluate_spline_double_prime_sigma(coefs,s)

% Can evaluate more than one s, but they must all be in the same
% span (i.e. sigma <= s(i) < sigma + 1 for all i).

% Precomputed matrices to speed up evaluation
global BS_sigma;
global G_sigma;

% make sure s is a row vector
s = reshape(s,length(s),1);

% find the span
sigma = unique(floor(s));

% make sure all of the s(i) are in the same span
if length(sigma) ~= 1,
  error('All s must be in the same span');
  return;
end;

% Find the relevant BS_sigma and G matrices
B_sigma = squeeze(BS_sigma(:,:,sigma+1));
G = squeeze(G_sigma(:,:,sigma+1));

% Evaluate the curve derivative
one = ones(size(s));
call = (B_sigma * G * coefs')';

x = sum([0*one,0*one,2*one,6*s] .* call,2);