function vrepSum = addMinkowskiVrep(vrepOne , vrepTwo, zerotol, reduce)
% calculate the Minkowski sum of two polytopes in V-representation
%
% Syntax: vrepSum = addMinkowskiVrep(vrepOne, vrepTwo)
%   or    vrepSum = addMinkowskiVrep(vrepOne, vrepTwo, zerotol)
%   or    vrepSum = addMinkowskiVrep(vrepOne, vrepTwo, zerotol, reduce)
%
% The parameter zerotol applies for scaleVrep only, which is used for
%   the reduction to extreme points (see below). Default: see elkZerotol.
%
% Consider that cddlib applies an additional zerotol that is stored in the 
%   function zerotolCddlib during initialization of this library 
%   (buildCddlib).
% 
% The parameter reduce lets you decide whether the points that are
%   generated for the sum are reduced to the set of extreme points. Chosing
%   reduce=0 to switch this reduction off may increase performance,
%   dependent on subsequent operations, e.g. if you perform a polytope
%   conversion afterwards where this reduction happens automatically. 
%   Default: 1 (on).
%
% See also: convertPolytope, obtainPolytope, addMinkowskiPolytope

% 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')
    zerotol = elkZerotol;
end

if ~exist('reduce', 'var')
    reduce = 1;
end

if size(vrepOne.V, 2) ~= size(vrepTwo.V, 2);
    error('elk:polytope:wrongType', ...
        'both summands must have the same dimension');
end

%% Determine some sizes
mOne = size(vrepOne.V, 1);
mTwo = size(vrepTwo.V, 1);
dim = size(vrepOne.V, 2);
mResult = mOne*mTwo;

%% Construct a matrix filled repetedly with vrepOne
% This vector already has the size of the non-reduced result.
vrepSum.V = zeros(mOne * mTwo, dim);
for iRow = 1:mTwo
    vrepSum.V((1:mOne) + (iRow-1)*mOne,:) = vrepOne.V;
end

%% Add vrepTwo
for iRow = 1:mOne
    mThisRowVector = (1:mOne:mResult) + (iRow-1);
    vrepSum.V(mThisRowVector,:) = vrepSum.V(mThisRowVector,:) + vrepTwo.V;
end

%% Calculate convex hull and eliminate unused points
if reduce
    [vrepSum, scale, x0] = scaleVrep(vrepSum, zerotol);
    vrepSum = reduceVrep(vrepSum);
    vrepSum = stretchVrep(vrepSum, scale);
    vrepSum = moveVrep(vrepSum, x0);
end
