compute the mean width of a polytope given H-rep and V-rep Syntax: meanWidth = computeMeanWidthCore(hrep, vrep, edgeStruct) This function is no user function! It is called by computeMeanWidthVrep and computeMeanWidthHrep
0001 function meanWidth = computeMeanWidthCore(hrep, vrep, edgeStruct) 0002 % compute the mean width of a polytope given H-rep and V-rep 0003 % 0004 % Syntax: meanWidth = computeMeanWidthCore(hrep, vrep, edgeStruct) 0005 % 0006 % This function is no user function! It is called by computeMeanWidthVrep 0007 % and computeMeanWidthHrep 0008 0009 % The elk-library: convex geometry applied to crystallization modeling. 0010 % Copyright (C) 2012 Alexander Reinhold 0011 % 0012 % This program is free software: you can redistribute it and/or modify it 0013 % under the terms of the GNU General Public License as published by the 0014 % Free Software Foundation, either version 3 of the License, or (at your 0015 % option) any later version. 0016 % 0017 % This program is distributed in the hope that it will be useful, but 0018 % WITHOUT ANY WARRANTY; without even the implied warranty of 0019 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0020 % General Public License for more details. 0021 % 0022 % You should have received a copy of the GNU General Public License along 0023 % with this program. If not, see <http://www.gnu.org/licenses/>. 0024 0025 %% Algorithm 0026 % cycle through edges and sum up the mean width. The basic equation to do 0027 % this can be found in [Schneider, 2008: "Convex Bodies: the 0028 % Brunn-Minkowski Theory", equation4.2.30]. It calculates the intrinsic 0029 % volume. For the three dimensional case it can be interpreted as 0030 % calculating 3xV(K,K,B) (a mixed volume) so that calculating the mean 0031 % width multiplied by 2pi can be reformulated to: 0032 % Sum up over all edges: [length of the edge] x [angle between the 0033 % adjacent facets] / [full angle] x [area of the unit disc] 0034 % which then is well interpretable with the behavior of mixed volumes or 0035 % the Steiner formula, respectively. For the general n-dimensional case 0036 % the formula reads: 0037 % Sum over all edges: [length of edge] x [external angle of edge] 0038 % where the external angle is defined in a way that does not easily 0039 % extend in higher dimensions as calculations on ball-sections are 0040 % required (Schneider, 2008 - page 100) 0041 meanWidth = 0; 0042 for iEdge = 1:size(edgeStruct.vertexIndexMatrix, 1) 0043 thisEdgeLength = norm(vrep.V(edgeStruct.vertexIndexMatrix(iEdge, 1), :) - ... 0044 vrep.V(edgeStruct.vertexIndexMatrix(iEdge, 2), :) ... 0045 ); 0046 % scalar product 0047 scalarProduct = hrep.A(edgeStruct.facetIndexMatrix(iEdge, 1), :) * ... 0048 hrep.A(edgeStruct.facetIndexMatrix(iEdge, 2), :)'; 0049 % prevent imaginary numbers from acos 0050 if scalarProduct > 1 0051 scalarProduct = 1; 0052 elseif scalarProduct < -1 0053 scalarProduct = -1; 0054 end 0055 % opening anlge 0056 thisAngle = acos(scalarProduct); 0057 % [length] * [angle] / 2 / pi * pi: 0058 meanWidth = meanWidth + thisEdgeLength * thisAngle / 2; 0059 end 0060 meanWidth = meanWidth / 2 / pi;