Создайте случайный единичный вектор внутри определенной конической области

Question

Я ищу простой способ создания случайного единичного вектора, ограниченного конической областью. Происхождение всегда равно [0,0,0].

Мое решение до сих пор:

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 

Мой код петель до тех пор, случайным образом генерируется единичный вектор не находится внутри определенного конуса. Есть ли лучший способ сделать это?

Результирующее изображение из тестового кода ниже

результирующего распределения частот с использованием коды Ahmed Fasih (в комментарии). Интересно, как получить прямоугольное или нормальное распределение.

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

код

Тестирование:

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

Вот решение. Это основано на замечательном ответе на https://math.stackexchange.com/a/205589/81266. Я нашел этот ответ путем поиска «случайных точек на шаровой шапочке», после того как я узнал на Mathworld, что a spherical cap - это разрез 3-сферы с плоскостью.

Вот функция:

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 -z y; z 0 -x; -y x 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 

Этот код только реализует joriki’s answer on math.stackexchange, заполняя все детали, которые joriki опущенными.

Вот сценарий, который показывает, как его использовать.

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') 

Вот 3D сюжет 10'000 точек от 30 ° сферической крышки с центром вокруг [1; 1; 1] вектора:

Вот 120 ° сферическая крышка:

---

Теперь, если лет u посмотрите на гистограмму углов между этими случайными точками на coneDir = [1;1;1], вы увидите, что распределение искажено. Вот распределение:

Код для создания этого:

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)) 

Ну, это имеет смысл! Если вы создадите точки из сферического колпачка 120 ° вокруг coneDir, , конечно,, у крышки 1 ° будет очень мало таких образцов, тогда как полоса между колпачками 10 ° и 11 ° будет иметь гораздо больше точек. Поэтому мы хотим нормализовать количество точек под заданным углом по площади поверхности сферического колпачка при , что угол.

Вот функция, которая дает нам площадь поверхности сферического купола с радиусом R и угол в радианах theta (уравнение 16 на Mathworld’s spherical cap статьи):

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

Тогда мы можем нормализовать количество гистограммы для каждого бен (n выше) по разнице в площади поверхности сферических колпачков по краям бункера в:

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

на рисунке:

Это говорит нам «число случайных векторов, деленных на площади поверхности сферического сегмента между гистограммой бин краями». Это единообразно!

(NB Если вы сделаете это нормированную гистограмму для векторов, порожденных исходного кодом, используя выборку отклонения, то же самое имеет место.. Нормированная гистограмма равномерна Это просто, что отбор пробы отказа является дорогостоящим по сравнению с этим)

(Примечание 2: обратите внимание, что наивный способ выбора случайных точек на шаре - сначала генерирование азимута/углов возвышения, а затем преобразование этих сферических координат в декартовы координаты - не является хорошим, потому что он сгущает точки вблизи полюсов: Mathworld, example, example 2. Один из способов выбрать точки на всей сфере - это выборка из нормального распределения 3D: таким образом, вы не получите гроздь возле полюсов. Поэтому я считаю, что ваша оригинальная техника отлично что дает вам приятные, равномерно распределенные точки на сфере без каких-либо группировок. Этот алгоритм, описанный выше, также делает «правильную вещь» и должен избегать группировки. Внимательно оцените любые предлагаемые алгоритмы, чтобы избежать проблемы с группировкой рядом с полюсами.)

Answer 2

лучше использовать сферические координаты и преобразовать его в декартовых координатах:

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) 

, если все точки должны быть длиной 1, заданы r = единицы (N, 1);

Редактировать: Поскольку пересечение конуса с сферой образует круг, мы сначала создаем случайные точки внутри круга с raduis (45/2) в полярных координатах и ​​как @Ahmed Fasih прокомментировал, чтобы предотвратить концентрацию точек вблизи полюса. должны сначала преобразовать эти случайные точки, а затем преобразовать полярны декартовы координаты 2D с образованием x0 и y0

мы можем использовать x0 и y0, как фи & тета угол сферических координат и добавить coneDirtheta & coneDirphi как смещения для этих координат. затем преобразовать сферические в декартовые 3D-координаты