sample from n-dim cone (disc base) Syntax: [pointMatrix, probabilityVector] = sampleNormal(nSample, ray, radius) THIS IS NO USER FUNCTION
0001 function [positionMatrix probabilityVector] = sampleRay(... 0002 nSample, ray, radius) 0003 % sample from n-dim cone (disc base) 0004 % 0005 % Syntax: [pointMatrix, probabilityVector] = sampleNormal(nSample, ray, radius) 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 % A point of the cone can be sampled by first determine the level at 0027 % the ray, next the distance from the ray and finally a point at the 0028 % boundary for the direction. This topic is handled for n-dimensional 0029 % cases. 0030 % There must be linearly increasing more points at the base of the cone 0031 % than at the tip. The probability distribution for points at the ray 0032 % (l=0 -> tip, l=1 -> base) can be given by: 0033 % P(l) = 2*l 0034 % so that the cumulative density function is given by: 0035 % PC(l) = l^2 0036 % from which can be sampled using a uniform distribution for PC to 0037 % get l: 0038 % l = sqrt(PC) 0039 % The distance from the ray shall also be chosen in consideration of 0040 % uniformly distributed points in the cone. The sphere at the 0041 % corresponding radius increases its surface polynomial with r^(n-2) 0042 % so that the radius is chosen from the uniformly distributed 0043 % variable PC: 0044 % r = PC^(1/(n-2)) 0045 % Special considerations are required for n=1. 0046 % The direction that is given by a point of the boundary of the 0047 % crosssectional hypersphere is detemined by a formula that can be 0048 % found in D. E. Knuth "The art of computer Programming, vol. 2 0049 % Seminumerical algorithms". Therefore, (n-1) normally distributed 0050 % variables x are used to obtain the points of a (n-1)-sphere 0051 % embedded in the n-dimensional space: 0052 % y = [x / norm(x); 0]; 0053 % We that need a rotation matrix that transforms the n-th axis to the 0054 % ray that will be U and obtain the final random points: 0055 % pivotPosition = U*(r*l)*y + l*ray; 0056 0057 %% input handling 0058 % adjust to column vector 0059 ray = ray(:); 0060 % problem dimension 0061 nDim = length(ray); 0062 0063 %% required random numbers 0064 % radial distance from cone 0065 r = rand(1, nSample).^(1/(nDim-1)); 0066 % position on ray (length) 0067 l = rand(1, nSample).^(1/nDim); 0068 % normal distributed points (equally distributed over angle of cone if 0069 % projected to unit length) 0070 x = randn(nDim-1, nSample); 0071 0072 %% building a cone (direction to x-axis 0073 % ray along x-axis, radius 1 0074 positionMatrix = [l*norm(ray); bsxfun(@times, radius*(r)./sqrt(sum(x.*x, 1)).*l, x)]; 0075 rayUnit = zeros(nDim, 1); 0076 rayUnit(1) = 1; 0077 0078 %% rotation matrix (U) 0079 [U, S, V] = svd([ray, zeros(nDim, nDim-1)]); 0080 % correct direction of rotation - svd does not care about the sign of new 0081 % coordinates, but here we do. Tolerance is not essential, however. 0082 if ~all(isZero(ray - U*rayUnit*norm(ray), 1e-10*norm(ray))) 0083 U = -1*U; 0084 end 0085 positionMatrix = U*positionMatrix; 0086 0087 %% density 0088 volume = computeConeVolume(ray, radius); 0089 probabilityVector = positionMatrix(1,:)*0 + 1/volume;