Create a random unit vector inside a specific conical area

Question

I am looking for an easy way to create a random unit vector bounded by a conical region. The origin is always [0,0,0].

My solution so far:

function v = GetRandomVectorInsideCone(coneDir,coneAngleDegree) coneDir = normc(coneDir); ang = coneAngleDegree + 1; while ang > coneAngleDegree v = [randn(1); randn(1); randn(1)]; v = v + coneDir; v = normc(v); ang = atan2(norm(cross(v,coneDir)), dot(v,coneDir))*180/pi; end 

My code loops until a random generated unit vector is inside a specific cone. Is there a better way to do this?

The resulting image from the test code below

The resulting frequency distribution using the Ahmed Fasih code (in the comments). I wonder how to get a rectangular or normal distribution.

 c = [1;1;1]; angs = arrayfun(@(i) subspace(c, GetRandomVectorInsideCone(c, 30)), 1:1e5) * 180/pi; figure(); hist(angs, 50); 

Test Code:

 clearvars; clc; close all; coneDir = [randn(1); randn(1); randn(1)]; coneDir = [0 0 1]'; coneDir = normc(coneDir); coneAngle = 45; N = 1000; vAngles = zeros(N,1); vs = zeros(3,N); for i=1:N vs(:,i) = GetRandomVectorInsideCone(coneDir,coneAngle); vAngles(i) = subspace(vs(:,i),coneDir)*180/pi; end maxAngle = max(vAngles); minAngle = min(vAngles); meanAngle = mean(vAngles); AngleStd = std(vAngles); fprintf('v angle\n'); fprintf('Direction: [%.3f %.3f %.3f]^T. Angle: %.2fº\n',coneDir,coneAngle); fprintf('Min: %.2fº. Max: %.2fº\n',minAngle,maxAngle); fprintf('Mean: %.2fº\n',meanAngle); fprintf('Standard Dev: %.2fº\n',AngleStd); %% Plot figure; grid on; rotate3d on; axis equal; axis vis3d; axis tight; hold on; xlabel('X'); ylabel('Y'); zlabel('Z'); % Plot all vectors p1 = [0 0 0]'; for i=1:N p2 = vs(:,i); plot3ex(p1,p2); end % Trying to plot the limiting cone, but no success here :( % k = [0 1]; % [X,Y,Z] = cylinder([0 1 0]'); % testsubject = surf(X,Y,Z); % set(testsubject,'FaceAlpha',0.5) % N = 50; % r = linspace(0, 1, N); % [X,Y,Z] = cylinder(r, N); % % h = surf(X, Y, Z); % % rotate(h, [1 1 0], 90); 

plot3ex.m:

 function p = plot3ex(varargin) % Plots a line from each p1 to each p2. % Inputs: % p1 3xN % p2 3xN % args plot3 configuration string % NOTE: p1 and p2 number of points can range from 1 to N % but if the number of points are different, one must be 1! % PVB 2016 p1 = varargin{1}; p2 = varargin{2}; extraArgs = varargin(3:end); N1 = size(p1,2); N2 = size(p2,2); N = N1; if N1 == 1 && N2 > 1 N = N2; elseif N1 > 1 && N2 == 1 N = N1 elseif N1 ~= N2 error('if size(p1,2) ~= size(p1,2): size(p1,2) and/or size(p1,2) must be 1 !'); end for i=1:N i1 = i; i2 = i; if i > N1 i1 = N1; end if i > N2 i2 = N2; end x = [p1(1,i1) p2(1,i2)]; y = [p1(2,i1) p2(2,i2)]; z = [p1(3,i1) p2(3,i2)]; p = plot3(x,y,z,extraArgs{:}); end 

Answer 1

Here is the solution. It is based on a wonderful answer at https://math.stackexchange.com/a/205589/81266 . I found this answer by searching for "random points on a spherical cap" after I learned at Mathworld that a spherical cap is a section of a 3-sphere with a plane.

Here's the function:

 function r = randSphericalCap(coneAngleDegree, coneDir, N, RNG) if ~exist('coneDir', 'var') || isempty(coneDir) coneDir = [0;0;1]; end if ~exist('N', 'var') || isempty(N) N = 1; end if ~exist('RNG', 'var') || isempty(RNG) RNG = RandStream.getGlobalStream(); end coneAngle = coneAngleDegree * pi/180; % Generate points on the spherical cap around the north pole [1]. % [1] See https://math.stackexchange.com/a/205589/81266 z = RNG.rand(1, N) * (1 - cos(coneAngle)) + cos(coneAngle); phi = RNG.rand(1, N) * 2 * pi; x = sqrt(1-z.^2).*cos(phi); y = sqrt(1-z.^2).*sin(phi); % If the spherical cap is centered around the north pole, we're done. if all(coneDir(:) == [0;0;1]) r = [x; y; z]; return; end % Find the rotation axis `u` and rotation angle `rot` [1] u = normc(cross([0;0;1], normc(coneDir))); rot = acos(dot(normc(coneDir), [0;0;1])); % Convert rotation axis and angle to 3x3 rotation matrix [2] % [2] See https://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle crossMatrix = @(x,y,z) [0 -zy; z 0 -x; -yx 0]; R = cos(rot) * eye(3) + sin(rot) * crossMatrix(u(1), u(2), u(3)) + (1-cos(rot))*(u * u'); % Rotate [x; y; z] from north pole to `coneDir`. r = R * [x; y; z]; end function y = normc(x) y = bsxfun(@rdivide, x, sqrt(sum(x.^2))); end 

This code simply implements jorikis response to math.stackexchange , filling in all the details that are missing joriki.

Heres a script that shows how to use it.

 clearvars coneDir = [1;1;1]; coneAngleDegree = 30; N = 1e4; sol = randSphericalCap(coneAngleDegree, coneDir, N); figure;plot3(sol(1,:), sol(2,:), sol(3,:), 'b.', 0,0,0,'rx'); grid xlabel('x'); ylabel('y'); zlabel('z') legend('random points','origin','location','best') title('Final random points on spherical cap') 

Here is a 3D graph of 10,000 points from a 30 ° spherical cap centered around a vector [1; 1; 1] [1; 1; 1] [1; 1; 1] :

Heres a 120 ° spherical cap:

---

Now, if you look at the histogram of the angles between these random points at coneDir = [1;1;1] , you will see that the distribution is distorted. Heres the distribution:

Code to generate:

 normc = @(x) bsxfun(@rdivide, x, sqrt(sum(x.^2))); mysubspace = @(a,b) real(acos(sum(bsxfun(@times, normc(a), normc(b))))); angs = arrayfun(@(i) mysubspace(coneDir, sol(:,i)), 1:N) * 180/pi; nBins = 16; [n, edges] = histcounts(angs, nBins); centers = diff(edges(1:2))*[0:(length(n)-1)] + mean(edges(1:2)); figure('color','white'); bar(centers, n); xlabel('Angle (degrees)') ylabel('Frequency') title(sprintf('Histogram of angles between coneDir and random points: %d deg', coneAngleDegree)) 

Well, that makes sense! If you create points from a 120 ° spherical cap around coneDir , of course, there will be very few such samples per 1 ° cap, while the strip between the 10 ° and 11 ° caps will have much more points. Therefore, we want to normalize the number of points at a given angle over the surface area of ​​a spherical cap at that angle.

Here is a function that gives us the surface area of ​​a spherical cap with radius R and angle in theta radians (equation 16 on a Mathworlds spherical cap ):

 rThetaToH = @(R, theta) R * (1 - cos(theta)); rThetaToS = @(R, theta) 2 * pi * R * rThetaToH(R, theta); 

Then we can normalize the histogram score for each bin ( n above) by the difference in the surface area of ​​the spherical caps at the edges of the bin:

 figure('color','white'); bar(centers, n ./ diff(rThetaToS(1, edges * pi/180))) 

Picture:

This tells us "the number of random vectors divided by the surface area of ​​the spherical segment between the edges of the histogram bin." This is uniform!

(NB If you make this normalized histogram for vectors generated by your source code using a reject selection, the same thing happens: the normalized histogram is homogeneous. Just so that the reject selection is expensive compared to this.)

(NB 2: note that the naive way of selecting random points on a sphere - first creating azimuth / elevation angles and then converting these spherical coordinates to Cartesian coordinates - is not good because it will thicken the points near the poles: Mathworld , example , example 2. One way to select points across the entire sphere is to select from the normal 3D distribution: you won’t get close to the poles, so I think your original technique is perfect, giving you nice, evenly distributed points on the spheres e without any groupings. This algorithm described above also does the “right thing” and should avoid grouping. Carefully evaluate any proposed algorithms to ensure that the problem of grouping near the poles is avoided.)

Answer 2

it is better to use spherical coordinates and convert them to Cartesian coordinates:

 coneDirtheta = rand(1) * 2 * pi; coneDirphi = rand(1) * pi; coneAngle = 45; N = 1000; %perfom transformation preventing concentration of points around the pole rpolar = acos(cos(coneAngle/2*pi/180) + (1-cos(coneAngle/2*pi/180)) * rand(N, 1)); thetapolar = rand(N,1) * 2 * pi; x0 = rpolar .* cos(thetapolar); y0 = rpolar .* sin(thetapolar); theta = coneDirtheta + x0; phi = coneDirphi + y0; r = rand(N, 1); x = r .* cos(theta) .* sin(phi); y = r .* sin(theta) .* sin(phi); z = r .* cos(phi); scatter3(x,y,z) 

if all points must have a length of 1, set r = ones (N, 1);

Edit: since the intersection of the cone with the sphere forms a circle, we first create random points inside the circle with raduis (45/2) in polar coordinates and, as @Ahmed Fasih commented, to prevent the concentration of points near the pole, we must first transform these random points , then convert the polar to Cartesian 2D coordinates with the formation of x0 and y0

we can use x0 and y0 in the form of phi and theta-angle of spherical coordinates and add the cones Dirtheta and coneDirphi as offsets to these coordinates. then convert spherical to cartesian 3D coordinates