function varargout = reduceVrepDimension(vrep, zerotol)
% project lower dimensional (flat) polytope into embedding space
%
% Syntax: vrep = reduceVrepDimension(vrep, zerotol)
%   or    [vrep, rot, x0] = reduceVrepDimension(vrep, zerotol)
% 
% If the polytope vrep has a lower dimension r<n than the space it is
%   embedded in, this fucntion moves and rotates the polytope such that the
%   values for the last (n-r) coordinate axes are zero for all points in
%   vrep.
%
% The zerotol determines how flat a polytope must be to be recognized.
%   Points can be approximately a distance of zerotol away from the flat
%   polytope. Default: see elkZerotol.
%
% With the second syntax, the rotation matrix rot and the offset x0 is
%   returned so that the original polytope can be obtained by:
%     removedCoordinates = zeros(size(vrep, 1), size(V, 1) - size(vrep, 2));
%     vrepOriginal = [vrep, removedCoordinates]*rot';
%     vrepOriginal = moveVrep(vrepOriginal, x0);
%
% See also: reduceVrep, scaleVrep

% 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
if ~exist('zerotol', 'var')
    zerotol = elkZerotol;
end

%% ensure zero mean of points
vrepDim = size(vrep.V, 2);
x0 = mean(vrep.V, 1);
vrepMoved = moveVrep(vrep, -1*x0);

%% singular value decomposition
[u, s, rot] = svd(vrepMoved.V, 0);
filter = abs(diag(s)) >= zerotol;

%% assign output
% and do nothing whenever there was no reduction
if sum(filter) == vrepDim
    x0 = zeros(1, vrepDim);
    rot = diag(ones(1, vrepDim));
    varargout{1} = vrep;
else
    varargout{1} = vrepMoved.V*rot;
    varargout{1} = struct('V', varargout{1}(:, filter));
end

if nargout == 3
    varargout{2} = rot;
    varargout{3} = x0;
end