identify the edges of a V-representation
Syntax: vertexIndexMatrix = computeEdgeSetVrep(vrep)
outStruct = computeEdgeSetVrep(vrep, 'output', 'full')
The matrix vertexIndexMatrix gives in each row the indices into vrep.V
for the two vertices that define an edge. This algorithm is valid for
arbitrary dimensions.
computeEdgeSetVrep(.., 'output', out) lets you specify the type of
output that is provided. With out='normal', only the index matrix to
the vertices in vrep.V is returned. For out='full' a structure is
returned with the fields:
vrep - the V-representation of the polytope
hrep - the H-representation of the polytope
nEdge - the total number of edges
vertexIndexMatrix - indices to the vertices in vrep.V for the edges
facetIndexMatrix - indices to the (n-1) facets that define the edge
incidenceMatrix - a matrix with a row for each extreme point and a
column for each facet that is assigned with 1 when the vertex
belongs to the facet and 0 otherwise.
computeEdgeSetVrep(.., 'hrep', hrep) hands over the precalculated
H-representation of the polytope so that it is not calculated twice if
it is already known.
computeEdgeSetVrep(.., 'zerotol', zerotol) allows to set the tolerance
for underlying algorithms (determining flat polytopes; default in
elkZerotol).
See also: obtainVrep, viewPolytope0001 function outVar = computeEdgeSetVrep(vrep, varargin) 0002 % identify the edges of a V-representation 0003 % 0004 % Syntax: vertexIndexMatrix = computeEdgeSetVrep(vrep) 0005 % outStruct = computeEdgeSetVrep(vrep, 'output', 'full') 0006 % 0007 % The matrix vertexIndexMatrix gives in each row the indices into vrep.V 0008 % for the two vertices that define an edge. This algorithm is valid for 0009 % arbitrary dimensions. 0010 % 0011 % computeEdgeSetVrep(.., 'output', out) lets you specify the type of 0012 % output that is provided. With out='normal', only the index matrix to 0013 % the vertices in vrep.V is returned. For out='full' a structure is 0014 % returned with the fields: 0015 % vrep - the V-representation of the polytope 0016 % hrep - the H-representation of the polytope 0017 % nEdge - the total number of edges 0018 % vertexIndexMatrix - indices to the vertices in vrep.V for the edges 0019 % facetIndexMatrix - indices to the (n-1) facets that define the edge 0020 % incidenceMatrix - a matrix with a row for each extreme point and a 0021 % column for each facet that is assigned with 1 when the vertex 0022 % belongs to the facet and 0 otherwise. 0023 % 0024 % computeEdgeSetVrep(.., 'hrep', hrep) hands over the precalculated 0025 % H-representation of the polytope so that it is not calculated twice if 0026 % it is already known. 0027 % 0028 % computeEdgeSetVrep(.., 'zerotol', zerotol) allows to set the tolerance 0029 % for underlying algorithms (determining flat polytopes; default in 0030 % elkZerotol). 0031 % 0032 % See also: obtainVrep, viewPolytope 0033 0034 % The elk-library: convex geometry applied to crystallization modeling. 0035 % Copyright (C) 2012 Alexander Reinhold 0036 % 0037 % This program is free software: you can redistribute it and/or modify it 0038 % under the terms of the GNU General Public License as published by the 0039 % Free Software Foundation, either version 3 of the License, or (at your 0040 % option) any later version. 0041 % 0042 % This program is distributed in the hope that it will be useful, but 0043 % WITHOUT ANY WARRANTY; without even the implied warranty of 0044 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0045 % General Public License for more details. 0046 % 0047 % You should have received a copy of the GNU General Public License along 0048 % with this program. If not, see <http://www.gnu.org/licenses/>. 0049 0050 %% Input handling 0051 [hrep, zerotol, output] = mapOption(varargin, ... 0052 'hrep', [], ... 0053 'zerotol', elkZerotol, ... 0054 'output', 'normal'); 0055 0056 if isempty(hrep) 0057 hrep = convertVrepToHrep(vrep); 0058 end 0059 0060 nVertex = size(vrep.V, 1); 0061 nFacet = size(hrep.A, 1); 0062 nDim = size(hrep.A, 2); 0063 0064 %% Incidence matrix 0065 % generate an incidence matrix - rows for vertices and columns for facets - 0066 % a 1 is assigned when the vertex belongs to the corresponding facet - 0 0067 % is assigned otherwise. 0068 incidenceMatrix = zeros(nVertex, nFacet, 'uint8'); 0069 for iVertex = 1:size(vrep.V, 1) 0070 for iFacet = 1:size(hrep.A, 1) 0071 incidenceMatrix(iVertex, iFacet) = abs(... 0072 hrep.A(iFacet, :)*vrep.V(iVertex,:)' - hrep.h(iFacet)) < zerotol; 0073 end 0074 end 0075 0076 %% Evaluation of incidence 0077 % search all vertex pairs for which the logical and-operation of the 0078 % corresponding rows leaves exactly two ones. Less than two ones implies 0079 % that these two vertices do not have two facets in common (necessary 0080 % condition to form an edge). If more that two facets are adjacent for 0081 % both vertices, then it does no matter which two to take. 0082 vertexIndexMatrix = []; 0083 facetIndexMatrix = []; 0084 for iVertexOne = 1:nVertex 0085 for iVertexTwo = (iVertexOne+1):nVertex 0086 thisFacetFilter = incidenceMatrix(iVertexOne, :) & ... 0087 incidenceMatrix(iVertexTwo, :); 0088 nCommonFacet = sum(thisFacetFilter); 0089 if nCommonFacet >= (nDim-1) 0090 vertexIndexMatrix(end+1, :) = [iVertexOne, iVertexTwo]; 0091 facetIndexMatrix(end+1, :) = find(thisFacetFilter, 2); 0092 end 0093 % if nCommonFacet > (nDim-1) 0094 % error('elk:polytope:numericalError', ... 0095 % 'You have two vertices that belong to the same (n-1) facets. This should not happen.'); 0096 % end 0097 end 0098 end 0099 0100 %% Output handling 0101 switch lower(output) 0102 case 'normal' 0103 outVar = vertexIndexMatrix; 0104 case 'full' 0105 outVar.hrep = hrep; 0106 outVar.vrep = vrep; 0107 outVar.vertexIndexMatrix = vertexIndexMatrix; 0108 outVar.facetIndexMatrix = facetIndexMatrix; 0109 outVar.incidenceMatrix = incidenceMatrix; 0110 outVar.nEdge = size(vertexIndexMatrix, 1); 0111 otherwise 0112 error('elk:polytope:inputError', ... 0113 'Please use ''normal'' or ''full'' for the output parameter.'); 0114 end