function [allx, ally] = quickcriticalpoints(in)

npts = length(in);
allx = zeros(1,npts/2);
ally = allx;
wasSmaller = in(2) > in(1);
wasPositive = in(1) > 0;
constStart = 0;
if in(2) == in(1); constStart = 1; end
idx = 1;

for i=2:npts-1
    isSmaller = in(i+1) > in(i);
    isPositive = in(i) > 0;
    isConstant = in(i+1) == in(i);
    
    if wasPositive ~= isPositive
        % check that the previous point was an extremum
        %       if (cp->n)
        %         assert(cp->y[cp->n-1] != 0);
        allx(idx) = i-0.5;
        ally(idx) = 0;
        idx = idx+1;
        wasPositive = isPositive;
    end
    
    % find extrema
    if constStart==0 & ~isConstant
        if wasSmaller ~= isSmaller
            diff = in(i+1)-in(i-1);
            den = (2*in(i)-in(i-1)-in(i+1));
            allx(idx) = i;
            ally(idx) = in(i);
            idx = idx+1;
            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
            allx(idx) = 0.5*(constStart+i);
            ally(idx) = in(i);
            idx=idx+1;
            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

% check for zero crossing on last point
isPositive = in(end) > 0;
if wasPositive ~= isPositive
    allx(idx) = length(in) + in(end)/(in(end)-in(end-1));
    ally(idx) = 0;
    idx=idx+1;
end

% clear unused points
allx(idx:end) = [];
ally(idx:end) = [];

if nargout < 2
    allx(2,:)=ally(1,:);
end