function hrep = convertVrepToHrepByConvhulln(vrep, zerotol)
% convert the V-representation of a polytope into a H-representation
%
% Syntax: hrep = convertVrepToHrepByConvhulln(vrep, zerotol)
%
% The input should already be scaled by calling function.
%
% THIS FUNCTION IS NO USER FUNCTION.

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2012 Alexander Reinhold
%
% This program 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 (at your 
%   option) any later version.
%
% This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.

%% handle low-dim cases
nDim = size(vrep.V, 2);
if nDim < 1
    hrep.A = [];
    hrep.h = [];
    return
elseif nDim == 1
    hrep.A = [1; -1];
    hrep.h = [max(vrep.V(:)); -1*min(vrep.V(:))];
    return
end

%% call convhulln
patchIndexMatrix = convhulln(vrep.V, {'Qt', 'Qs'});
nPatch = size(patchIndexMatrix, 1);

%% identify facets
facetNormalMatrix = nan(nPatch, nDim);
facetDistanceVector = nan(nPatch, 1);
for iPatch = 1:nPatch
    % generate matrix of vertices
    thisVertexMatrix = vrep.V(patchIndexMatrix(iPatch, :), :);
    % move so that one vertex equals the origin
    thisVertexMatrix = thisVertexMatrix - ones(nDim, 1) * thisVertexMatrix(1,:);
    
    % an additional check of dimensionality could be implemented here it
    %   should hold:
    %     rank(thisVertexMatrix) == nDim - 1
    
    % the normal vector v will fulfill the following equation:
    %     thisV * v = 0
    % Hence, it equals the eigenvector that corresponds to the eigenvalue
    %   of zero. While this is true, there are frequently multiple
    %   eigenvalues and the precission of the corresponsing eigenvectors is
    %   extremely bad.
    % To overcome this accuracy problem, the singular value decomposition
    %   is applied.
    [~, thisS, thisV] = svd(thisVertexMatrix);
    thisS = diag(thisS);
    thisVecIndex = find(abs(thisS) < zerotol);
    if numel(thisVecIndex) ~= 1       
%         error('elk:polytope:numericalError', ...
%             'This patch is empty. This should not happen.');
        facetNormalMatrix(iPatch, :) = zeros(1, nDim);
        facetDistanceVector(iPatch) = 0;
        continue;
    else
        facetNormalMatrix(iPatch, :) = thisV(:, thisVecIndex)';
    end
    vertexDistanceVector = facetNormalMatrix(iPatch, :) * vrep.V';
    facetDistanceVector(iPatch) = ...
        vertexDistanceVector(patchIndexMatrix(iPatch, 1));
    
    % The normal vector should be an outer normal. This cannot be decided
    %   on the facet vertices itself. A vertex that belongs to the polytope
    %   aside the facet is required. This vertes must lie in the negative 
    %   direction of the normal vector.
    vertexDistanceVector = vertexDistanceVector - facetDistanceVector(iPatch);
    maxDistance = max(vertexDistanceVector);
    minDistance = min(vertexDistanceVector);
    if maxDistance > zerotol
        % change sign, if an additional vertex lies in the direction of the
        %   normal vector
        facetDistanceVector(iPatch) = -1 * facetDistanceVector(iPatch);
        facetNormalMatrix(iPatch, :) = -1 * facetNormalMatrix(iPatch, :);
        % in this case minDistance should be around zero
        if ~(abs(minDistance) < zerotol)
            error('elk:polytope:numericalError', ...
                'The polytope has vertices inside and ouside of an identified halfspace');
        end
    elseif ~(minDistance < -zerotol)
        error('elk:polytope:numericalError', ...
            'The polytope has vertices inside and ouside of an identified halfspace');
    end
end

%% make unique normal vectors
% This happens, because the facets are triangulated and the normal vectors
%   for these simplices are identified.
[~, idx] = eliminateRows(facetNormalMatrix, ...
    [zeros(1, nDim); ...
     facetNormalMatrix], zerotol);
idx(idx == 1) = [];
idx = idx - 1;
idx = unique(idx);

%% write output
hrep.A = facetNormalMatrix(idx, :);
hrep.h = facetDistanceVector(idx);