appendUnifiedPartition

 add unified partition, calculate simplex partitions and perform book
   keeping. Up to here, I know the extreme rays of the unified
   partition.
 The strategy is:
   (1) Triangulate the vrep of the unified partition
   (2) Generate the simplex partitions (producing hreps from vreps is
       simple)
   (3) Filter the outer facets of the unified partition from the
       simplex partitions.

 THIS IS NO USER FUNCTION

• computeTriangulationOfCone compute the triangulation of a convex cone given by its extreme rays
• computeVolumeVrep compute the volume of a polytope in V-representation
• normalizeHrep normalize length of facet normals to 1

• addConfinementPartitionData add confined partitions
• addMinkowskiPartitionData add acucmulated data to decomposition structure
• computeExtendedEmbedding compute extended embedding from H- and original representation
• projectIncompleteDec project a complete proper representation

0001 function dec = appendUnifiedPartition(dec, seIndexVector, zerotol, ...
0002     isValid, isConfined)
0003 % add unified partition, calculate simplex partitions and perform book
0004 %   keeping. Up to here, I know the extreme rays of the unified
0005 %   partition.
0006 % The strategy is:
0007 %   (1) Triangulate the vrep of the unified partition
0008 %   (2) Generate the simplex partitions (producing hreps from vreps is
0009 %       simple)
0010 %   (3) Filter the outer facets of the unified partition from the
0011 %       simplex partitions.
0012 %
0013 % THIS IS NO USER FUNCTION
0014 
0015 % The elk-library: convex geometry applied to crystallization modeling.
0016 %   Copyright (C) 2012 Alexander Reinhold
0017 %
0018 % This program is free software: you can redistribute it and/or modify it
0019 %   under the terms of the GNU General Public License as published by the
0020 %   Free Software Foundation, either version 3 of the License, or (at your
0021 %   option) any later version.
0022 %
0023 % This program is distributed in the hope that it will be useful, but
0024 %   WITHOUT ANY WARRANTY; without even the implied warranty of
0025 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
0026 %   General Public License for more details.
0027 %
0028 % You should have received a copy of the GNU General Public License along
0029 %   with this program.  If not, see <http://www.gnu.org/licenses/>.
0030 
0031 %! Note: Why was the conversion vrep->hrep based on the conversion
0032 %           function?
0033 % Coming from the vrep, that could be given by:
0034 %     struct('V', [dec.seScaledToClose(seIndexVector, :); ...
0035 %                  zeros(1, dec.nC)]);
0036 %   is not good. Testing cs-decomposition throws errors for the
0037 %   conversion afterwards, which is:
0038 %     partition = convertVrepToHrep(vrep, zerotol);
0039 % It was implemented like above, but changed for the given reasons.
0040 %/
0041 
0042 %% initialization
0043 if ~isfield(dec, 'unifiedPartition') || ~exist('seIndexVector', 'var')
0044     dec.simplexPartition = cell(0);
0045     dec.unifiedPartition = cell(0);
0046     dec.simplexToUnifiedVector = [];
0047     dec.simplexIsValidVector = [];
0048     dec.simplexIsConfinedVector = [];
0049     dec.unifiedToSimplexVector = [];
0050     dec.unifiedIsValidVector = [];    
0051     dec.unifiedIsConfinedVector = [];
0052     dec.seIncidenceMatrix = eye(dec.nS, dec.nS);
0053     return
0054 end
0055 
0056 thisUnifiedIndex = length(dec.unifiedPartition) + 1;
0057 
0058 % triangulation of unified partition
0059 seSimplexIndexMatrix = computeTriangulationOfCone(...
0060     dec.seScaledToClose(seIndexVector, :));
0061 
0062 % prepare matrix of facet normals from all simplex partitions
0063 simplexNormalVectorMatrix = [];
0064 simplexVolumeSum = 0;
0065     
0066 % adding simplex partitions
0067 for iSimplex = 1:size(seSimplexIndexMatrix, 1)
0068     
0069     % that's the new index for the simplex partition
0070     thisSimplexIndex = length(dec.simplexPartition) + 1;
0071     
0072     % thats the matrix of structuring elements
0073     thisSeMatrix = dec.seScaledToClose(seIndexVector(...
0074         seSimplexIndexMatrix(iSimplex, :)), :);
0075     
0076     % sometimes, it happens that the triangulation is not full. We
0077     %   ignore these components.
0078     if rank(thisSeMatrix, zerotol) < dec.nC
0079         continue
0080     end
0081     
0082     % For simplex partitions, the hrep can be generated directly by inversion
0083     %   of the matrix of structuring elements. For a unified partition, this is
0084     %   not possible. How does this work?
0085     % The points of the cone is described by:
0086     %     {hc = se' * lam | lam >= 0}
0087     %   where the equation can be converted to:
0088     %     lam = inv(se') * hc
0089     %   to obtain the corresponding H-representation:
0090     %     {hc | -inv(se') * hc <= 0}
0091     
0092     % generate hrep of partition
0093     dec.simplexPartition{thisSimplexIndex}.A = -1*inv(thisSeMatrix');
0094     dec.simplexPartition{thisSimplexIndex}.h = zeros([dec.nC, 1]);
0095     dec.simplexPartition{thisSimplexIndex} = normalizeHrep(...
0096         dec.simplexPartition{thisSimplexIndex});
0097     
0098     % compose the matrix of all normal vectors
0099     simplexNormalVectorMatrix = [simplexNormalVectorMatrix; ...
0100         dec.simplexPartition{thisSimplexIndex}.A];
0101     
0102     % add structuring element indices
0103     dec.simplexPartition{thisSimplexIndex}.seIndexVector = ...
0104         sort(seIndexVector(seSimplexIndexMatrix(iSimplex, :)));
0105     
0106     % obtain scaled volume
0107     dec.simplexPartition{thisSimplexIndex}.volumeScaled = ...
0108         computeVolumeVrep(...
0109         struct('V', [thisSeMatrix; zeros([1, dec.nC])]) ...
0110         ) / dec.confinementCone.volume;
0111     simplexVolumeSum = simplexVolumeSum + ...
0112         dec.simplexPartition{thisSimplexIndex}.volumeScaled;
0113     
0114     % and we need to add this index
0115     dec.simplexToUnifiedVector(thisSimplexIndex) = thisUnifiedIndex;
0116     dec.simplexIsValidVector(thisSimplexIndex) = isValid;
0117     dec.simplexIsConfinedVector(thisSimplexIndex) = isConfined;
0118     
0119     % and we should add the mapping
0120     dec.simplexPartition{thisSimplexIndex}.mappingHcToLam = inv(dec.se(...
0121         dec.simplexPartition{thisSimplexIndex}.seIndexVector, : ...
0122         )');
0123 end
0124 
0125 % The simplex partitions were generated from a triangulation from which
0126 %   it is assumed that this was complete. As the unified partition is
0127 %   the union of all simplex partitions, we only have to filter for
0128 %   facets normals that have no opposing normal vector. These must be
0129 %   the outer normal vectors of the union.
0130 
0131 % generate unified hrep from vrep
0132 dec.unifiedPartition{thisUnifiedIndex}.A = eliminateRows(...
0133     simplexNormalVectorMatrix, -1*simplexNormalVectorMatrix, zerotol);
0134 dec.unifiedPartition{thisUnifiedIndex}.h = zeros(...
0135     [size(dec.unifiedPartition{thisUnifiedIndex}.A, 1), 1]);
0136 
0137 % add structuring element indices
0138 dec.unifiedPartition{thisUnifiedIndex}.seIndexVector = sort(seIndexVector);
0139 
0140 % add mapping to simplex partition (simply the last simplex mapping entry)
0141 dec.unifiedToSimplexVector(thisUnifiedIndex) = length(dec.simplexToUnifiedVector);
0142 dec.unifiedIsValidVector(thisUnifiedIndex) = isValid;
0143 dec.unifiedIsConfinedVector(thisUnifiedIndex) = isConfined;
0144 
0145 % obtain scaled volume
0146 dec.unifiedPartition{thisUnifiedIndex}.volumeScaled = ...
0147     computeVolumeVrep(...
0148     struct('V', [...
0149     dec.seScaledToClose(seIndexVector, :); ...
0150     zeros([1, dec.nC])...
0151     ]) ...
0152     ) / dec.confinementCone.volume;
0153 
0154 % add incidence information
0155 for iSeIndex = 1:length(seIndexVector)
0156     for jSeIndex = 1:length(seIndexVector)
0157         dec.seIncidenceMatrix(seIndexVector(iSeIndex), ...
0158             seIndexVector(jSeIndex)) = 1;
0159     end
0160 end
0161 
0162 end