sample from n-dim polygon Syntax: [pointMatrix, probabilityVector] = sampleNormal(nSample, poly) THIS IS NO USER FUNCTION
0001 function [positionMatrix probabilityVector] = samplePolygon(... 0002 nSample, poly) 0003 % sample from n-dim polygon 0004 % 0005 % Syntax: [pointMatrix, probabilityVector] = sampleNormal(nSample, poly) 0006 % 0007 % THIS IS NO USER FUNCTION 0008 0009 % The elk-library: convex geometry applied to crystallization modeling. 0010 % Copyright (C) 2013 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 % Well, the polygon is triangulated and the volume for each simplex is 0027 % computed. For each sample point, a continuous random variable between 0 0028 % and 1 is sampled to select the simplex proportional to the volume. Once 0029 % the simplex is selected, nDim uniform random variables are selected 0030 % between 0 and 1. They are sorted and 0 is prepended giving a vector: 0031 % [r_0=0, r_1, ..., r_nDim]. 0032 % Now, we have the differences: 0033 % s_i = r_i - r_(i-1) 0034 % being uniformly distributed in the simplex: 0035 % A = { [x_1, ..., x_nDim]: x_i >= 0, sum(x_i) <= 1 } 0036 % A coordinate transformation converting the base vectors of the 0037 % coordinate system to the base vectors [v_(i) - v_(nDim+1)] gives the 0038 % right shaped simplex. And shifting it to v_(nDim+1) also provides the 0039 % right position. 0040 % Source: Luc Devroye, "Non-Uniform Random Variate Generation", 1986, 0041 % Springer, New York: section 2.1 on page 207. 0042 0043 %% input handling 0044 % ensure we have a vrep to work with 0045 infoPoly = identifyPolytope(poly); 0046 if strcmp(infoPoly.rep, 'Vrep') 0047 vrep = poly; 0048 else 0049 vrep = convertHrepToVrep(poly); 0050 end 0051 % obtain triangulation 0052 vrepTri = triangulateVrep(vrep); 0053 nSimplex = length(vrepTri); 0054 % obtain volumes 0055 vrepVolumeVector = zeros(1, nSimplex); 0056 for iSimplex = 1:nSimplex 0057 vrepVolumeVector(iSimplex) = computeVolumeVrep(vrepTri{iSimplex}); 0058 end 0059 0060 vrepTotalVolume = sum(vrepVolumeVector); 0061 vrepSelectionVector = cumsum(vrepVolumeVector) / vrepTotalVolume; 0062 0063 %% obtain simplex indices 0064 % uniform random variables 0065 randVector = rand(1, nSample); 0066 % accumulate indices 0067 simplexIndexVector = ones(1, nSample); 0068 for iSimplex = 2:nSimplex 0069 thisFilter = randVector > vrepSelectionVector(iSimplex-1); 0070 simplexIndexVector(thisFilter) = simplexIndexVector(thisFilter) + 1; 0071 end 0072 0073 %% obtain uniform points in stanard simplex 0074 % uniform random variables (columns for samples) 0075 randomMatrix = rand(infoPoly.dim, nSample); 0076 % sort rows and prepend 0 0077 randomMatrix = [zeros(1, nSample); 0078 sort(randomMatrix, 1)]; 0079 % obtain uniform points in standard simplex by difference 0080 positionMatrix = randomMatrix(2:end, :) - randomMatrix(1:(end-1), :); 0081 0082 %% perform transformations 0083 % as nSample is usually high while nSimplex shoud be low, we iterate over 0084 % the different simplices and, hence, over the different coordinate 0085 % transformations 0086 for iSimplex = 1:nSimplex 0087 % identify affected sample points 0088 thisFilter = simplexIndexVector == iSimplex; 0089 0090 % get coordinate transformation 0091 thisSimplexCoord = vrepTri{iSimplex}.V; 0092 targetBase = thisSimplexCoord(1:(end-1), :) - ... 0093 ones(infoPoly.dim, 1)*thisSimplexCoord(end, :); 0094 thisTrafo = targetBase'; 0095 0096 % apply trafo 0097 positionMatrix(:, thisFilter) = thisTrafo * positionMatrix(:, thisFilter); 0098 0099 % moving sample points properly 0100 positionMatrix(:, thisFilter) = positionMatrix(:, thisFilter) + ... 0101 thisSimplexCoord(end, :)' * ones(1, sum(thisFilter)); 0102 end 0103 0104 %% density 0105 probabilityVector = ones(1, nSample) / vrepTotalVolume;