%------------------------------------------------
% 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 [xx, yy] = LatexFilterHull(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;
    
    xx = nan(size(x, 1), NumberOfPoints*2+2);
    yy = nan(size(y, 1), NumberOfPoints*2+2);

    xxx = nan(size(x, 1), NumberOfPoints*2);
    yyy = nan(size(y, 1), NumberOfPoints*2);

    % Search maximum and minimum in each interval of length IntLength
    a = 1;
    b = 2; % always minimum first! -> easier to reorder
    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);
        xxx(:, 2*k+a) = x(:, k*IntLength+mi_idx);
        xxx(:, 2*k+b) = x(:, k*IntLength+ma_idx);
        yyy(:, 2*k+a) = mi;
        yyy(:, 2*k+b) = 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);
        xxx(:, 2*(k+NumberOfPoints-NrLongInt)+a) = x(:, sIdx+k*IntLengthLong+mi_idx);
        xxx(:, 2*(k+NumberOfPoints-NrLongInt)+b) = x(:, sIdx+k*IntLengthLong+ma_idx);
        yyy(:, 2*(k+NumberOfPoints-NrLongInt)+a) = mi;
        yyy(:, 2*(k+NumberOfPoints-NrLongInt)+b) = ma;
    end
    
    % always use first and last value
    xx(:, 1) = x(:, 1);
    yy(:, 1) = y(:, 1);
    xx(:, NumberOfPoints+2) = x(:, end);
    yy(:, NumberOfPoints+2) = y(:, end);
    % re-order points for use as a hull in latex
    for k=0:NumberOfPoints-1
        % minima
        xx(:, k+1+1) = xxx(:, 2*k+1);
        yy(:, k+1+1) = yyy(:, 2*k+1);
        % maxima
        xx(:, 2*NumberOfPoints-k+2) = xxx(:, 2*k+2);
        yy(:, 2*NumberOfPoints-k+2) = yyy(:, 2*k+2);
    end
    
%     fprintf(' Output - %i .\n', length(xx));
% 
%     figure;
%     plot(x, y, 'r');
%     hold on;
%     plot(xx, yy);

end