compute data bounding cone in scaled space THIS IS NO USER FUNCTION
0001 function [rayVectorMatrix, normalVectorMatrix] = ... 0002 computeDataBoundingCone(grassStruct, positionMatrix) 0003 % compute data bounding cone in scaled space 0004 % 0005 % THIS IS NO USER FUNCTION 0006 0007 % The elk-library: convex geometry applied to crystallization modeling. 0008 % Copyright (C) 2013 Alexander Reinhold 0009 % 0010 % This program is free software: you can redistribute it and/or modify it 0011 % under the terms of the GNU General Public License as published by the 0012 % Free Software Foundation, either version 3 of the License, or (at your 0013 % option) any later version. 0014 % 0015 % This program is distributed in the hope that it will be useful, but 0016 % WITHOUT ANY WARRANTY; without even the implied warranty of 0017 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0018 % General Public License for more details. 0019 % 0020 % You should have received a copy of the GNU General Public License along 0021 % with this program. If not, see <http://www.gnu.org/licenses/>. 0022 0023 % To have a quick way of distinguishing valid grassmann coordinates for 0024 % which the range of the projection intersects with the data points, we 0025 % enclose the data points by a cone and use the extreme rays for 0026 % identifyiers. Note that the first PCA component contains the 'size' 0027 % of the particles s% To have a quick way of distinguishing valid grassmann coordinates for 0028 % which the range of the projection intersects with the data points, we 0029 % enclose the data points by a cone and use the extreme rays for 0030 % identifyiers. Note that the first PCA component contains the 'size' 0031 % of the particles so that no points are on the negative axis of this 0032 % PCA-coordinate. Hence, we first fit a polytope in the subspace of 0033 % xPca1 = 1 with all data points scaled to into that subspace. 0034 positionMatrixInPcaSpace = grassStruct.mapToPcaSpace * positionMatrix; 0035 % identify whether all data should be in positive or negative halfspace of 0036 % xPca: 0037 if sum(positionMatrixInPcaSpace(1, :) < 0) < ... 0038 sum(positionMatrixInPcaSpace(1, :) > 0) 0039 xPcaDirection = 1; 0040 else 0041 xPcaDirection = -1; 0042 end 0043 0044 % filter zero (or negative xPca1 samples) 0045 positionMatrixInPcaSpace(:, ... 0046 xPcaDirection*positionMatrixInPcaSpace(1,:) <= 0) = []; 0047 % scale all to xPca1 = +/-1 0048 pointMatrixForPolytope = bsxfun(@times, positionMatrixInPcaSpace, ... 0049 1./positionMatrixInPcaSpace(1, :)); 0050 % remove first coordinate (all are zero) to obtain points in subspace with 0051 % xPca1 = 1 0052 pointMatrixForPolytope(1, :) = []; 0053 0054 % Now, we generate some facet normal vectors. The simplest case is: 0055 % [eye(nC-1); -1*eye(nC-1)] 0056 % which gives 2*(nC-1) facets, representing an (nC-1)-hypercube giving 0057 % 2^(nC-1) extreme points. The (nC-1) hypercross gives 2*(nC-1) 0058 % additional extreme points with 2^(nC-1) hyperplanes. 0059 % These constraints are hopefully sufficient to enclose the data. 0060 hrep.A = [eye(grassStruct.nDim-1); -1*eye(grassStruct.nDim-1)]; 0061 % To add the normal vectors for the cross, we use a binary encoding where 0062 % the binary digits encodes the + and - signs for the normal vector. 0063 crossNormalMatrix = ones(grassStruct.nDim-1, 2^(grassStruct.nDim-1)); 0064 for iNormal = 1:(2^(grassStruct.nDim-1)) 0065 crossNormalMatrix(bitand(2.^(0:(grassStruct.nDim-1)), ... 0066 iNormal-1) ~= 0, iNormal) = -1; 0067 end 0068 hrep.A = [hrep.A; crossNormalMatrix']; 0069 % normalize normals to unit length 0070 hrep.A = bsxfun(@times, hrep.A, 1./sqrt(sum(hrep.A.^2, 2))); 0071 0072 % compute required facet distances 0073 hrep.h = max(hrep.A * pointMatrixForPolytope, [], 2); 0074 % compute extreme points 0075 vrep = convertHrepToVrep(hrep); 0076 % lift back to scaled-PCA space 0077 extremePointMatrix = [xPcaDirection*ones(1, size(vrep.V, 1)); ... 0078 vrep.V']; 0079 coneNormalVectorMatrix = [-1*xPcaDirection*hrep.h'; hrep.A']; 0080 % normalize new normal vectors 0081 coneNormalVectorMatrix = bsxfun(@times, coneNormalVectorMatrix, ... 0082 1./sqrt(sum(coneNormalVectorMatrix.^2, 1))); 0083 0084 % transform to original space 0085 rayVectorMatrix = grassStruct.mapToRealSpace * extremePointMatrix; 0086 rayVectorMatrix = bsxfun(@times, rayVectorMatrix, ... 0087 1./sqrt(sum(rayVectorMatrix.^2, 1))); 0088 0089 normalVectorMatrix = (coneNormalVectorMatrix' * ... 0090 grassStruct.mapToPcaSpace)';