computeCollectiveShapeStruct

 compute information to quickly operate on lambda-shapes

 THIS IS NO USER FUNCTION

• computeEdgeLengthFunctional compute function handles for edge lengths and vrep
• computeVolumeVrep compute the volume of a polytope in V-representation
• convertHrepToVrep convert the H-representation of a polytope into V-representation
• reduceHrep reduce a H-representation to non-redundant support hyperplanes

• fitShapeBnd fit a rounded shape to boundary data

0001 function collectiveStruct = computeCollectiveShapeStruct(hrep, zerotol)
0002 % compute information to quickly operate on lambda-shapes
0003 %
0004 % THIS IS NO USER FUNCTION
0005 
0006 % The elk-library: convex geometry applied to crystallization modeling.
0007 %   Copyright (C) 2013 Alexander Reinhold
0008 %
0009 % This program is free software: you can redistribute it and/or modify it
0010 %   under the terms of the GNU General Public License as published by the
0011 %   Free Software Foundation, either version 3 of the License, or (at your
0012 %   option) any later version.
0013 %
0014 % This program is distributed in the hope that it will be useful, but
0015 %   WITHOUT ANY WARRANTY; without even the implied warranty of
0016 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
0017 %   General Public License for more details.
0018 %
0019 % You should have received a copy of the GNU General Public License along
0020 %   with this program.  If not, see <http://www.gnu.org/licenses/>.
0021 
0022 %% initial structure
0023 collectiveStruct(1).hrep = hrep;
0024 collectiveStruct(1).lambda = 0;
0025 
0026 while ~isempty(collectiveStruct(end).hrep.A)
0027     %% ensure nonredundant representation
0028     collectiveStruct(end).hrep = reduceHrep(collectiveStruct(end).hrep);
0029     
0030     %% check emptyness
0031     % compute vrep with large infinity box
0032     collectiveStruct(end).vrep = ...
0033         convertHrepToVrep(collectiveStruct(end).hrep, [], 1e5);
0034     
0035     % polytope is not sufficiently closed if a vertex exceeds the following
0036     %   threshold; it is also considered final when its volume does not
0037     %   contain a pixel anymore
0038     if any(abs(collectiveStruct(end).vrep.V(:)) > 9e4) || ...
0039        (computeVolumeVrep(collectiveStruct(end).vrep) < 1e-3)
0040         % if not closed, there is no hrep and vrep equals the origin
0041         collectiveStruct(end).hrep = struct('A', [], 'h', []);
0042         collectiveStruct(end).vrep.V = [0 0];
0043     end
0044     
0045     %% compute vrep and edgeLength functionals
0046     [collectiveStruct(end).edgeLengthFunctional, ...
0047      collectiveStruct(end).vrepFunctional, ...
0048      collectiveStruct(end).incidenceMatrix] = ...
0049      computeEdgeLengthFunctional(collectiveStruct(end).hrep, ...
0050                                  collectiveStruct(end).vrep, zerotol);
0051     
0052     %% find next lambda
0053     % now we investigate the edge lengths of P(h-lam) which are calculated
0054     %   by:
0055     %     e = M*(h-lam)
0056     %   while we demand e>=0 and actually search for the largest lambda
0057     %   that fulfills this inequality:
0058     %     <maximize lam> with (M*1)*lam<=M*h.
0059     %   We can directly calculate the candidate lambda values for each edge
0060     %   and subequently have to take the smallest lambda
0061     nEdge = size(collectiveStruct(end).hrep.A, 1);
0062     leftVector = collectiveStruct(end).edgeLengthFunctional(ones(nEdge, 1));
0063     rightVector = collectiveStruct(end).edgeLengthFunctional(...
0064         collectiveStruct(end).hrep.h);
0065     lambdaVector = rightVector ./ leftVector;
0066     [thisLambda thisRemIndex] = min(lambdaVector);
0067     
0068     %% assign next hrep
0069     if  ~isempty(collectiveStruct(end).hrep.A)
0070         collectiveStruct(end+1).hrep = collectiveStruct(end).hrep;
0071         collectiveStruct(end).hrep.A(thisRemIndex, :) = [];
0072         collectiveStruct(end).hrep.h(thisRemIndex) = [];
0073         collectiveStruct(end).hrep.h = collectiveStruct(end).hrep.h - thisLambda;
0074         collectiveStruct(end).lambda = collectiveStruct(end-1).lambda + ...
0075             thisLambda;
0076     end
0077     
0078 end