Utwórz wektor losowy w określonym obszarze stożkowym

Question

Poszukuję prostego sposobu na utworzenie losowego wektora jednostkowego ograniczonego przez obszar stożkowy. Pochodzenie jest zawsze [0,0,0].

Moje rozwiązanie do tej pory:

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 

Mój kod pętle aż losowo generowany wektor jednostkowy jest wewnątrz zdefiniowanego stożka. Czy jest lepszy sposób to zrobić?

Wynikowy obraz z kodem testowym poniżej

wynikowego rozkładu częstotliwości z wykorzystaniem Ahmed Fasih kod (w komentarzach). Zastanawiam się, jak uzyskać rozkład prostokątny lub normalny.

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

kod

Testowanie:

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

Oto rozwiązanie. Opiera się na cudownej odpowiedzi pod adresem https://math.stackexchange.com/a/205589/81266. Znalazłem tę odpowiedź, wyszukując "przypadkowe punkty na czapce kulistej", po tym, jak nauczyłem się w Mathworld, że a spherical cap jest cięciem 3-kuli z płaszczyzną.

Oto funkcja:

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 

Ten kod zaledwie narzędzia joriki’s answer on math.stackexchange, wypełnianie wszystkich szczegółów, które joriki pominięte.

Oto skrypt, który pokazuje, jak z niego korzystać.

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

Tutaj jest wykresem 3D 10'000 punktów od 30 ° kulistej skupione wokół [1; 1; 1] wektorze:

oto 120 ° Kuliste:

---

Teraz, jeśli yo Spójrz na histogram kątów między tymi przypadkowymi punktami na coneDir = [1;1;1], zobaczysz, że rozkład jest przekrzywiony. Oto podział:

kod, aby wygenerować to:

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

No to sens! Jeśli wygenerujesz punkty z kulistej nasadki 120 ° wokół coneDir, oczywiście, czapeczka 1 ° będzie miała bardzo niewiele z tych próbek, podczas gdy pasek pomiędzy czapeczkami 10 ° i 11 ° będzie miał o wiele więcej punktów. Tak więc chcemy znormalizować liczbę punktów pod danym kątem przez powierzchnię czaszy kulistej pod kątem pod kątem.

Oto funkcja, która daje nam pole powierzchni kulistej o promieniu R i kąta w radianach theta (równanie 16 na Mathworld’s spherical cap artykułu):

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

Następnie możemy unormować liczbę histogramu dla każdego pojemnik (n wyżej) różnicy w powierzchni kuliste nasadki na krawędziach bin:

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

do rysunku:

To mówi nam, „liczba losowych wektorów podzielonej przez pole powierzchni sferycznej pomiędzy krawędziami segmentu bin Histogram”. To jest jednolite!

(NB Jeśli zrobisz to znormalizowany histogram dla wektorów generowanych przez oryginalnego kodu, stosując próbkowanie odrzucenia, to samo.. Znormalizowany histogram jest jednolita Tyle, że próbkowanie odrzucenie jest drogie w porównaniu do tego)

(Uwaga 2: zauważ, że naiwny sposób wybierania losowych punktów na kuli - przez pierwsze generowanie kątów azymutu/elewacji, a następnie przekształcenie tych współrzędnych sferycznych na współrzędne kartezjańskie - nie jest dobry, ponieważ skupia punkty w pobliżu biegunów: Mathworld, example, example 2.Jednym sposobem, aby wybrać punkty na całej kuli jest próbkowanie z rozkładu normalnego 3D: w ten sposób nie będziesz się wiązał z bliskimi biegunami, więc wierzę, że twoja oryginalna technika jest idealnie odpowiednie, dając ładne, równomiernie rozmieszczone punkty na kuli bez żadnego pęcznienia. Algorytm opisany powyżej również robi "właściwą rzecz" i powinien unikać skupiania. Dokładnie oceń wszystkie proponowane algorytmy, aby uniknąć problemów z pęczkiem w pobliżu bieguna.)

Answer 2

lepiej jest użyć współrzędnych sferycznych i przekonwertować go do współrzędnych kartezjańskich:

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) 

jeśli wszystkie punkty powinny mieć długość 1 zbiór r = ones (N, 1);

Edytować: ponieważ przecięcie stożka ze sferą tworzy koło najpierw tworzymy losowe punkty wewnątrz okręgu z raduią (45/2) we współrzędnych biegunowych i jako @Ahmed Fasih skomentowaliśmy, aby zapobiec koncentracji punktów w pobliżu bieguna należy najpierw przekształcić to losowe punkty, a następnie przekonwertować do polarnego współrzędne kartezjańskie 2D tworząc x0 i y0

możemy użyć x0 i y0 jako phi & kąta teta współrzędnych sferycznych i dodać coneDirtheta & coneDirphi jak zmiany tych współrzędnych. następnie konwertuj współrzędne sferyczne na kartezjańskie 3D