compute the volume of a polytope in V-representation Syntax: volume = computeVolumeVrep(vrep) Compute the volume based on MATLAB function convhulln that itself relies on qhull. computeVolumeVrep(..., zerotol) allows to set the tolerance for underlying algorithms (e.g. determining flat polytopes) (Default in elkZerotol). See also: computeVolumeHrep, computeFacetAreaVrep, computeProjectionAreaVrep, computeFeretVrep, computeSupportValueVrep
0001 function volume = computeVolumeVrep(vrep, zerotol) 0002 % compute the volume of a polytope in V-representation 0003 % 0004 % Syntax: volume = computeVolumeVrep(vrep) 0005 % 0006 % Compute the volume based on MATLAB function convhulln that itself relies 0007 % on qhull. 0008 % 0009 % computeVolumeVrep(..., zerotol) allows to set the tolerance for 0010 % underlying algorithms (e.g. determining flat polytopes) (Default in 0011 % elkZerotol). 0012 % 0013 % See also: computeVolumeHrep, computeFacetAreaVrep, 0014 % computeProjectionAreaVrep, computeFeretVrep, computeSupportValueVrep 0015 0016 % The elk-library: convex geometry applied to crystallization modeling. 0017 % Copyright (C) 2012 Alexander Reinhold 0018 % 0019 % This program is free software: you can redistribute it and/or modify it 0020 % under the terms of the GNU General Public License as published by the 0021 % Free Software Foundation, either version 3 of the License, or (at your 0022 % option) any later version. 0023 % 0024 % This program is distributed in the hope that it will be useful, but 0025 % WITHOUT ANY WARRANTY; without even the implied warranty of 0026 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0027 % General Public License for more details. 0028 % 0029 % You should have received a copy of the GNU General Public License along 0030 % with this program. If not, see <http://www.gnu.org/licenses/>. 0031 0032 %% Input handling 0033 if ~exist('zerotol', 'var') 0034 zerotol = elkZerotol; 0035 end 0036 0037 %% Collect 1-dim and 2-dim 0038 % they are not properly treated by convhulln. 0039 if size(vrep.V, 2) == 0 0040 volume = 0; 0041 return; 0042 elseif size(vrep.V, 2) == 1 0043 volume = max(vrep.V) - min(vrep.V); 0044 return 0045 elseif size(vrep.V, 1) <= size(vrep.V, 2) 0046 % a polytope requires at least (nDim+1) extreme points 0047 volume = 0; 0048 return 0049 end 0050 0051 %% Collect obviously degenerates and simplices 0052 if size(vrep.V, 1) <= size(vrep.V, 2) 0053 volume = 0; 0054 return 0055 elseif size(vrep.V, 1) == size(vrep.V, 2) + 1 0056 volume = abs(1/factorial(size(vrep.V, 2)) * ... 0057 det(bsxfun(@minus, vrep.V(2:end, :), vrep.V(1,:)))... 0058 ); 0059 return 0060 end 0061 0062 %% Scaling 0063 [vrepScaled, scale, ~] = scaleVrep(vrep, zerotol, 1); 0064 0065 %% Checking, whether Vrep is full-dimensional 0066 % Sometimes it is hard for qhull to recognise flat polytopes. In n 0067 % dimensions n+1 'proper' vertices are required to form a simplex. One of 0068 % them can be moved to the origin at all times which leaves n verties 0069 % that simply must have the full rank to form a full-dimesnional 0070 % polytope. 0071 vrepTmp = moveVrep(vrepScaled, -1*vrepScaled.V(1,:)); 0072 if rank(vrepTmp.V, zerotol) < size(vrepTmp.V, 2) 0073 volume = 0; 0074 return 0075 end 0076 0077 %% Volume calculation 0078 %! debug: the enclosed try/catch is debug code, you may uncomment it to 0079 % have a look on the polytope that causes troubles. It should not % 0080 % be flat polytope. 0081 % unfortunatelly, convhulln sometimes throws errors on non-singular 0082 % polytopes so that the command is enclosed in a try/catch environment. 0083 try 0084 % Option 'Qs' selected because generic_paracetamol_15 throws an error 0085 % during decomposition otherwise. From this error, the following line 0086 % was contained: "use 'Qs' to search all points for the initial 0087 % simplex" 0088 [~, volume] = convhulln(vrepScaled.V, {'Qt', 'Qs'}); 0089 catch errMsg 0090 % try 0091 % % joggling input for qhull. This may be required for some 0092 % % examples. (http://www.qhull.org/html/qh-optq.htm#QJn) 0093 % [~, volume] = convhulln(vrep, {['QJ1e' num2str(log(zerotol)/log(10))]}); 0094 % catch 0095 0096 % if this does not work, save the case for further analysis 0097 if exist('log_computeVolumeVrep.mat', 'file') 0098 load log_computeVolumeVrep.mat 0099 logVrep{end+1} = vrep; 0100 else 0101 logVrep{1} = vrep; 0102 end 0103 save log_computeVolumeVrep.mat logVrep 0104 error('elk:polytope:computeMeasure', ... 0105 ['Something is wrong. Polytope stored in ' cd filesep 'log_computeVolumeVrep.mat']); 0106 0107 % end %/ 0108 end 0109 0110 %% consider scaling 0111 volume = volume * scale^(size(vrep.V, 2)); 0112 0113 end 0114