convert the H-representation of a polytope into V-representation Syntax: vrep = convertHrepToVrep(hrep) or vrep = convertHrepToVrep(hrep, zerotol) or vrep = convertHrepToVrep(hrep, zerotol, infinityBoxSize) The parameter zerotol applies for scaleHrep only (default is elkZerotol). Consider that cddlib applies an additional zerotol that is stored in the function zerotolCddlib during initialization of this library (buildCddlib). The parameter infinityBoxSize can be provided and is used when the given H-representation is unbounded. In that case the V-representation of the original object intersected with an infinity box is calculated. The infinity box is a hypercube, centered at the origin and extended to [-1, 1] for each dimension. Default: 1e3. There are two underlying methods to perform this conversion. The direct approach uses the double-description method (cddlib). An alternative approach is to utilize duality and use convertVrepToHrep with the underlying triangulation (convhulln). The default is choosen by speed and reliability. The alternative is used if this approach fails. convert(.., prefereCddlib) gives the possibility to change the default standard algorithm to be used. Default: 0 See also: convertPolytope, convertVrepToHrep, scaleHrep
0001 function vrep = convertHrepToVrep(hrep, zerotol, infinityBoxSize, useCddlib) 0002 % convert the H-representation of a polytope into V-representation 0003 % 0004 % Syntax: vrep = convertHrepToVrep(hrep) 0005 % or vrep = convertHrepToVrep(hrep, zerotol) 0006 % or vrep = convertHrepToVrep(hrep, zerotol, infinityBoxSize) 0007 % 0008 % The parameter zerotol applies for scaleHrep only (default is elkZerotol). 0009 % 0010 % Consider that cddlib applies an additional zerotol that is stored in the 0011 % function zerotolCddlib during initialization of this library 0012 % (buildCddlib). 0013 % 0014 % The parameter infinityBoxSize can be provided and is used when the given 0015 % H-representation is unbounded. In that case the V-representation of the 0016 % original object intersected with an infinity box is calculated. The 0017 % infinity box is a hypercube, centered at the origin and extended to 0018 % [-1, 1] for each dimension. Default: 1e3. 0019 % 0020 % There are two underlying methods to perform this conversion. The direct 0021 % approach uses the double-description method (cddlib). An alternative 0022 % approach is to utilize duality and use convertVrepToHrep with the 0023 % underlying triangulation (convhulln). The default is choosen by speed 0024 % and reliability. The alternative is used if this approach fails. 0025 % 0026 % convert(.., prefereCddlib) gives the possibility to change the default 0027 % standard algorithm to be used. Default: 0 0028 % 0029 % See also: convertPolytope, convertVrepToHrep, scaleHrep 0030 0031 % The elk-library: convex geometry applied to crystallization modeling. 0032 % Copyright (C) 2012 Alexander Reinhold 0033 % 0034 % This program is free software: you can redistribute it and/or modify it 0035 % under the terms of the GNU General Public License as published by the 0036 % Free Software Foundation, either version 3 of the License, or (at your 0037 % option) any later version. 0038 % 0039 % This program is distributed in the hope that it will be useful, but 0040 % WITHOUT ANY WARRANTY; without even the implied warranty of 0041 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0042 % General Public License for more details. 0043 % 0044 % You should have received a copy of the GNU General Public License along 0045 % with this program. If not, see <http://www.gnu.org/licenses/>. 0046 0047 %% validate input 0048 if ~exist('zerotol', 'var') || isempty(zerotol) 0049 zerotol = elkZerotol; 0050 end 0051 if ~exist('infinityBoxSize', 'var') || isempty(infinityBoxSize) 0052 infinityBoxSize = 1e3; 0053 end 0054 if ~exist('useCddlib', 'var') || isempty(useCddlib) 0055 useCddlib = 1; 0056 end 0057 if size(hrep.A, 1) ~= size(hrep.h, 1) 0058 error('elk:polytope:wrongType', ... 0059 'given H-representation is invalid, check if facet normals are a column vector'); 0060 end 0061 nDim = size(hrep.A, 2); 0062 0063 %% Shortcut for simplices 0064 % might be possible 0065 % if size(hrep.A, 1) == (dim + 1) && rank(hrep.A - hrep.A(1, 2) 0066 % end 0067 0068 %% reduce dimension & scaling 0069 [hrepReduced reduceTrafo reduceOffset] = reduceHrepDimension(hrep, zerotol); 0070 if isempty(hrepReduced) 0071 % empty polytope 0072 vrep = []; 0073 return 0074 end 0075 nDimReduced = size(hrepReduced.A, 2); 0076 if nDimReduced == 0 0077 % 0-dimensional polytope (single point) 0078 vrep.V = reduceOffset'; 0079 return 0080 end 0081 % The remaining polytope is at least one-dimensional. 0082 0083 %% scaling 0084 [hrepScaled, scale, scaleOffset] = scaleHrep(hrepReduced, zerotol); 0085 if scale == -inf 0086 % empty polytope 0087 vrep.V = zeros(0, nDim); 0088 return 0089 end 0090 if scale == 0 0091 % flat polytope 0092 error('elk:polytope:numericalError', ['Polytope appears to be flat ' ... 0093 'after reduceHrepDimension is aplied']); 0094 end 0095 if scale == inf 0096 % unbounded polytope >> apply infinityBox, retry scaling and continue 0097 boundingBoxHrep.A = [eye(nDimReduced); -1*eye(nDimReduced)]; 0098 boundingBoxHrep.h = infinityBoxSize*ones(2*nDimReduced,1); 0099 hrepReduced.A = [hrepReduced.A; boundingBoxHrep.A]; 0100 hrepReduced.h = [hrepReduced.h; boundingBoxHrep.h]; 0101 % hrepReduced.h = [hrepReduced.h; boundingBoxHrep.h] 0102 [hrepScaled, scale, scaleOffset] = scaleHrep(hrepReduced, 0.5*zerotol); 0103 if scale == inf || scale == 0 || scale == -inf 0104 error('elk:polytopeOperation:numericalError', ['Oh, come on! ' ... 0105 'first assumed to be unbounded and now the polytope appears ' ... 0106 'empty flat or even still unbounded. I am giving up. Not your ' ... 0107 'fault, though. This should not happen.']); 0108 end 0109 end 0110 %! DEV: the code below was used beforehand to additionally identify 0111 % unbounded polytopes. This should not be necessary. 0112 % [u singularValueVector v] = svd(hrep.A); 0113 % if size(singularValueVector, 1) == 1 || size(singularValueVector, 2) == 1 0114 % singularValueVector = singularValueVector(1); 0115 % else 0116 % singularValueVector = diag(singularValueVector); 0117 % end 0118 % if (scale == inf) || (~all(abs(singularValueVector) > zerotol)) 0119 % % also unbounded case 0120 % end 0121 %/ 0122 0123 %% normal processing 0124 % outer conversion structure 0125 useAlternative = 0; 0126 while useCddlib > -1 0127 try 0128 if useCddlib 0129 %% direct approach (cddlib) 0130 vrepScaled = convertHrepToVrepByCddlib(hrepScaled, 1); 0131 else 0132 %% dual approach (convhulln) 0133 % map to dual polytope in H-representationcomputeDualPolytope(hrepScaled) 0134 dualVrep = computeDualPolytope(hrepScaled); 0135 0136 % convert to dual Vrep 0137 dualHrep = convertVrepToHrep(dualVrep, zerotol, inf, 0); 0138 0139 % convert back to normal Hrep 0140 vrepScaled = computeDualPolytope(dualHrep); 0141 end 0142 % assign to exit loop: 0143 useCddlib = -1; 0144 catch thisError 0145 if ~useAlternative 0146 % we do not rethrow the original error, as the cddlib error does 0147 % not contain any hint. 0148 warning('elk:polytope:conversionMethod', ... 0149 ['Using alternative approach to convert polytope ' ... 0150 '(cddlib returned an error).']); 0151 useAlternative = 1; 0152 useCddlib = ~useCddlib; 0153 % try again 0154 continue 0155 end 0156 % safe input polytope 0157 if exist('log_convertHrepToVrep.mat', 'file') 0158 load log_convertHrepToVrep.mat 0159 logHrep{end+1} = hrep; 0160 else 0161 logHrep{1} = hrep; 0162 end 0163 save log_convertHrepToVrep.mat logHrep 0164 warning('elk:polytope:method', ... 0165 ['Given hrep cannot be converted to vrep (usually cddlib). ' ... 0166 'Polytope stored in: ' cd filesep 'log_convertHrepToVrep.mat']); 0167 rethrow(thisError); 0168 end 0169 end 0170 0171 %% Create result 0172 % revert scaling 0173 vrep = stretchVrep(vrepScaled, scale, zeros(nDimReduced, 1)); 0174 vrep = moveVrep(vrep, scaleOffset); 0175 % revert dimension cutoff 0176 if nDimReduced < nDim 0177 vrep.V = [vrep.V zeros(size(vrep.V, 1), nDim - nDimReduced)]; 0178 vrep.V = vrep.V * reduceTrafo; 0179 end