triangulate a polytope in V-representattion into simplices THIS FUNCTION IS PUBLISHED BUT TESTED AND MAINTAINED ONLY IN SCOPE OF APPLICATIONS (decomposition case study; technical reports) Syntax: vrepList = triangulateVrep(vrep) Generates a cell array of simplices in V-representation that presents a triangulation of the input polytope. To obtain this, the polytope is scaled and projected into an (n+1)-dimensional space to use qhulln. triangulateVrep(.., zerotol) lets you specify the zerotol that is used to identify values that are close to zero by zero. Default: elkZerotol. triangulateVrep(.., zerotol, isCone) is a specific logical parameter that identifies the input as a convex cone that contains the origin as one point. All other points must be at a hyperplane. In this case (isCone=1), the polytope is not elevated to a higher dimension. The polytope is also not scaled. See also: identifyPolytopeList, subtractPolytope
0001 function [vrepList, volumeVector] = triangulateVrep(vrep, zerotol, isCone) 0002 % triangulate a polytope in V-representattion into simplices 0003 % 0004 % THIS FUNCTION IS PUBLISHED BUT TESTED AND MAINTAINED ONLY IN SCOPE 0005 % OF APPLICATIONS (decomposition case study; technical reports) 0006 % 0007 % Syntax: vrepList = triangulateVrep(vrep) 0008 % 0009 % Generates a cell array of simplices in V-representation that presents a 0010 % triangulation of the input polytope. To obtain this, the polytope is 0011 % scaled and projected into an (n+1)-dimensional space to use qhulln. 0012 % 0013 % triangulateVrep(.., zerotol) lets you specify the zerotol that is used to 0014 % identify values that are close to zero by zero. Default: elkZerotol. 0015 % 0016 % triangulateVrep(.., zerotol, isCone) is a specific logical parameter that 0017 % identifies the input as a convex cone that contains the origin as one 0018 % point. All other points must be at a hyperplane. In this case 0019 % (isCone=1), the polytope is not elevated to a higher dimension. The 0020 % polytope is also not scaled. 0021 % 0022 % See also: identifyPolytopeList, subtractPolytope 0023 0024 % The elk-library: convex geometry applied to crystallization modeling. 0025 % Copyright (C) 2012 Alexander Reinhold 0026 % 0027 % This program is free software: you can redistribute it and/or modify it 0028 % under the terms of the GNU General Public License as published by the 0029 % Free Software Foundation, either version 3 of the License, or (at your 0030 % option) any later version. 0031 % 0032 % This program is distributed in the hope that it will be useful, but 0033 % WITHOUT ANY WARRANTY; without even the implied warranty of 0034 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0035 % General Public License for more details. 0036 % 0037 % You should have received a copy of the GNU General Public License along 0038 % with this program. If not, see <http://www.gnu.org/licenses/>. 0039 0040 %% Input handling 0041 if ~exist('zerotol', 'var') || isempty(zerotol) 0042 zerotol = elkZerotol; 0043 end 0044 if ~exist('isCone', 'var') || isempty(isCone) 0045 isCone = 0; 0046 end 0047 0048 %% reduce dimension if necessary 0049 %! ToDo: the following does only make sense for cones, for polytopes one 0050 % might look to computeVolumeVrep where this is done properly. 0051 reduce_dim = 0; 0052 if rank(vrep.V, zerotol) ~= size(vrep.V, 2), reduce_dim = 1; end 0053 0054 % NOTE: DEGENERATED POLYTOPE 0055 % This part was once here but no function should call triangulateVrep with 0056 % a degtenerated polytope. This is problematic when the polytope is 0057 % accidently degenerate. In this case, the lines below assign the volume 0058 % calculated in the lower-dimensional space. These are unpredictable 0059 % results. Therefore, we trigger an error here to see which function 0060 % actually might input such polytopes to this function. 0061 if reduce_dim 0062 vrepList = {}; 0063 volumeVector = []; 0064 return 0065 % error('elk:triangulation:wrongInput', ['The input polytope is ' ... 0066 % 'degenerate it can only be triangulated, if it is provided ' ... 0067 % ' as in a projected version into the appropriate subspace']); 0068 end 0069 % old code: 0070 % if reduce_dim 0071 % vrep_old = vrep; 0072 % [vrep, rot, x0] = reduceVrepDimension(vrep, zerotol); 0073 % end 0074 %/ 0075 0076 nDim = size(vrep.V, 2); 0077 0078 %% Catch trivial examples 0079 if nDim <= 1 0080 vrepList = {vrep}; 0081 volumeVector = computeVolumeVrep(vrep); 0082 return 0083 end 0084 0085 %% Vrep preparation 0086 if isCone 0087 % find origin 0088 originIndex = find(all(abs(vrep.V) < zerotol, 2)); 0089 if numel(originIndex) ~= 1 || isempty(originIndex) 0090 error('elk:decomposition:numericalProblem', ... 0091 'The cone does contain no or more than one origin. This should not happen!'); 0092 end 0093 vrepCone = vrep; 0094 else 0095 % scale polytope 0096 [vrepScaled, ~, ~] = scaleVrep(vrep, zerotol); 0097 0098 % elevate to cone 0099 vrepCone.V = [ones([size(vrepScaled.V, 1), 1]), vrepScaled.V; ... 0100 zeros([1, 1+size(vrepScaled.V, 2)])... 0101 ]; 0102 end 0103 0104 %% Triangulate cone 0105 patchIndexMatrix = convhulln(vrepCone.V, {'Qt', 'Qs'}); 0106 0107 %% Throw away patches that contain the origin 0108 if isCone 0109 patchIndexMatrix(sum(patchIndexMatrix == originIndex, 2) > 0 , :) = []; 0110 else 0111 filter = sum(patchIndexMatrix == max(patchIndexMatrix(:)), 2) > 0; 0112 patchIndexMatrix(filter, :) = []; 0113 end 0114 0115 %% Generate vrep of simplices 0116 vrepList = cell([1, size(patchIndexMatrix, 1)]); 0117 for iSimplex = 1:length(vrepList) 0118 if isCone 0119 vrepList{iSimplex} = struct('V', ... 0120 [vrepCone.V(patchIndexMatrix(iSimplex, :), :); ... 0121 zeros([1 size(vrepCone.V, 2)])... 0122 ]); 0123 else 0124 vrepList{iSimplex} = struct('V', ... 0125 vrep.V(patchIndexMatrix(iSimplex, :), :)... 0126 ); 0127 end 0128 end 0129 0130 %% Generate simplex volumes 0131 volumeVector = zeros([1, length(vrepList)]); 0132 for iSimplex = 1:length(vrepList) 0133 volumeVector(iSimplex) = computeVolumeVrep(vrepList{iSimplex}, zerotol); 0134 end 0135 0136 %% Throw away empty polytopes 0137 % .. that happens for convhulln 0138 % .. I don't want to mess with the implementation for cones, but the 0139 % implementation for polytopes seems to require some zerotol, here 0140 if isCone 0141 vrepList(volumeVector == 0) = []; 0142 volumeVector(volumeVector == 0) = []; 0143 else 0144 vrepList(volumeVector < zerotol) = []; 0145 volumeVector(volumeVector < zerotol) = []; 0146 end 0147 0148 %% Put back coordinates that were deleted through dimension reduction 0149 % Following lines are commented out, see note "DEGENERATED POLYTOPE", 0150 % above. 0151 % if reduce_dim 0152 % disp('.. done (reduced dimension).') 0153 % for iSimplex = 1:length(vrepList) 0154 % removedCoordinates = zeros(size(vrepList{iSimplex}.V, 1), size(vrep_old.V, 2) - size(vrep.V, 2)); 0155 % vrepList{iSimplex}.V = [vrepList{iSimplex}.V, removedCoordinates]*rot'; 0156 % vrepList{iSimplex} = moveVrep(vrepList{iSimplex}, x0); 0157 % end 0158 % end 0159 %/ 0160 0161 %% Output handling 0162 varagout{1} = vrepList; 0163 if nargout >= 2 0164 varargout{2} = volumeVector; 0165 end