function [ListOfVertices ListOfEdges ListOfFaces] = ...
			CombinatorialStructureOfCyclicZonohedron(NumOfGens);
% Purpose		Find the combinatorial relationships between the
%				vertices, edges, and faces of the zonohedron generated by a set of
%				three-dimensional vectors in cyclic form, meaning that no vector is in the
%				convex cone of the remaining vectors.  The
%				generating vectors must be in general position: no set of 
%				three (or fewer) input vectors can be linearly dependent.
%
% 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.  The coefficients of the linear
%				combinations that make up the vertices are all either 0 or 1.  
%
%				While finding the vertices, edges, and faces of a zonohedron is a difficult
%				problem in general, that problem is much simpler for cyclic zonohedra.
%				Pp. 113-115 of [Centore2013] gives a general procedure.  Assume that the
%				generating vectors are numbered cyclically, either clockwise or
%				counterclockwise.  Then the vertices of the zonohedron can be divided into
%				levels, depending on the number of generating vectors that are summed up
%				to make a vertex.  Furthermore, each vertex is the sum of a consecutive
%				sequence of generators, allowing for wraparound.  The zeroth level is just
%				the origin.  The first level consists of the generators themselves, each
%				of which is a vertex.  The second level consists of sums of pairs of
%				consecutive vectors (generator 1 + generator 2, generator 2 + generator 3,
%				etc.), and so on for higher levels.  The highest, or nth level, is a single
%				vertex, consisting of the sum of all n generators.  The edges and faces of
%				a cyclic zonohedron can also be easily determined, following Eq. (8) in
%				[Centore2013].  This routine follows that formulation, stepping through
%				the levels, and constructing the vertices, edges, and faces that result
%				when each level is added.
%
%				In addition to its geometric structure, which depends on metric properties
%				like distances and angles, a polyhedron also has a combinatorial structure,
%				which consists of containment relationships between vertices, edges, and
%				faces.  The combinatorial structure of a cyclic zonohedron depends only on
%				the number of generators.  Each vertex is the sum of a consecutive sequence
%				of generators, and that sequence will always be a vertex, regardless of
%				the coordinates of the individual generators, as long as the generators have
%				been numbered cyclically.
%
%				This routine uses three output lists to express a cyclic zonohedron s
%				combinatorial structure.  The first list gives the vertices.  Each entry
%				in this list is a vector of indices, and the sum of the generators with
%				those indices is a vertex of the zonohedron.  The first entry in the list
%				is the empty vector, which corresponds to the origin, which is also a
%				vertex.  The second list gives the edges.  Each entry in this list consists
%				of a two-element vector.  The two elements are indices into the first list,
%				of vertices; the zonohedron has an edge between the two vertices with those
%				indices.  The third list gives the parallelogram faces.  Each entry is
%				a four-element vector, where the elements are indices in the list of
%				vertices.  If the fourth vertex is joined to the first, then those four 
%				vertices form a parallelogram face of the cyclic zonohedron.
%
%				[Centore2013] 	Paul Centore, A Zonohedral Approach to Optimal Colours,
%								Color Research and Application, Vol. 38, No. 2, pp. 110-119,
%								April 2013.
%
%				NumOfGens	A positive integer, giving the number of generating vectors, which
%							are assumed to be indexed in cyclic order.
%
%				ListOfVertices	A list of vectors.  The entries of each vector are indices
%								to the generators; summing up the generators with those
%								indices gives a vertex of the zonohedron.
%
%				ListOfEdges		A list of two-element vectors.  The entries of each vector 
%								are indices into the list of vertices.  If two vertices
%								appear in such a vector, then the zonohedron has an edge
%								joining them.
%
%				ListOfFaces		A list of four-element vectors.  The entries of each vector 
%								are indices into the list of vertices.  If four vertices
%								appear in such a vector, then they are the corners of a
%								parallelogram face on the zonohedron.
%
% Author		Paul Centore (January 28, 2016)
%
% Copyright 2016 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 <http://www.gnu.org/licenses/>.

% Initialize output variables
ListOfVertices = {}	;
ListOfEdges    = {}	;
ListOfFaces    = {}	;

% Handle the degenerate case of only one generator.  In that case, the
% zonohedron consists of one line segment, whose endpoints are the origin and that generator.
% Place that information in the output variables directly, and exit the routine.
if NumOfGens == 1
	ListOfVertices{1} = []		;
	ListOfVertices{2} = [1]		;
	
	ListOfEdges{1}    = [1 2]	;
	return
end

% Handle a second degenerate case, in which there are only two generators.  Then the
% zonohedron consists of one parallelogram, whose terminal point is the sum of the two
% generators.  Place that information in the output variables directly, and exit the routine.
if NumOfGens == 2
	ListOfVertices{1} = []		;
	ListOfVertices{2} = [1]		;
	ListOfVertices{3} = [2]		;
	ListOfVertices{4} = [1 2]	;
	
	ListOfEdges{1} = [1 2]	;
	ListOfEdges{2} = [1 3]	;
	ListOfEdges{3} = [2 4]	;
	ListOfEdges{4} = [3 4]	;
	
	ListOfFaces{1} = [1 2 4 3]	;
	return 
end

% Make combinatorial structure of cyclic zonohedron, assuming generators are given in
% cyclic order.  Step through the levels of the zonohedron, where each level consists
% of vertices with the same number of summands.  The first level and last level have
% different structures from the rest, so they are handled separately.  

% The first vertex is always the origin, which is the sum of no generators, so it is
% represented by the empty set.  
ListOfVertices{1} = []	;	

% At the first level, where the number of summands is 1, each generator appears as a 
% vertex, and is joined directly to the origin by an edge.  Record this information.
for FirstGenerator = 1:NumOfGens
	ListOfVertices{end + 1} = [FirstGenerator]		;
	ListOfEdges{end+1}      = [1, FirstGenerator+1]	;
end

for NumberOfSummands = 2:(NumOfGens-1)
	for FirstGenerator = 1:NumOfGens
		% In a zonohedron, each new vertex is constructed by adding a generator
		% to an already-constructed vertex (which is itself a sum of generators).  In a
		% cyclic zonohedron, with the generators numbered in cyclic order, each vertex is
		% a sum of consecutive generators (allowing for wraparound).  To make a new
		% vertex, then, take an already-constructed previous vertex, and add the next generator,
		% after the generators that are already listed for that vertex. 
		NextVertexIndex = length(ListOfVertices) + 1	; % Index for vertex being added
		PreviousVertexIndex = NextVertexIndex-NumOfGens				;
		PreviousVertex      = ListOfVertices{PreviousVertexIndex}	;	
		
		% An additional generator will be appended to the end of the vector of indices
		% for the previous vertex.  Since the indices are sequential, the index of the
		% new generator is just one more than the last entry of the previous vertex; we
		% must, however, allow for wraparound when going from the nth generator to the
		% first.
		if PreviousVertex(end) == NumOfGens - 1
			AdditionalGenerator = NumOfGens	;
		else
			AdditionalGenerator = mod(PreviousVertex(end)+1,NumOfGens)	;
		end
		% Make a new vertex by appending the next generator to the sequential list of
		% generators for the previous vertex
		ListOfVertices{NextVertexIndex} = [PreviousVertex, AdditionalGenerator]	;
		
		% Each time a new vertex is added, additional edges must be added that connect
		% the new vertex to already-constructed vertices.  Two new edges are needed for
		% each new vertex.
		ListOfEdges{end+1} = [PreviousVertexIndex, NextVertexIndex]	;
		% For the second edge, an adjustment for wraparound might be needed.
		if FirstGenerator == NumOfGens
			VertexAfterPreviousVertex = PreviousVertexIndex + 1 - NumOfGens	;
		else
			VertexAfterPreviousVertex = PreviousVertexIndex+1	;
		end
		ListOfEdges{end+1} = [VertexAfterPreviousVertex, NextVertexIndex]	;
		
		% Each time a new vertex is added, an additional face must be added that connects
		% the new vertex to already-constructed vertices.  That face always has the shape
		% of a parallelogram. Again, wraparound must be checked.
		if PreviousVertexIndex <= NumOfGens + 1 	
			% All vertices in first row connect to origin
			OriginatingVertexIndex = 1	;
		elseif PreviousVertexIndex == 1 + ((NumberOfSummands-1) * NumOfGens)
			OriginatingVertexIndex = PreviousVertexIndex - 2*NumOfGens + 1	;
		else
			OriginatingVertexIndex = PreviousVertexIndex - NumOfGens + 1	;
		end
		ListOfFaces{end+1} = [OriginatingVertexIndex, PreviousVertexIndex, ...
							  NextVertexIndex, VertexAfterPreviousVertex]	;
	end
end

% At the final, or nth level, there is one terminal vertex, which is the sum of all the 
% generators.  This vertex is joined by an edge to each vertex in the (n-1)th level, and
% there are additional parallelogram faces.
NumberOfVerticesSoFar = length(ListOfVertices) 	;

% Add the terminal vertex, which is the sum of all generators.
ListOfVertices{end+1} = [1:NumOfGens]			; 	

% There is an edge connecting the terminal vertex to each of the previous round of vertices
for ctr = 1:NumOfGens
	ListOfEdges{end+1} = [NumberOfVerticesSoFar - NumOfGens + ctr, length(ListOfVertices)]	;
end

% There are also parallelogram faces connecting the terminal vertex to previous vertices.
% Wraparound adjustments are sometimes required when adding these faces.
for ctr = 1:NumOfGens
	PreviousVertexIndex = NumberOfVerticesSoFar - NumOfGens + ctr	;
	if PreviousVertexIndex == NumberOfVerticesSoFar
		VertexAfterPreviousVertex = PreviousVertexIndex + 1 - NumOfGens	;
	else
		VertexAfterPreviousVertex = PreviousVertexIndex + 1	;
	end
	if PreviousVertexIndex == NumberOfVerticesSoFar
		OriginatingVertexIndex = PreviousVertexIndex - 2*NumOfGens + 1	;
	else
		OriginatingVertexIndex = PreviousVertexIndex - NumOfGens + 1	;
	end
	ListOfFaces{end+1} = [OriginatingVertexIndex, PreviousVertexIndex, ...
						  length(ListOfVertices), VertexAfterPreviousVertex]	;
end