function [rayVectorMatrix, normalVectorMatrix] = ...
    computeDataBoundingCone(grassStruct, positionMatrix)
% compute data bounding cone in scaled space
%
% THIS IS NO USER FUNCTION

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2013 Alexander Reinhold
%
% This program is free software: you can redistribute it and/or modify it 
%   under the terms of the GNU General Public License as published by the 
%   Free Software Foundation, either version 3 of the License, or (at your 
%   option) any later version.
%
% This program is distributed in the hope that it will be useful, but 
%   WITHOUT ANY WARRANTY; without even the implied warranty of 
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
%   General Public License for more details.
%
% You should have received a copy of the GNU General Public License along 
%   with this program.  If not, see <http://www.gnu.org/licenses/>.

% To have a quick way of distinguishing valid grassmann coordinates for
%   which the range of the projection intersects with the data points, we
%   enclose the data points by a cone and use the extreme rays for
%   identifyiers. Note that the first PCA component contains the 'size'
%   of the particles s% To have a quick way of distinguishing valid grassmann coordinates for
%   which the range of the projection intersects with the data points, we
%   enclose the data points by a cone and use the extreme rays for
%   identifyiers. Note that the first PCA component contains the 'size'
%   of the particles so that no points are on the negative axis of this
%   PCA-coordinate. Hence, we first fit a polytope in the subspace of
%   xPca1 = 1 with all data points scaled to into that subspace.
positionMatrixInPcaSpace = grassStruct.mapToPcaSpace * positionMatrix;
% identify whether all data should be in positive or negative halfspace of
%   xPca:
if sum(positionMatrixInPcaSpace(1, :) < 0) < ...
        sum(positionMatrixInPcaSpace(1, :) > 0)
    xPcaDirection = 1;
else
    xPcaDirection = -1;
end

% filter zero (or negative xPca1 samples)
positionMatrixInPcaSpace(:, ...
    xPcaDirection*positionMatrixInPcaSpace(1,:) <= 0) = [];
% scale all to xPca1 = +/-1
pointMatrixForPolytope = bsxfun(@times, positionMatrixInPcaSpace, ...
    1./positionMatrixInPcaSpace(1, :));
% remove first coordinate (all are zero) to obtain points in subspace with
%   xPca1 = 1
pointMatrixForPolytope(1, :) = [];

% Now, we generate some facet normal vectors. The simplest case is:
%     [eye(nC-1); -1*eye(nC-1)]
%   which gives 2*(nC-1) facets, representing an (nC-1)-hypercube giving
%   2^(nC-1) extreme points. The (nC-1) hypercross gives 2*(nC-1)
%   additional extreme points with 2^(nC-1) hyperplanes.
% These constraints are hopefully sufficient to enclose the data.
hrep.A = [eye(grassStruct.nDim-1); -1*eye(grassStruct.nDim-1)];
% To add the normal vectors for the cross, we use a binary encoding where
%   the binary digits encodes the + and - signs for the normal vector.
crossNormalMatrix = ones(grassStruct.nDim-1, 2^(grassStruct.nDim-1));
for iNormal = 1:(2^(grassStruct.nDim-1))
    crossNormalMatrix(bitand(2.^(0:(grassStruct.nDim-1)), ...
        iNormal-1) ~= 0, iNormal) = -1;
end
hrep.A = [hrep.A; crossNormalMatrix'];
% normalize normals to unit length
hrep.A = bsxfun(@times, hrep.A, 1./sqrt(sum(hrep.A.^2, 2)));

% compute required facet distances
hrep.h = max(hrep.A * pointMatrixForPolytope, [], 2);
% compute extreme points
vrep = convertHrepToVrep(hrep);
% lift back to scaled-PCA space
extremePointMatrix = [xPcaDirection*ones(1, size(vrep.V, 1)); ...
    vrep.V'];
coneNormalVectorMatrix = [-1*xPcaDirection*hrep.h'; hrep.A'];
% normalize new normal vectors
coneNormalVectorMatrix = bsxfun(@times, coneNormalVectorMatrix, ...
    1./sqrt(sum(coneNormalVectorMatrix.^2, 1)));

% transform to original space
rayVectorMatrix = grassStruct.mapToRealSpace * extremePointMatrix;
rayVectorMatrix = bsxfun(@times, rayVectorMatrix, ...
    1./sqrt(sum(rayVectorMatrix.^2, 1)));

normalVectorMatrix = (coneNormalVectorMatrix' * ...
    grassStruct.mapToPcaSpace)';