decompose the given polytope into Minkowski summands Syntax: hrepList = decomposeMinkowskiHrep(hrep) THIS FUNCTION IS NO USER FUNCTION.
0001 function summandMatrix = decomposeMinkowskiHrep(hrep, varargin) 0002 % decompose the given polytope into Minkowski summands 0003 % 0004 % Syntax: hrepList = decomposeMinkowskiHrep(hrep) 0005 % 0006 % THIS FUNCTION IS NO USER FUNCTION. 0007 0008 % ..event though, this function was originally designed to be. It might be 0009 % included in the polytope subsection in some later release. The main 0010 % missing feature is the generation of the summands itself. 0011 0012 % Draft documentation: 0013 % 0014 % Note that Minkowski decomposition is not unique. The result of this 0015 % algorithm gives indecomposable polytopes in hrepList so that it holds: 0016 % hrep = addMinkowskiPolytope(hrepList{1}, ..., hrepList{n}). 0017 % And because every hrepVector contains the original full matrix of facet 0018 % normals this is equivalent to: 0019 % hrep.h = hrepList{1}.h + .. + hrepList{n}.h 0020 % 0021 % The following options can be appended: 0022 % zerotol - tolerance for singularities and close to zero operations. 0023 % Default: 1e-13 0024 % verbose - switch on (1) and off (0) some verbosity. Default: 0 0025 % matrix - instead of a lost of H-representations a matrix with the 0026 % facet distances is returned (1). One row for each summand, 0027 % one column for each facet in hrep.A. Default: 0 0028 % result - This algorithm is also the main component for the 0029 % decomposition of the constrained h-space so that with the 0030 % option set to 'indecomposable', all indecomposable summands 0031 % are returned. Default: 'summand' 0032 % 0033 % See also: addMinkowskiPolytope, addMinkowskiHrep 0034 0035 % The elk-library: convex geometry applied to crystallization modeling. 0036 % Copyright (C) 2012 Alexander Reinhold 0037 % 0038 % This program is free software: you can redistribute it and/or modify it 0039 % under the terms of the GNU General Public License as published by the 0040 % Free Software Foundation, either version 3 of the License, or (at your 0041 % option) any later version. 0042 % 0043 % This program is distributed in the hope that it will be useful, but 0044 % WITHOUT ANY WARRANTY; without even the implied warranty of 0045 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0046 % General Public License for more details. 0047 % 0048 % You should have received a copy of the GNU General Public License along 0049 % with this program. If not, see <http://www.gnu.org/licenses/>. 0050 0051 %% sources 0052 % The original sketch (and proof) of the algorithm can be found in Meyer, 0053 % 1974: "Indecomposable Polytopes". It is modified and presented again in 0054 % Mount and Silverman, 1991: "Combinatorial and Computational Aspects of 0055 % Minkowski Decomposition". The algorithm here is adopted for constrained 0056 % Minkowski decomposition. It also has an easy to implement approach to 0057 % eliminate the degrees of freedom by translation. In Meyer 1974 the 0058 % Steiner-point was used. 0059 0060 %% Input 0061 % some numbers 0062 nDim = size(hrep.A, 2); 0063 nFacet = size(hrep.A, 1); 0064 0065 % map options 0066 [zerotol, verbose, outputResult, outputMatrix, constraintMatrix] = mapOption(varargin, ... 0067 'zerotol', elkZerotol, ... 0068 'verbose', 0, ... 0069 'result', 'summand',... 0070 'matrix', 0, ... 0071 'constraint', zeros(1, nFacet)); 0072 0073 % check polytope vs. constraint 0074 if ~isAllZero(constraintMatrix * hrep.h, zerotol) 0075 error('elk:decomposition:inputError', ... 0076 'The input polytope must fulfill the provided constraints'); 0077 end 0078 0079 %% Fix position & Scaling 0080 % get additional constraint and move hrep 0081 [positionConstraintMatrix, mapFacetDistanceToOffset] = ... 0082 fixPosition(hrep.A, constraintMatrix, zerotol); 0083 positionFixOffset = mapFacetDistanceToOffset * hrep.h; 0084 hrepPositionFixed = moveHrep(hrep, positionFixOffset); 0085 0086 % Scaling 0087 % normal scaling cannot be applied according to position fixing. The only 0088 % thing that can be done is scaling the h-vector. We simply scale the 0089 % largest facet distance to one. 0090 scaleFactor = 1/max(abs(hrepPositionFixed.h)); 0091 hrepScaled = hrepPositionFixed; 0092 hrepScaled.h = hrepScaled.h * scaleFactor; 0093 0094 % check final polytope vs. input and position constraint 0095 if ~isAllZero(constraintMatrix * hrepScaled.h, zerotol) 0096 error('elk:decomposition:wrongInput', ... 0097 'The input polytope must fulfill the given constraints'); 0098 end 0099 if ~isAllZero(positionConstraintMatrix * hrepScaled.h, zerotol) 0100 error('elk:decomposition:numericalProblem', ... 0101 'The scaled polytope does not fulfill the fixed position constraint. This should not happen!'); 0102 end 0103 0104 if verbose,disp('Generating incidence matrix');end 0105 0106 %% Generate incidence matrix 0107 % Facets, incident to a vertex are stored in an incidence matrix. Rows for 0108 % vertices, columns for facets. 0109 0110 % extreme points and initialization 0111 vrepScaled = convertHrepToVrep(hrepScaled, zerotol); 0112 nExtreme = size(vrepScaled.V, 1); 0113 incidenceVerticesToFacets = zeros(nExtreme, nFacet, 'uint8'); 0114 0115 for iExtreme = 1:nExtreme 0116 incidenceVerticesToFacets(iExtreme, :) = isZero( ... 0117 hrepScaled.A * vrepScaled.V(iExtreme, :)' - hrepScaled.h, ... 0118 zerotol); 0119 end 0120 0121 %% Equations and Inequalities 0122 % The following loop goes through all extreme point / facet 0123 % combinations. While it follows three states: 0124 % (1) finding the representative facets 0125 % (2) adding equality constraints 0126 % (3) adding inequality constraints 0127 % Both matrices are stored in the form that each row is a equation 0128 % (or inequality) and each column corresponds to a facet. 0129 0130 if verbose,disp('Extracting equations and inequalities');end 0131 0132 % initialize equality/inequality matrices 0133 nEquation = sum(sum(incidenceVerticesToFacets, 2) - nDim); 0134 nInequality = sum(sum(~incidenceVerticesToFacets, 2)); 0135 equationMatrix = zeros(nEquation, nFacet); 0136 inequalityMatrix = zeros(nInequality, nFacet); 0137 thisEquationIndex = 0; 0138 thisInequalityIndex = 0; 0139 0140 % loop over extreme points 0141 for iExtreme = 1:nExtreme 0142 0143 thisAdjacentFacetVector = find(incidenceVerticesToFacets(iExtreme, :)); 0144 thisNonadjacentFacetVector = find(~incidenceVerticesToFacets(iExtreme, :)); 0145 0146 % find representative facets 0147 thisRank = 0; 0148 repFacetMatrix = zeros(nDim, nDim); 0149 repFacetIndexVector = zeros(1, nDim); 0150 for iFacet = thisAdjacentFacetVector 0151 repFacetMatrix(thisRank+1, :) = ... 0152 hrepScaled.A(iFacet, :); 0153 if rank(repFacetMatrix) > thisRank 0154 thisRank = thisRank + 1; 0155 repFacetIndexVector(thisRank) = iFacet; 0156 if thisRank == nDim, break;end 0157 end 0158 end 0159 0160 % error if representatives not found 0161 if thisRank < nDim 0162 error('elk:decomposition:numericalProblem', ... 0163 'Failed to find representative facets. This should not happen!'); 0164 end 0165 0166 % the inversion of the representative facet normals 0167 % this is used together with the fourth facet normal (a4). For all 0168 % four one can find Li that satisy (L4 = 1): 0169 % A' * [L1, L2, L3]' = -L4*a4 0170 % by using: 0171 % [L1, L2, L3] = -L4 * a4' * inv(A) 0172 invertedRepresentativeNormalMatrix = inv(repFacetMatrix); 0173 0174 % eliminating representatives from adjacency vector 0175 thisAdjacentFacetVector = eliminateRows(thisAdjacentFacetVector', ... 0176 repFacetIndexVector', eps)'; 0177 0178 % adding equations 0179 % Knowing L1, L2, L3 and L4 we gather the equation: 0180 % [L1 L2 L3 L4] * [h1; h2; h3; h4] = 0 0181 for iFacet = thisAdjacentFacetVector 0182 thisEquationIndex = thisEquationIndex + 1; 0183 equationMatrix(thisEquationIndex, iFacet) = 1; 0184 lambdaVector = -1 * hrepScaled.A(iFacet, :) * ... 0185 invertedRepresentativeNormalMatrix; 0186 equationMatrix(thisEquationIndex, ... 0187 repFacetIndexVector) = lambdaVector; 0188 end 0189 0190 % adding inequalities 0191 % Knowing L1, L2, L3 and L4 we gather the equation: 0192 % -1 * [L1 L2 L3 L4] * [h1; h2; h3; h4] <= 0 0193 for iFacet = thisNonadjacentFacetVector 0194 thisInequalityIndex = thisInequalityIndex + 1; 0195 inequalityMatrix(thisInequalityIndex, iFacet) = -1; 0196 lambdaVector = 1 * hrepScaled.A(iFacet, :) * ... 0197 invertedRepresentativeNormalMatrix; 0198 inequalityMatrix(thisInequalityIndex, ... 0199 repFacetIndexVector) = lambdaVector; 0200 end 0201 0202 end 0203 0204 if thisEquationIndex < nEquation 0205 error('elk:decomposition:numericalProblem', ... 0206 'Failed to add all equations. This should not happen!'); 0207 end 0208 if thisInequalityIndex < nInequality 0209 error('elk:decomposition:numericalProblem', ... 0210 'Failed to add all inequalities. This should not happen!'); 0211 end 0212 0213 if verbose,disp(['Generated ' num2str(nEquation) ' equations and ' ... 0214 num2str(nInequality) ' inequalities']);end 0215 0216 if verbose,disp('Solving equations and inequalities');end 0217 0218 %% Solve equations and inequalities 0219 summandMatrix = computeExtremePointsOfCone(... 0220 [equationMatrix; constraintMatrix; positionConstraintMatrix], ... 0221 inequalityMatrix, zerotol); 0222 nSummand = size(summandMatrix, 1); 0223 0224 %% Output handling 0225 0226 % generate summands or not 0227 switch lower(outputResult) 0228 case 'indecomposable' 0229 %! ToDo: there should be some scaling. But this does not seem to be 0230 % trivial with fixed position and in hrep. 0231 % Update: Scaling is now possible in a reasonable manner, using 0232 % the mean width (implemented in scaleToReal) 0233 %/ 0234 case 'summand' 0235 error('this option is currently switched off'); 0236 % get scalar factors 0237 % For any number of summands we search a solution for: 0238 % summandMatrix' * lam = h with lam >= 0. 0239 % We simply use an optimizer with the function to be minimized: 0240 % (1:nSummands) * lam 0241 % cddlibOutput = callCddLib('solve_lp', struct(... 0242 % 'A', [summandMatrix'; eye(nSummand)], ... 0243 % 'B', [hrepScaled.h; zeros(nSummand, 1)], ... 0244 % 'lin', 1:nFacet, ... 0245 % 'obj', nSummand:-1:1 ... 0246 % )); 0247 % scalarMultiples = cddlibOutput.xopt; 0248 0249 %! ToDo: This part does currently not work. I could need to triangulate the 0250 % cone to find an appropriate simplex partition for which the 0251 % condition lam > 0 is automatically fulfilled. 0252 0253 % summandMatrixReduced = (projectionMatrix * summandMatrix')'; 0254 % facetDistanceVectorReduced = (projectionMatrix * hrepScaled.h)'; 0255 0256 % triangulationIndexMatrix = computeTriangulationOfCone(summandMatrixReduced); 0257 % for iSimplex = 1:size(triangulationIndexMatrix, 1) 0258 % thisRayMatrix = summandMatrixReduced(triangulationIndexMatrix(iSimplex, :), :); 0259 % thisMultipleVector = inv(thisRayMatrix') * facetDistanceVectorReduced'; 0260 % if any(thisMultipleVector < 0) 0261 % scalarMultipleVector = []; 0262 % continue 0263 % else 0264 % scalarMultipleVector = zeros(size(summandMatrixReduced, 1), 1); 0265 % scalarMultipleVector(triangulationIndexMatrix(iSimplex, :)) = ... 0266 % thisMultipleVector; 0267 % break; 0268 % end 0269 % end 0270 % 0271 % if isempty(scalarMultipleVector) 0272 % error('elk:decomposition:numericalProblem', ... 0273 % 'Could not determine scalar coefficients for summands. This should not happen.'); 0274 % end 0275 % 0276 % for iSummand = 1:nSummand 0277 % summandMatrix(iSummand, :) = summandMatrix(iSummand, :) * ... 0278 % scalarMultipleVector(iSummand); 0279 % end 0280 % 0281 % % throw away almost zero rows 0282 % summandMatrix(... 0283 % prod(... 0284 % double(isZero(summandMatrix, zerotol))... 0285 % , 2) > 0, :) = []; 0286 % nSummand = size(summandMatrix, 1); 0287 otherwise 0288 error('elk:decomposition:inputIssue', ... 0289 'The option ''result'' is either ''summand'' or ''indecomposable'''); 0290 end 0291 0292 % revert scaling 0293 summandMatrix = summandMatrix / scaleFactor; 0294 0295 % revert position fix 0296 % this should only be done for the first summand and only for result type 0297 % 'summand' 0298 if strcmpi(outputResult, 'summand') 0299 summandMatrix(1,:) = (summandMatrix(1,:)' - hrep.A * positionFixOffset)'; 0300 end 0301 0302 switch outputMatrix 0303 case 1 0304 % nothing to do. 0305 0306 case 0 0307 % convert to a cell array of h-representations 0308 hrepTemplate.A = hrep.A; 0309 hrepList = cell(nSummand, 1); 0310 for iSummand = 1:nSummand 0311 hrepList{iSummand} = hrepTemplate; 0312 hrepList{iSummand}.h = summandMatrix(iSummand, :)'; 0313 end 0314 summandMatrix = hrepList; 0315 0316 otherwise 0317 error('elk:decomposition:wrongInput', ... 0318 'The option ''matrix'' is either 1 or 0'); 0319 end 0320