merge polytopes to decrease total number THIS FUNCTION IS NOT PUBLISHED
0001 function vrepList = mergeVrepList(vrepList, zerotol) 0002 % merge polytopes to decrease total number 0003 % 0004 % THIS FUNCTION IS NOT PUBLISHED 0005 0006 % Ducmentation draft for publishing: 0007 % 0008 % THIS FUNCTION IS PUBLISHED BUT TESTED AND MAINTAINED ONLY IN SCOPE 0009 % OF APPLICATIONS (decomposition case study; technical reports) 0010 % 0011 % Syntax: vrepList = mergeVrepList(vrepList) 0012 % 0013 % Reduce the set of simplices in vrepList by subsequent merging. The 0014 % algorithm that is used here is only implemented to generate some 0015 % figures in technical reports. 0016 % 0017 % See also: triangulateVrep, identifyPolytopeList, subtractPolytope 0018 0019 % The elk-library: convex geometry applied to crystallization modeling. 0020 % Copyright (C) 2012 Alexander Reinhold 0021 % 0022 % This program is free software: you can redistribute it and/or modify it 0023 % under the terms of the GNU General Public License as published by the 0024 % Free Software Foundation, either version 3 of the License, or (at your 0025 % option) any later version. 0026 % 0027 % This program is distributed in the hope that it will be useful, but 0028 % WITHOUT ANY WARRANTY; without even the implied warranty of 0029 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0030 % General Public License for more details. 0031 % 0032 % You should have received a copy of the GNU General Public License along 0033 % with this program. If not, see <http://www.gnu.org/licenses/>. 0034 0035 %% input 0036 if ~exist('zerotol', 'var') 0037 zerotol = elkZerotol; 0038 end 0039 0040 % The idea is the following: 0041 % 1) cycle through all polytopes until the last is reached 0042 % 2) merge this polytope with as much polytopes as possible: 0043 % 2a) select a facet (or the next facet) 0044 % 2b) find a neighbouring polytope that shares this facet 0045 % 2c) if the convex hull contains vertices of other polytope, include 0046 % them in the current list 0047 % 2d) if the volume of the convex hull of the current list equals the 0048 % sum of the volumes of the polytopes in the list, merge polytope: 0049 % * substitute current polytope with merged polytope 0050 % * delete all merged polytopes in complete list 0051 % 2e) update information and go to 2a 0052 % 2f) if no merge was found, go to 2a anyway (without updating) 0053 % 2g) if no facet is left, loop finishes 0054 % 3) done 0055 0056 % counters and loop control 0057 iPoly = 1; 0058 iFacet = nan; 0059 filterNeighbour = 0; 0060 updateThisPoly = 1; 0061 while iPoly < length(vrepList) 0062 disp(['mergeVrepList: nPart = ' num2str(length(vrepList))]); 0063 % catch change of current polytope 0064 if updateThisPoly 0065 % polytope was merged, or the next polytope in list was chosen 0066 vrepMain = vrepList{iPoly}; 0067 hrepMain = convertVrepToHrep(vrepMain); 0068 updateThisPoly = 0; 0069 iFacet = 0; 0070 filterNeighbour = 0; 0071 continue; 0072 end 0073 0074 % catch change of facet, when neighbourList is empty 0075 if sum(filterNeighbour) <= 0 0076 % no merging happened and no neighbours on facet remain 0077 iFacet = iFacet + 1; 0078 % if all facets are searched, chose next polytope 0079 if iFacet > length(hrepMain.h) 0080 iPoly = iPoly + 1; 0081 updateThisPoly = 1; 0082 continue 0083 end 0084 0085 % determine new vrepNeighbourList 0086 % a neighbour must: 0087 % * must have n extreme points that can be expressed as a linear 0088 % combination of the points of the current facet 0089 facetVertexMatrix = vrepMain.V(... 0090 isZero(hrepMain.A(iFacet, :) * vrepMain.V' - hrepMain.h(iFacet), zerotol), ... 0091 :); 0092 % singular value decomposition for inner-loop calculations 0093 [u,s,v] = svd(facetVertexMatrix, 0); 0094 % uSum = sum(u, 1); 0095 % loop init 0096 filterNeighbour = ones(1, length(vrepList)) == 1; 0097 filterNeighbour(iPoly) = 0; 0098 for iNeighbour = find(filterNeighbour) 0099 thisVrepNeighbour = vrepList{iNeighbour}; 0100 0101 thisFilterVector = isZero(... 0102 hrepMain.A(iFacet, :) * thisVrepNeighbour.V' - hrepMain.h(iFacet), ... 0103 zerotol); 0104 0105 % quick check, that at least n points of the neighbour must 0106 % reside on the supporting hyperplane of the facet 0107 if sum(thisFilterVector) < size(vrepMain.V, 2) 0108 filterNeighbour(iNeighbour) = 0; 0109 % disp('first drop out') 0110 continue 0111 end 0112 % disp('accepted') 0113 % hrepMain.A(iFacet, :) * thisVrepNeighbour.V' - hrepMain.h(iFacet) 0114 0115 % for the vertices of the neighbour that are on the facet of 0116 % the current polytope, it must hold: 0117 % lam * vFacet = v; lam > 0 0118 % where lam is [lam1, .., lamN], vFacet are all facets of the 0119 % current polytope and v (column vector) is a single vertex 0120 % that already resides on the facet. 0121 % accorindg to singular value decomposition, we might have: 0122 % vFacet = U*S*(V') 0123 % for which we introduce: 0124 % gam = lam*U OR lam = gam*(U') 0125 % to obtain: 0126 % gam * S = v*V; gam*(U') > 0 0127 % or finally: 0128 % v*V*(S^-1)*(U') > 1 0129 nVertexOnFacet = 0; 0130 % for the vertices that are already on the supporting 0131 % hyperplane, count the ones that are on the facet 0132 for iVertex = find(thisFilterVector) 0133 thisVertex = thisVrepNeighbour.V(iVertex, :); 0134 if all(thisVertex*v*(diag(1./diag(s)))*(u') > -zerotol) 0135 nVertexOnFacet = nVertexOnFacet + 1; 0136 end 0137 end 0138 0139 if nVertexOnFacet >= size(vrepMain.V, 2) 0140 % disp('taken') 0141 filterNeighbour(iNeighbour) = 1; 0142 else 0143 % disp('second drop out') 0144 filterNeighbour(iNeighbour) = 0; 0145 end 0146 0147 end 0148 0149 % displayNeighbourFilter = filterNeighbour 0150 0151 % enter next part while reentering the loop 0152 continue; 0153 end 0154 0155 % check merging last entry in neighbour list 0156 if sum(filterNeighbour) > 0 && iFacet <= length(hrepMain.h) 0157 % skip that one for now: 0158 % filterNeighbour = 0; 0159 % continue 0160 0161 iNeighbour = find(filterNeighbour, 1, 'last'); 0162 0163 % add single polytope 0164 vrepMerged.V = [vrepMain.V; vrepList{iNeighbour}.V]; 0165 hrepMerged = convertVrepToHrep(vrepMerged); 0166 0167 % check for other polytopes 0168 filterPolyAlso = 1:length(vrepList) ~= iNeighbour; 0169 filterPolyAlso(iPoly) = 0; 0170 for iPolyAlso = find(filterPolyAlso) 0171 if any(prod(double(... 0172 bsxfun(@minus, ... 0173 hrepMerged.A * vrepList{iPolyAlso}.V', ... 0174 hrepMerged.h) < -zerotol*10), 1)) || ... 0175 all(prod(double(... 0176 bsxfun(@minus, ... 0177 hrepMerged.A * vrepList{iPolyAlso}.V', ... 0178 hrepMerged.h) < zerotol*10), 1)) 0179 % disp('added polytope') 0180 filterPolyAlso(iPolyAlso) = 1; 0181 vrepMerged.V = [vrepMerged.V; vrepList{iPolyAlso}.V]; 0182 else 0183 filterPolyAlso(iPolyAlso) = 0; 0184 end 0185 end 0186 0187 % the very final check for the volume 0188 totalVolume = 0; 0189 for iPart = [iPoly, iNeighbour, find(filterPolyAlso)] 0190 totalVolume = totalVolume + computeVolumeVrep(vrepList{iPart}); 0191 end 0192 0193 if abs(totalVolume - computeVolumeVrep(vrepMerged)) < zerotol 0194 % accept merging 0195 % listLengthBefore = length(vrepList) 0196 % mergedPolys = 2 + sum(filterPolyAlso) 0197 vrepList{iPoly} = reduceVrep(vrepMerged); 0198 filterDelete = filterPolyAlso; 0199 filterDelete(iNeighbour) = 1; 0200 vrepList(filterDelete) = []; 0201 % listLengthAfter = length(vrepList) 0202 updateThisPoly = 1; 0203 % disp('polytopes merged'); 0204 else 0205 if false 0206 totalVolume - computeVolumeVrep(vrepMerged) 0207 figure(66), hold on 0208 viewPolytope(vrepList{iPoly}, 'color', [1 0 0]) 0209 viewPolytope(vrepList{iNeighbour}, 'color', [0 1 0]) 0210 for iTmp = find(filterPolyAlso) 0211 % pause 0212 viewPolytope(vrepList{iTmp}) 0213 end 0214 % figure(67) 0215 % viewPolytope(vrepMerged) 0216 pause 0217 close(66) 0218 % close(67) 0219 end 0220 % neighbour did not succeed, delete it and try next one 0221 filterNeighbour(iNeighbour) = 0; 0222 % disp('polytopes rejected'); 0223 end 0224 else 0225 error('elk:polytope:numericalError', ... 0226 'Neighbour list is empty, this should not happen at this stage'); 0227 end 0228 0229 end