calculates the set difference for a triangulated V-representation and any H-representation THIS FUNCTION IS PUBLISHED BUT TESTED AND MAINTAINED ONLY IN SCOPE OF APPLICATIONS (decomposition case study; technical reports) Syntax: vrepList = subtractPolytope(polyMinuend, polySubtrahend) The input minuend can be a V-representation, H-representation or a cell array of these representations. While in the latter case it is expected that this list contains only simplices. The subtrahend should be a H-representation of V-representation only. subtractPolytope(.., zerotol) lets you specify the zerotol that is used to identify values that are close to zero by zero. Default: elkZerotol. subtractPolytope(.., 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. This parameter is passed to triangulateVrep. See also: triangulateVrep, identifyPolytope
0001 function [resultVrepList resultVolumeVector resultIndexVector] = subtractPolytope(... 0002 polyMinuend, polySubtrahend, varargin) 0003 % calculates the set difference for a triangulated V-representation and any 0004 % H-representation 0005 % 0006 % THIS FUNCTION IS PUBLISHED BUT TESTED AND MAINTAINED ONLY IN SCOPE 0007 % OF APPLICATIONS (decomposition case study; technical reports) 0008 % 0009 % Syntax: vrepList = subtractPolytope(polyMinuend, polySubtrahend) 0010 % 0011 % The input minuend can be a V-representation, H-representation or a cell 0012 % array of these representations. While in the latter case it is expected 0013 % that this list contains only simplices. The subtrahend should be a 0014 % H-representation of V-representation only. 0015 % 0016 % subtractPolytope(.., zerotol) lets you specify the zerotol that is used 0017 % to identify values that are close to zero by zero. Default: elkZerotol. 0018 % 0019 % subtractPolytope(.., zerotol, isCone) is a specific logical parameter that 0020 % identifies the input as a convex cone that contains the origin as one 0021 % point. All other points must be at a hyperplane. This parameter is 0022 % passed to triangulateVrep. 0023 % 0024 % See also: triangulateVrep, identifyPolytope 0025 0026 % The elk-library: convex geometry applied to crystallization modeling. 0027 % Copyright (C) 2012 Alexander Reinhold 0028 % 0029 % This program is free software: you can redistribute it and/or modify it 0030 % under the terms of the GNU General Public License as published by the 0031 % Free Software Foundation, either version 3 of the License, or (at your 0032 % option) any later version. 0033 % 0034 % This program is distributed in the hope that it will be useful, but 0035 % WITHOUT ANY WARRANTY; without even the implied warranty of 0036 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0037 % General Public License for more details. 0038 % 0039 % You should have received a copy of the GNU General Public License along 0040 % with this program. If not, see <http://www.gnu.org/licenses/>. 0041 0042 %% Input handling 0043 [zerotol, isCone, minuendVolumeVector, scalingVolume, infoMinuend] = ... 0044 mapOption(varargin, ... 0045 'zerotol', elkZerotol, ... 0046 'isCone', 0, ... 0047 'volumeVector', [], ... 0048 'scalingVolume', 1, ... 0049 'infoMinuend', []); 0050 0051 % minuend (ensure triangulated vrep) 0052 if isempty(infoMinuend) 0053 infoMinuend = identifyPolytope(polyMinuend); 0054 end 0055 if strcmpi(infoMinuend.rep, 'vrep') 0056 [minuendVrepList minuendVolumeVector] = triangulateVrep(... 0057 polyMinuend, zerotol, isCone); 0058 infoMinuend = identifyPolytope(minuendVrepList); 0059 elseif strcmpi(infoMinuend.rep, 'hrep') 0060 polyMinuend = convertHrepToVrep(polyMinuend); 0061 [minuendVrepList minuendVolumeVector] = triangulateVrep(... 0062 polyMinuend, zerotol, isCone); 0063 infoMinuend = identifyPolytope(minuendVrepList); 0064 elseif strcmpi(infoMinuend.rep, 'vrepList') 0065 if ~all(infoMinuend.m == infoMinuend.m(1)) 0066 error('elk:polytope:wrongInput', ... 0067 'The given list of polytopes in V-representation is not a triangulation'); 0068 end 0069 minuendVrepList = polyMinuend; 0070 elseif strcmpi(infoMinuend.rep, 'hrepList') 0071 if ~all(infoMinuend.m == infoMinuend.m(1)) 0072 error('elk:polytope:wrongInput', ... 0073 'The given list of polytopes in V-representation is not a triangulation'); 0074 end 0075 for iSimplex = 1:infoMinuend.nPoly 0076 polyMinuend{iSimplex} = convertHrepToVrep(polyMinuend{iSimplex}, zerotol); 0077 end 0078 minuendVrepList = polyMinuend; 0079 infoMinuend = identifyPolytope(minuendVrepList); 0080 else 0081 error('elk:polytope:wrongInput', ... 0082 'The input minuend should be vrep, hrep, vrepTri or hrepTri'); 0083 end 0084 0085 % subtrahend 0086 infoSubtrahend = identifyPolytope(polySubtrahend); 0087 if strcmpi(infoSubtrahend.rep, 'vrep') 0088 if ~isCone 0089 hrepSubtrahend = convertVrepToHrep(polySubtrahend, zerotol); 0090 else 0091 hrepSubtrahend = convertVrepToHrep(struct(... 0092 'V', [polySubtrahend.V; polySubtrahend.V(1,:)*0]), zerotol); 0093 closingFacetFilter = abs(hrepSubtrahend.h) > zerotol; 0094 hrepSubtrahend.A(closingFacetFilter, :) = []; 0095 hrepSubtrahend.h(closingFacetFilter) = []; 0096 end 0097 infoSubtrahend = identifyPolytope(hrepSubtrahend); 0098 elseif strcmpi(infoSubtrahend.rep, 'hrep') 0099 hrepSubtrahend = polySubtrahend; 0100 else 0101 error('elk:polytope:wrongInput', ... 0102 'The input subtrahend should be vrep or hrep'); 0103 end 0104 0105 %% Initial variables 0106 nDim = infoMinuend.dim; 0107 nMinuend = infoMinuend.nPoly; 0108 0109 resultVrepList = cell(0); 0110 resultIndexVector = []; 0111 resultVolumeVector = []; 0112 remainingVrepList = cell(0); 0113 remainingVolumeVector = []; 0114 remainingIndexVector = []; 0115 0116 0117 %% Algorithm overview 0118 0119 % This simple form of the set difference is based on the following term 0120 % that represents the set difference of vrepMinuend by hrepSubtrahend: 0121 % Union over all {x | x in vrepMinuend * x <= 0 AND 0122 % hrepSubtrahend(i, :) * x > = 0}. 0123 % The inner terms of the union are all convex, contain only points of the 0124 % minuend but no point of the subtrahend (except boundaries). However, 0125 % they might have full-dimensional intersections. This problem can be 0126 % resolved. 0127 % The outer algorithm works as follows: 0128 % 1) From the simplices in in vrepMinuend and the last facet 0129 % hrepSubtrahend(end,:), calculate a list of simplices that are 0130 % definitely contained in the difference and a remaining list. The 0131 % simplices in vrepResultList and vrepRemainingList have no full- 0132 % dimensional intersecntions and their union equals vrepMinuend. 0133 % 2) Delete the last facet from hrepSubtrahend. 0134 % 3) If hrepSubtrahend is empty, the algorithm terminates with 0135 % vrepResultList. 0136 % 4) resultList = {resultList thisFunction(vrepRemainingList, 0137 % hrepSubtrahend)} 0138 % 4) Back to step 1 with vrepMinuendList=remainingList and 0139 % hrepSubtrahend. 0140 0141 % For step one, each simplex in vrepMinuendList is handled separately. For 0142 % each vertex of such a simplex, there are three possibilities according 0143 % to the facet in hrepSubtrahend(end, :): 0144 % ( 0) - The vertex is on the facet 0145 % (-1) - The vertex is on the opposite site of the facet and, hence, 0146 % contained in the difference 0147 % ( 1) - The vertex is in the described halfspace and, hence, might be 0148 % contained in the difference, dependent on the other facets. 0149 % For each simplex, there are three resulting cases: 0150 % (0) All vertices are case-(0) vertices. This should not happen, as this 0151 % implies that the simplex is not full-dimensional 0152 % (1) All vertices are case-(-1) or case-(0) vertices. The simplex is 0153 % contained in the vrepResultList 0154 % (2) All vertices are case-(1) or case(0) vertices. The simplex might 0155 % be contained in the difference. They are stored in 0156 % vrepRemainingList. 0157 % (3) The vertices have mixed cases. This will be discussed in the 0158 % following. 0159 % For case (3) of a simplex, we need to determine the extreme points of the 0160 % intersection between the plane and the simplex. These extreme points 0161 % can be described by: 0162 % xNew = lam*x_(-1) + (1-lam)*x_(1) 0163 % for any combination of simplex vertices x_(-1) of case-(-1) and x_(1) 0164 % of case-(1). These points must fulfill: 0165 % hrepSubtrahend.A(end, :) * xNew = hrepSubtrahend.h(end) 0166 % From here, the unknown lam can be determined: 0167 % a*(lam*x_(-1) + (1-lam)*x_(1)) = h 0168 % lam*(a*x_(-1) - a*x_(1)) = h - a*x_(1) 0169 % lam = (h - a*x_(1)) / (a*x_(-1) - a*x_(1)). 0170 % All points xNew and the points x_(-1) are the V-representation of a 0171 % polytope that is definitely contained in the difference. All points 0172 % xNew and the points x_(1) are the V-representation of a polytope that 0173 % might be contained in the sum, dependent on the remaining facets. Both 0174 % polytopes are triangulated and the corresponding simplizes are added to 0175 % the appropriate lists. 0176 0177 % Along with the resultList, we store to which original simplex (from 0178 % minuendList) each simplex in resultList belongs. This makes it possible 0179 % to merge original simplices again that were splitted during the 0180 % algorithm. 0181 0182 0183 % %% Determine parts of the difference 0184 % % every combination of a hrep from hrepMinuend and a facet of 0185 % % hrepSubtrahend might create a new part. We initialize a cell array that 0186 % % contains the different hrep parts and a matrix that keeps the nonempty 0187 % % entries: 0188 % hrepPartList = cell(nMinuend, nSubtrahendFacet); 0189 % hrepPartNonemptyMatrix = zeros(nMinuend, nSubtrahendFacet, 'uint8'); 0190 0191 %% Treat last facet of hrepSubtrahend 0192 thisNormal = hrepSubtrahend.A(end,:); 0193 thisDistance = hrepSubtrahend.h(end); 0194 % start loop for minuends 0195 % disp(['subtractPolytope: nMinuend = ' num2str(nMinuend)]) 0196 for iMinuend = 1:nMinuend 0197 thisSimplex = minuendVrepList{iMinuend}; 0198 0199 %% Classify vertices according to last facet 0200 distanceVector = thisNormal * thisSimplex.V' - thisDistance; 0201 filterMinus = distanceVector < (-1 * zerotol); 0202 filterPlus = distanceVector > zerotol; 0203 filterZero = ~filterMinus & ~filterPlus; 0204 0205 %% Case (0) - Invalid 0206 if all(filterZero) 0207 error('elk:decomposition:numericalProblem', ... 0208 'A flat polytope (partition) was generated. This should not happen.'); 0209 end 0210 0211 %% Case (1) - Simplex fully contained 0212 if all(filterPlus | filterZero) 0213 resultVrepList{end+1} = thisSimplex; 0214 resultIndexVector(end+1) = iMinuend; 0215 if ~isempty(minuendVolumeVector) 0216 resultVolumeVector(end+1) = minuendVolumeVector(iMinuend); 0217 else 0218 resultVolumeVector(end+1) = computeVolumeVrep(thisSimplex, zerotol); 0219 end 0220 end 0221 0222 %% Case (2) - Simplex fully a "maybe" 0223 if all(filterMinus | filterZero) 0224 remainingVrepList{end+1} = thisSimplex; 0225 remainingIndexVector(end+1) = iMinuend; 0226 if ~isempty(minuendVolumeVector) 0227 remainingVolumeVector(end+1) = minuendVolumeVector(iMinuend); 0228 else 0229 remainingVolumeVector(end+1) = computeVolumeVrep(thisSimplex, zerotol); 0230 end 0231 end 0232 0233 %% Case (3) - Partially contained 0234 if any(filterPlus) && any(filterMinus) 0235 % prepare mid-point matrix 0236 midPointMatrix = nan([sum(filterMinus)*sum(filterPlus), nDim]); 0237 iMidPoint = 0; 0238 0239 % determine mid-points 0240 for iCaseMinus = find(filterMinus) 0241 for iCasePlus = find(filterPlus) 0242 iMidPoint = iMidPoint + 1; 0243 0244 % lam = (h - a*x_(1)) / (a*x_(-1) - a*x_(1)) 0245 lam = (thisDistance - thisNormal * thisSimplex.V(iCasePlus, :)') / ... 0246 (thisNormal * thisSimplex.V(iCaseMinus, :)' - ... 0247 thisNormal * thisSimplex.V(iCasePlus, :)'); 0248 0249 % xNew = lam*x_(-1) + (1-lam)*x_(1) 0250 midPointMatrix(iMidPoint, :) = ... 0251 lam * thisSimplex.V(iCaseMinus, :) + ... 0252 (1-lam) * thisSimplex.V(iCasePlus, :); 0253 end 0254 end 0255 0256 % generate result simplex/simplices 0257 [plusVrepList plusVolumeVector] = triangulateVrep(struct('V', ... 0258 [thisSimplex.V(filterPlus, :); ... 0259 thisSimplex.V(filterZero, :); ... 0260 midPointMatrix... 0261 ]), zerotol, isCone); 0262 % Note: DEGENERATE SPLITTING 0263 % When a boundary cuts the simplex only slighty because of 0264 % numerical issues, then the simplex is split into a part that 0265 % is roughly itself and a part that is flat/degenerate. Hence, 0266 % the following error is not reasonable. 0267 % if isempty(thisPartPlus) 0268 % error('elk:decomposition:numericalProblem', ... 0269 % 'Polytope in set difference not splitted properly. This should not happen.'); 0270 % end 0271 % But we still can ignore these "flat" polytopes: 0272 plusVrepList(plusVolumeVector < zerotol*scalingVolume) = []; 0273 plusVolumeVector(plusVolumeVector < zerotol*scalingVolume) = []; 0274 %/ 0275 nPlus = length(plusVrepList); 0276 [resultVrepList{end + (1:nPlus)}] = plusVrepList{:}; 0277 resultIndexVector(end + (1:nPlus)) = iMinuend; 0278 resultVolumeVector(end + (1:nPlus)) = plusVolumeVector(:); 0279 0280 % generate remaining simplex/simplices 0281 [minusVrepList minusVolumeVector] = triangulateVrep(struct('V', ... 0282 [thisSimplex.V(filterMinus, :); ... 0283 thisSimplex.V(filterZero, :); ... 0284 midPointMatrix... 0285 ]), zerotol, isCone); 0286 0287 % Commented out, see note "DEGENERATE SPLITTING", above. 0288 % if isempty(thisPartMinus) 0289 % error('elk:decomposition:numericalProblem', ... 0290 % 'Polytope in set difference not splitted properly. This should not happen.'); 0291 % end 0292 minusVrepList(minusVolumeVector < zerotol*scalingVolume) = []; 0293 minusVolumeVector(minusVolumeVector < zerotol*scalingVolume) = []; 0294 %/ 0295 nMinus = length(minusVrepList); 0296 [remainingVrepList{end + (1:nMinus)}] = minusVrepList{:}; 0297 remainingIndexVector(end + (1:nMinus)) = iMinuend; 0298 remainingVolumeVector(end + (1:nMinus)) = minusVolumeVector(:); 0299 end 0300 end 0301 0302 %% Recursively handle facet (end-1) from hrepSubtrahend 0303 if length(hrepSubtrahend.h) > 1 && ~isempty(remainingVrepList) 0304 % delete last facet (already treated) 0305 hrepSubtrahend.A(end, :) = []; 0306 hrepSubtrahend.h(end) = []; 0307 0308 % only the remaining part needs a decision whether they are contained 0309 % or not. This is a recursive call of this function. 0310 [vrepNextList nextVolumeVector nextIndexVector] = ... 0311 subtractPolytope(remainingVrepList, hrepSubtrahend, ... 0312 'zerotol', zerotol, 'isCone', isCone, ... 0313 'volumeVector', remainingVolumeVector, ... 0314 'scalingVolume', scalingVolume); 0315 0316 % assign result 0317 [resultVrepList{end + (1:length(vrepNextList))}] = vrepNextList{:}; 0318 resultIndexVector(end + (1:length(vrepNextList))) = ... 0319 remainingIndexVector(nextIndexVector); 0320 resultVolumeVector(end + (1:length(vrepNextList))) = nextVolumeVector(:); 0321 end 0322 0323 % %% recover uneccesary splitting of simplices 0324 % % the following is only performed, when the volume vector of the 0325 % % minuendList is provided 0326 0327 %! DEV: reprogram the following such that the lists/vectors are reduced in 0328 % only one step. This step seems to take far too long. Same for the lists 0329 % above. 0330 %/ 0331 if ~isempty(minuendVolumeVector) 0332 filterRemove = false(1, length(minuendVrepList)); 0333 for iMinuend = 1:length(minuendVrepList) 0334 thisFilter = (resultIndexVector == iMinuend); 0335 0336 if sum(thisFilter) > 1 && ... 0337 abs(sum(resultVolumeVector(thisFilter)) - ... 0338 minuendVolumeVector(iMinuend)) < zerotol*scalingVolume 0339 0340 % get real indices 0341 indexVector = find(thisFilter); 0342 0343 % we substitute the first element with the non-splitted original 0344 resultVrepList{indexVector(1)} = minuendVrepList{iMinuend}; 0345 resultVolumeVector(indexVector(1)) = minuendVolumeVector(iMinuend); 0346 % .. index to original remains unchanged 0347 0348 % remove all other simplices 0349 filterRemove(indexVector(2:end)) = 1; 0350 end 0351 end 0352 resultVrepList(filterRemove) = []; 0353 resultVolumeVector(filterRemove) = []; 0354 resultIndexVector(filterRemove) = []; 0355 end