function [vrepList, volumeVector] = triangulateVrep(vrep, zerotol, isCone)
% triangulate a polytope in V-representattion into simplices
%
% THIS FUNCTION IS PUBLISHED BUT TESTED AND MAINTAINED ONLY IN SCOPE 
%   OF APPLICATIONS (decomposition case study; technical reports)
%
% Syntax: vrepList = triangulateVrep(vrep)
%
% Generates a cell array of simplices in V-representation that presents a
%   triangulation of the input polytope. To obtain this, the polytope is
%   scaled and projected into an (n+1)-dimensional space to use qhulln.
%
% triangulateVrep(.., zerotol) lets you specify the zerotol that is used to 
%   identify values that are close to zero by zero. Default: elkZerotol.
%
% triangulateVrep(.., zerotol, isCone) is a specific logical parameter that
%   identifies the input as a convex cone that contains the origin as one
%   point. All other points must be at a hyperplane. In this case 
%   (isCone=1), the polytope is not elevated to a higher dimension. The
%   polytope is also not scaled.
%
% See also: identifyPolytopeList, subtractPolytope

% 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/>.

%% Input handling
if ~exist('zerotol', 'var') || isempty(zerotol)
    zerotol = elkZerotol;
end
if ~exist('isCone', 'var') || isempty(isCone)
    isCone = 0;
end

%% reduce dimension if necessary
%! ToDo: the following does only make sense for cones, for polytopes one
%    might look to computeVolumeVrep where this is done properly.
reduce_dim = 0;
if rank(vrep.V, zerotol) ~= size(vrep.V, 2), reduce_dim = 1; end

% NOTE: DEGENERATED POLYTOPE
% This part was once here but no function should call triangulateVrep with
%   a degtenerated polytope. This is problematic when the polytope is
%   accidently degenerate. In this case, the lines below assign the volume
%   calculated in the lower-dimensional space. These are unpredictable
%   results. Therefore, we trigger an error here to see which function
%   actually might input such polytopes to this function.
if reduce_dim
    vrepList = {};
    volumeVector = [];
    return
%     error('elk:triangulation:wrongInput', ['The input polytope is ' ...
%         'degenerate  it can only be triangulated, if it is provided ' ...
%         ' as in a projected version into the appropriate subspace']);
end
% old code:
% if reduce_dim
%     vrep_old = vrep;
%     [vrep, rot, x0] = reduceVrepDimension(vrep, zerotol);
% end
%/

nDim = size(vrep.V, 2);

%% Catch trivial examples
if nDim <= 1
    vrepList = {vrep};
    volumeVector = computeVolumeVrep(vrep);
    return
end

%% Vrep preparation
if isCone
    % find origin
    originIndex = find(all(abs(vrep.V) < zerotol, 2));
    if numel(originIndex) ~= 1 || isempty(originIndex)
        error('elk:decomposition:numericalProblem', ...
            'The cone does contain no or more than one origin. This should not happen!');
    end
    vrepCone = vrep;
else
    % scale polytope
    [vrepScaled, ~, ~] = scaleVrep(vrep, zerotol);

    % elevate to cone
    vrepCone.V = [ones([size(vrepScaled.V, 1), 1]), vrepScaled.V; ...
                  zeros([1, 1+size(vrepScaled.V, 2)])...
                 ];
end

%% Triangulate cone
patchIndexMatrix = convhulln(vrepCone.V, {'Qt', 'Qs'});

%% Throw away patches that contain the origin
if isCone
    patchIndexMatrix(sum(patchIndexMatrix == originIndex, 2) > 0 , :) = [];
else
    filter = sum(patchIndexMatrix == max(patchIndexMatrix(:)), 2) > 0;
    patchIndexMatrix(filter, :) = [];
end

%% Generate vrep of simplices
vrepList = cell([1, size(patchIndexMatrix, 1)]);
for iSimplex = 1:length(vrepList)
    if isCone   
        vrepList{iSimplex} = struct('V', ...
            [vrepCone.V(patchIndexMatrix(iSimplex, :), :); ...
            zeros([1 size(vrepCone.V, 2)])...
            ]);
    else
        vrepList{iSimplex} = struct('V', ...
            vrep.V(patchIndexMatrix(iSimplex, :), :)...
            );
    end
end

%% Generate simplex volumes
volumeVector = zeros([1, length(vrepList)]);
for iSimplex = 1:length(vrepList)
    volumeVector(iSimplex) = computeVolumeVrep(vrepList{iSimplex}, zerotol);
end

%% Throw away empty polytopes
% .. that happens for convhulln
% .. I don't want to mess with the implementation for cones, but the
%    implementation for polytopes seems to require some zerotol, here
if isCone
    vrepList(volumeVector == 0) = [];
    volumeVector(volumeVector == 0) = [];
else
    vrepList(volumeVector < zerotol) = [];
    volumeVector(volumeVector < zerotol) = [];
end

%% Put back coordinates that were deleted through dimension reduction
% Following lines are commented out, see note "DEGENERATED POLYTOPE",
%   above.
% if reduce_dim
%     disp('.. done (reduced dimension).')
%     for iSimplex = 1:length(vrepList)
%         removedCoordinates = zeros(size(vrepList{iSimplex}.V, 1), size(vrep_old.V, 2) - size(vrep.V, 2));
%         vrepList{iSimplex}.V = [vrepList{iSimplex}.V, removedCoordinates]*rot';
%         vrepList{iSimplex} = moveVrep(vrepList{iSimplex}, x0);
%     end
% end
%/

%% Output handling
varagout{1} = vrepList;
if nargout >= 2
    varargout{2} = volumeVector;
end