function varargout = convhull(varargin)
%MATLAB Code Generation Library Function

%   Notes:
%   * Supports only double precision 2D inputs

%   Copyright 2019 The MathWorks, Inc.
%#codegen

narginchk(1,4);

opts = parseOptions(varargin{:});
[pointsOffHull, ptsAreCollinear] = qhullReject2D(opts.pts, 2, opts.simplify);
coder.internal.errorIf(ptsAreCollinear, 'MATLAB:mathcgeo_catalog:EmptyConvhull2DErrId');

ptsLikelyOnHull = opts.pts(~pointsOffHull,:);
npts = size(ptsLikelyOnHull, 1);
[orientedPts, orientPerm] = orient(ptsLikelyOnHull);
[~, cornerPos] = min(orientedPts(:,2));

traversalOrder = coder.internal.indexInt([cornerPos:npts 1:(cornerPos-1)]);
traversalOrder = coder.internal.indexInt(orientPerm(traversalOrder));

hull = grahamComputeHull(ptsLikelyOnHull, traversalOrder, true(0,1), opts.simplify);

ptsPosOnHull = find(~pointsOffHull);
hull = ptsPosOnHull(hull);

% starting from minPointId to conform with MATLAB built-in convhull
[~, minPtId] = min(hull);
if minPtId ~= 1
    hull = [hull(minPtId:end-1); hull(1:minPtId-1); hull(minPtId)];
end
varargout{1} = hull;

if nargout > 1
    varargout{2} = computeConvhullArea(opts.pts, hull);
end

end

%% convhull function argument parser
function opts = parseOptions(varargin)
opts = struct( ...
    'pts', zeros(0,2), ...
    'simplify', false);

coder.varsize('opts.pts');

if nargin == 1
    validateattributes(varargin{1}, {'double'}, {'2d', 'finite'}, 'convhull', 'x', 1);
    coder.internal.errorIf(size(varargin{1}, 2) ~= 2, 'Coder:toolbox:ConvhullUnsupportedDim');    
    opts.pts = varargin{1};
elseif nargin == 2
    coder.internal.errorIf(~isnumeric(varargin{1}) || ~isnumeric(varargin{2}), 'MATLAB:mathcgeo_catalog:BadArgsErrId');
    validateattributes(varargin{1}, {'double'}, {'vector', 'finite'}, 'convhull', 'x', 1);
    validateattributes(varargin{2}, {'double'}, {'vector', 'finite'}, 'convhull', 'y', 2);
    sizex = size(varargin{1});
    sizey = size(varargin{2});
    coder.internal.errorIf(numel(sizex) ~= 2 || numel(sizey) ~=2, 'MATLAB:mathcgeo_catalog:BadPointSpecErrId');
    coder.internal.errorIf(sizex(1) ~= sizey(1) || sizex(2) ~= sizey(2), 'MATLAB:mathcgeo_catalog:MixDimDataPtCoordsErrId');
    opts.pts = [varargin{1}(:) varargin{2}(:)];
elseif nargin == 3
    validateattributes(varargin{1}, {'double'}, {'2d', 'finite'}, 'convhull', 'x', 1);
    coder.internal.errorIf(size(varargin{1}, 2) ~= 2, 'Coder:toolbox:ConvhullUnsupportedDim');        
    opts.pts = varargin{1};
    coder.internal.errorIf(~ischar(varargin{2}), 'MATLAB:mathcgeo_catalog:BadArgsErrId');
    validatestring(varargin{2}, {'simplify'}, 'parseOptions');
    coder.internal.errorIf(~islogical(varargin{3}), 'MATLAB:mathcgeo_catalog:BadArgsErrId');
    opts.simplify = varargin{3};
elseif nargin == 4
    coder.internal.errorIf(~isnumeric(varargin{1}) || ~isnumeric(varargin{2}), 'MATLAB:mathcgeo_catalog:BadArgsErrId');
    validateattributes(varargin{1}, {'double'}, {'vector', 'finite'}, 'convhull', 'x', 1);
    validateattributes(varargin{2}, {'double'}, {'vector', 'finite'}, 'convhull', 'y', 2);    
    sizex = size(varargin{1});
    sizey = size(varargin{2});
    coder.internal.errorIf(numel(sizex) ~= 2 || numel(sizey) ~=2, 'MATLAB:mathcgeo_catalog:BadPointSpecErrId');
    coder.internal.errorIf(sizex(1) ~= sizey(1) || sizex(2) ~= sizey(2), 'MATLAB:mathcgeo_catalog:MixDimDataPtCoordsErrId');
    opts.pts = [varargin{1}(:) varargin{2}(:)];
    coder.internal.errorIf(~ischar(varargin{3}), 'MATLAB:mathcgeo_catalog:BadArgsErrId');
    validatestring(varargin{3}, {'simplify'}, 'parseOptions');
    coder.internal.errorIf(~islogical(varargin{4}), 'MATLAB:mathcgeo_catalog:BadArgsErrId');
    opts.simplify = varargin{4};
end

npts = size(opts.pts, 1);
coder.internal.errorIf(npts < 3, 'MATLAB:mathcgeo_catalog:NotEnoughPtsConvhullErrId');
end

%% Methods to orient the points
function [orientedPts, orientPerm] = orient(pts)
npts = size(pts, 1);
[~, corner_pt_pos] = min(pts(:,2));

corner_pt = pts(corner_pt_pos, :);

slopes = zeros(npts, 1);
parfor i = 1:npts
    relativePt = pts(i,:) - corner_pt(1,:);  %#ok<PFBNS>
    slopes(i) = relativePt(2) / relativePt(1);
    if isnan(slopes(i))
        slopes(i) = 0;
    end
end

[sortedSlopes, orientPerm] = sort(slopes);
lineStretches = equalityStretches(sortedSlopes);
numStretches = size(lineStretches, 1);

for i = 1:numStretches
    orientPerm = coder.internal.introsort(orientPerm, ...
        coder.internal.indexInt(lineStretches(i,1)), ...
        coder.internal.indexInt(lineStretches(i,2)), @(i,j)comparePoints(i,j,pts));
end

for i = 1:numStretches
    if lineStretches(i,1) > 1
        lineEndPt1 = orientPerm(lineStretches(i,1));
        lineEndPt2 = orientPerm(lineStretches(i,2));
        prevPt = orientPerm(lineStretches(i,1)-1);
        
        if ccw(pts(prevPt,:), pts(lineEndPt1,:), pts(lineEndPt2,:)) < 0
            orientPerm(lineStretches(i,1):lineStretches(i,2)) = ...
                orientPerm(lineStretches(i,2):-1:lineStretches(i,1));
        end
    end
    
end

orientedPts = pts(orientPerm, :);
end


function isccw = ccw(p1, p2, p3)
isccw = (p2(1) - p1(1))*(p3(2) - p1(2)) - (p2(2) - p1(2))*(p3(1) - p1(1));
end


% compare function for row sorting of points
function cmp = comparePoints(i,j,pts)
if pts(j,2) > pts(i,2)
    cmp = true;
elseif pts(j,2) < pts(i,2)
    cmp = false;
else
    cmp = (pts(j,1) > pts(i,1));
end
end

% helper function to return the intervals in a vector
% that have consecutive equal elements.
function stretches = equalityStretches(v)
npts = numel(v);
stretches = zeros(ceil(npts/2), 2);
numStretch = 0;
numSkip = 0;
for i = 1:npts-1
    if numSkip ~= 0 && i <= numSkip
        continue
    end
    
    if v(i) == v(i+1)
        stretchEnd = i+1;
        numStretch = numStretch + 1;
        stretches(numStretch,1) = i;
        while stretchEnd <= npts-1 && v(stretchEnd) == v(stretchEnd+1)
            stretchEnd = stretchEnd + 1;
        end
        
        if stretchEnd == npts
            stretches(numStretch,2) = npts;
            break
        else
            stretches(numStretch,2) = stretchEnd;
            numSkip = stretchEnd;
        end
    end
end

stretches = stretches(1:numStretch,:);
end


%% Method to compute 2D convhull area
function area = computeConvhullArea(pts, hull)
n = numel(hull) - 1;
area = pts(hull(1),1)*(pts(hull(2),2) - pts(hull(n),2));
for k = 2:n
    area = area + pts(hull(k),1)*(pts(hull(k + 1),2) - pts(hull(k - 1),2));
end
area = 0.5*abs(area);
end