function [positionMatrix probabilityVector] = samplePolygon(...
    nSample, poly)
% sample from n-dim polygon
%
% Syntax: [pointMatrix, probabilityVector] = sampleNormal(nSample, poly)
%
% 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/>.

%% Algorithm
% Well, the polygon is triangulated and the volume for each simplex is
%   computed. For each sample point, a continuous random variable between 0
%   and 1 is sampled to select the simplex proportional to the volume. Once
%   the simplex is selected, nDim uniform random variables are selected
%   between 0 and 1. They are sorted and 0 is prepended giving a vector:
%      [r_0=0, r_1, ..., r_nDim].
%   Now, we have the differences:
%       s_i = r_i - r_(i-1)
%   being uniformly distributed in the simplex:
%       A = { [x_1, ..., x_nDim]: x_i >= 0, sum(x_i) <= 1 }
%   A coordinate transformation converting the base vectors of the
%   coordinate system to the base vectors [v_(i) - v_(nDim+1)] gives the
%   right shaped simplex. And shifting it to v_(nDim+1) also provides the
%   right position.
% Source: Luc Devroye, "Non-Uniform Random Variate Generation", 1986, 
%         Springer, New York: section 2.1 on page 207.

%% input handling
% ensure we have a vrep to work with
infoPoly = identifyPolytope(poly);
if strcmp(infoPoly.rep, 'Vrep')
    vrep = poly;
else
    vrep = convertHrepToVrep(poly);
end
% obtain triangulation
vrepTri = triangulateVrep(vrep);
nSimplex = length(vrepTri);
% obtain volumes
vrepVolumeVector = zeros(1, nSimplex);
for iSimplex = 1:nSimplex
    vrepVolumeVector(iSimplex) = computeVolumeVrep(vrepTri{iSimplex});
end

vrepTotalVolume = sum(vrepVolumeVector);
vrepSelectionVector = cumsum(vrepVolumeVector) / vrepTotalVolume;

%% obtain simplex indices
% uniform random variables
randVector = rand(1, nSample);
% accumulate indices
simplexIndexVector = ones(1, nSample);
for iSimplex = 2:nSimplex
    thisFilter = randVector > vrepSelectionVector(iSimplex-1);
    simplexIndexVector(thisFilter) = simplexIndexVector(thisFilter) + 1;
end

%% obtain uniform points in stanard simplex
% uniform random variables (columns for samples)
randomMatrix = rand(infoPoly.dim, nSample);
% sort rows and prepend 0
randomMatrix = [zeros(1, nSample);
                sort(randomMatrix, 1)];
% obtain uniform points in standard simplex by difference
positionMatrix = randomMatrix(2:end, :) - randomMatrix(1:(end-1), :);

%% perform transformations
% as nSample is usually high while nSimplex shoud be low, we iterate over
%   the different simplices and, hence, over the different coordinate
%   transformations
for iSimplex = 1:nSimplex
    % identify affected sample points
    thisFilter = simplexIndexVector == iSimplex;
    
    % get coordinate transformation
    thisSimplexCoord = vrepTri{iSimplex}.V;
    targetBase = thisSimplexCoord(1:(end-1), :) - ...
        ones(infoPoly.dim, 1)*thisSimplexCoord(end, :);
    thisTrafo = targetBase';
    
    % apply trafo
    positionMatrix(:, thisFilter) = thisTrafo * positionMatrix(:, thisFilter);
    
    % moving sample points properly
    positionMatrix(:, thisFilter) = positionMatrix(:, thisFilter) + ...
        thisSimplexCoord(end, :)' * ones(1, sum(thisFilter));    
end

%% density
probabilityVector = ones(1, nSample) / vrepTotalVolume;