function [AllParallelepipeds] = AllParallelepipedsInZonohedron(NumberOfGeneratingVectors)
% Purpose Find all the parallelepipeds in the zonohedron generated by an input set of
% three-dimensional vectors, all of whose entries are non-negative. The
% generating vectors must be in general position: no input
% vector can be identically 0, no two input vectors can be parallel, and no
% three input vectors can be linearly dependent. The parallelepipeds are
% those whose vertices are all nodes, that is, linear combinations of the
% generating vectors with all coefficients either 0 or 1.
%
% Description A zonohedron is the set of all linear combinations of a set of generating
% vectors, such that the coefficients of any linear combination is between
% 0 and 1, inclusive. A zonohedron is a convex polytope. Each edge of the
% zonohedron is a translation of one of the generating vectors. Each face
% of the zonohedron is a parallelogram. This structure makes it natural to
% dissect a zonohedron into parallelepipeds. To help in such dissections,
% this routine finds all the parallelepipeds whose vertices are nodes (a
% node is the sum of a discrete subset of the generating vectors; in other
% words, a node is a linear combination of generating vectors, each of whose
% coefficients is either a 0 or a 1).
%
% Each parallelepiped has an originating or offset vector. The originating
% vector is a node, so it is a sum of generating vectors. All the sides of
% a parallelepiped are translations of three generating vectors, which are
% different from the generating vectors that make up the originating vector.
% A point in the parallelepiped is a linear combination of the originating
% vector and the three side vectors, such that the combination coefficients
% for the side vectors are between 0 and 1 inclusive.
%
% The routine uses a simple exhaustive algorithm. First, all possible nodes
% are calculated (these correspond to the power set of the set of generating
% vectors). Each node is considered as a possible originating vector. To
% be an originating vector, there must be at least three other generating
% vectors, which are not used in making the originating vector. For each node,
% there is a set of remaining vectors, which are not used in making the node.
% If that node is an originating vector, then the three side vectors must be
% drawn from the set of remaining vectors. All combinations of three vectors,
% drawn from the set of remaining vectors, are found; each threesome gives a
% possible parallelepiped, when combined with the originating vectors.
%
% This routine finds all possible parallelepipeds, with no regard to whether
% the parallelepipeds dissect the zonohedron. In fact, it is expected that
% many parallelepipeds will overlap.
%
% NumberOfGeneratingVectors The number of generating vectors.
%
% AllParallelepipeds A matrix which has a column for each generating
% vector, and a row for each parallelepiped. All the entries
% are either 0, 1, or -1. Each row defines a parallelepiped.
% Three entries of each row are -1. These identify the three
% generating vectors that make up the sides of the parallelepiped.
% The remaining entries are either 0 or 1, and are the coefficients
% of the linear combination that makes up the originating vector.
%
% Author Paul Centore (March 30, 2014)
% Revised Paul Centore (May 3, 2014)
% ---Updated documentation.
%
% Copyright 2014 Paul Centore
%
% This file is part of MunsellAndKubelkaMunkToolbox.
%
% MunsellAndKubelkaMunkToolbox 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.
%
% MunsellAndKubelkaMunkToolbox 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 MunsellAndKubelkaMunkToolbox. If not, see .
% Find all possible nodes. Each subset in the power set of generating vectors gives a node.
CoefficientsOfNodes = powerset(1:NumberOfGeneratingVectors) ;
% Make a matrix, with one row for each node. Each matrix entry is either a 0 or 1. The
% node for each row is made by treating the row entries as coefficients in a linear
% combination of generating vectors. Equivalently, we could just sum up the generating
% vectors whose columns have a 1.
NumberOfNodes = size(CoefficientsOfNodes,2) ;
Nodes = [] ;
for ctr = 1:NumberOfNodes
LineForNodes = zeros(1,NumberOfGeneratingVectors) ;
LineForNodes(CoefficientsOfNodes{1,ctr}) = 1;
Nodes =[Nodes; LineForNodes] ;
end
% Consider each node as a possible originating vector for parallelepipeds. The side-generating
% vectors would then be threesomes of the vectors which do not contribute to the originating
% vector. Find all such threesomes that exist, for each originating vector.
AllParallelepipeds = [] ;
for ctr = 1:NumberOfNodes
ElementsInSubset = CoefficientsOfNodes{1,ctr} ;
RemainingElements = setdiff(1:NumberOfGeneratingVectors, ElementsInSubset) ;
NumberOfRemainingElements = length(RemainingElements) ;
% Find threesomes, if any exist, of the vectors that do not contribute to the node.
if NumberOfRemainingElements >= 3
Triplets = combinator(NumberOfRemainingElements,3,'c') ;
NumberOfTriplets = size(Triplets,1) ;
% For each threesome of side vectors, add a row to the list of all parallelepipeds.
% Place a -1 in the entries for the side vectors.
for idx = 1:NumberOfTriplets
LineInAllParallelepipeds = Nodes(ctr,:) ;
LineInAllParallelepipeds(RemainingElements(Triplets(idx,:))) = -1 ;
AllParallelepipeds = [AllParallelepipeds; LineInAllParallelepipeds] ;
end
end
end