function [mx_x, mx_y, mn_x, mn_y] = findextrema(in)

npts = length(in);
mx_x = zeros(1,npts/4);
mx_y = mx_x;
mn_x = mx_x;
mn_y = mx_x;
wasSmaller = in(2) > in(1);
constStart = 0;
if in(2) == in(1); constStart = 1; end
xidx = 1;
nidx = 1;

for i=2:npts-1
    isSmaller = in(i+1) > in(i);
    isConstant = in(i+1) == in(i);
    
    if constStart==0 & ~isConstant
        if wasSmaller ~= isSmaller
            diff = in(i+1)-in(i-1);
            den = (2*in(i)-in(i-1)-in(i+1));
            x = i + diff / (2*den);
            y = in(i) + (diff*diff) / (8*den);
            if wasSmaller
                mx_x(xidx) = x;
                mx_y(xidx) = y;
                xidx=xidx+1;
            else
                mn_x(nidx) = x;
                mn_y(nidx) = y;
                nidx=nidx+1;
            end
            wasSmaller = isSmaller;
        end 
    elseif constStart~=0 & isConstant
        continue; % Don't touch wasSmaller!
    elseif constStart~=0 & ~isConstant
        % We haven't touched wasSmaller since the constant zone started!
        if wasSmaller ~= isSmaller
            % put an "extremum" at the middle point
            x = 0.5*(constStart+i);
            y = in(i);
            if wasSmaller
                mx_x(xidx) = x;
                mx_y(xidx) = y;
                xidx=xidx+1;
            else
                mn_x(nidx) = x;
                mn_y(nidx) = y;
                nidx=nidx+1;
            end
            wasSmaller = isSmaller;
        end % otherwise it's just a spurious constant value along a sloped region.
        constStart = 0;
    elseif constStart==0 & isConstant
        constStart = i;
    end
end

% delete extra points
mx_x(xidx:end) = [];
mx_y(xidx:end) = [];
mn_x(nidx:end) = [];
mn_y(nidx:end) = [];

if nargout == 2
    mx_x(2,:) = mx_y(1,:);
    mx_y = [mn_x; mn_y]; % actually minima
end