%------------------------------------------------
% 2011/01/02 version 1.0
%
% [xxx, yyy] = LatexFilterHullCurve(x, y, NumberOfPoints)
%
%------------------------------------------------
% output:
%  xxx, yyy: filtered data
% input:
%  x, y: data to be filtered
%  tol: tolerance

% author: Thomas Koenig
% date: 2011/01/02
%
% Copyright 2010 Thomas König, Alexander Michel
%
% This file is part of NumericPlots.
%
% NumericPlots is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% any later version.
% 
% NumericPlots is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with NumericPlots.  If not, see <http://www.gnu.org/licenses/>.

function [xxx, yyy] = LatexFilterMinMax(x, y, NumberOfPoints)

%     fprintf('LatexFilter: Input - %i ...', size(x,2));

    if (NumberOfPoints - round(NumberOfPoints))~=0
        error('NumberOfPoints must be an integer.');
    end
    if NumberOfPoints<3
        error('NumberOfPoints<3 does not make much sense.');
    end
    if NumberOfPoints>size(x, 2)/2
        % This would mean, that the filtered data would have more points than
        % the input data. This filter is not meant to do that.
        error('NumberOfPoints>length(x)');
    end

    IntLength = floor(size(x, 2)/NumberOfPoints);
    % The last NrLongInt intervals are enlarged by one to get all points in
    % x.
    NrLongInt = size(x,2)-(IntLength*NumberOfPoints);
    IntLengthLong = IntLength+1;
    
    xxx = nan(size(x, 1), NumberOfPoints*2+2);
    yyy = nan(size(y, 1), NumberOfPoints*2+2);

    % always use first and last value
    xxx(:,1) = x(:,1);
    yyy(:,1) = y(:,1);
    xxx(:,end) = x(:,end);
    yyy(:,end) = y(:,end);
    % Search maximum and minimum in each interval of length IntLength
    for k=0:NumberOfPoints-NrLongInt-1
        [mi, mi_idx] = min(y(:,(k*IntLength+1):(k+1)*IntLength), [], 2);
        [ma, ma_idx] = max(y(:,(k*IntLength+1):(k+1)*IntLength), [], 2);
        if mi_idx<ma_idx
            a=1;
            b=2;
        else
            a=2;
            b=1;
        end
        xxx(:, 2*k+a+1) = x(:, k*IntLength+mi_idx);
        xxx(:, 2*k+b+1) = x(:, k*IntLength+ma_idx);
        yyy(:, 2*k+a+1) = mi;
        yyy(:, 2*k+b+1) = ma;
    end
    sIdx = (NumberOfPoints-NrLongInt)*IntLength;
    for k=0:NrLongInt-1
        [mi, mi_idx] = min(y(:,(sIdx+k*IntLengthLong+1):(sIdx+(k+1)*IntLengthLong)), [], 2);
        [ma, ma_idx] = max(y(:,(sIdx+k*IntLengthLong+1):(sIdx+(k+1)*IntLengthLong)), [], 2);
        if mi_idx<ma_idx
            a=1;
            b=2;
        else
            a=2;
            b=1;
        end
        xxx(:, 2*(k+NumberOfPoints-NrLongInt)+a+1) = x(:, sIdx+k*IntLengthLong+mi_idx);
        xxx(:, 2*(k+NumberOfPoints-NrLongInt)+b+1) = x(:, sIdx+k*IntLengthLong+ma_idx);
        yyy(:, 2*(k+NumberOfPoints-NrLongInt)+a+1) = mi;
        yyy(:, 2*(k+NumberOfPoints-NrLongInt)+b+1) = ma;
    end
    
%     fprintf(' Output - %i .\n', length(xxx));
% 
%     figure;
%     plot(x, y, 'r');
%     hold on;
%     plot(xxx, yyy);
    

end