splineMaximaMinima.m

function [maxima, minima] = splineMaximaMinima(ppSpline)
%Get a list of all the local maxima and minima in a piecewise spline fit generated by spline in a
%command such as ppSpline = spline(x,Y); ppSpline must be of order 4 (ppSpline.order == 4)
%Example:
% ppSpline = spline(1:100,rand(100,1));
% %To plot the results use the following:
% figure;
% plot(1:0.1:100,ppval(ppSpline,1:0.1:100));
% [maxima, minima] = splineMaximaMinima(ppSpline);
% hold on
% plot(maxima,ppval(ppSpline,maxima),'r*');
% plot(minima,ppval(ppSpline,minima),'g*');
%Author: M Arthington
%Date: 09/12/2009
%Updated: 03/07/2010 - Changed fnval to ppval to remove the example's dependency on the Spline
%toolbox. Thanks to Tom Clark for pointing this out

if ppSpline.order ~=4
	error('The input to splineMaximaMinima must be a 4th order pp-form spline (a cubic).')
end

%Extract spline parameters
[breaks,coefs]=unmkpp(ppSpline);

%Are there any maxima in a piece
discriminant = (4*coefs(:,2).^2-12.*coefs(:,1).*coefs(:,3));
logicalPiecesWithStationary = discriminant>=0;

%The spline is made up of pieces relative to the first break
t1 = breaks(logicalPiecesWithStationary)' + (-2*coefs(logicalPiecesWithStationary,2)+sqrt(discriminant(logicalPiecesWithStationary)))./(6*coefs(logicalPiecesWithStationary,1));
t2 = breaks(logicalPiecesWithStationary)' + (-2*coefs(logicalPiecesWithStationary,2)-sqrt(discriminant(logicalPiecesWithStationary)))./(6*coefs(logicalPiecesWithStationary,1));

%Find the locations of stationary points that lie within their breaks
t1IsStationary = t1<breaks([false; logicalPiecesWithStationary])' & t1>breaks(logicalPiecesWithStationary)';
t2IsStationary = t2<breaks([false; logicalPiecesWithStationary])' & t2>breaks(logicalPiecesWithStationary)';

%The +sqrt value is always the location of the minimum (for reference, the +sqrt part is not
%always the upper x value, if the cubed coefficient is negative it is the lower x value)
minima = t1(t1IsStationary);
maxima = t2(t2IsStationary);
	
end

% %The above method is fast, because it uses a vectorised formula.
% %The matlab method (using the Spline toolbox) of achieving this functionality is, AFAIK,  to do the following (but this is very slow):
% %Differentiate the spline
% derivativeSpline = fnder(ppSpline);
% %Find the maxima and minima everywhere
% maximaMinima = mean(fnzeros(derivativeSpline));
% %Find the second differential of the spline
% secondDerivativeSpline = fnder(derivativeSpline);
% %Determine whether maximaMinima values are at maxima or minima positions
% minima = maximaMinima(find(ppval(secondDerivativeSpline,maximaMinima)>0));
% maxima = maximaMinima(find(ppval(secondDerivativeSpline,maximaMinima)<0));