Erstellen zufälliger Einheitsvektor innerhalb einer definierten konischen Bereich

Question

Ich bin auf der Suche nach einer einfachen Möglichkeit zum Erstellen eines zufälligen Einheitsvektors durch eine konische Region beschränkt. Der Ursprung ist immer der [0,0,0].

Meine Lösung bis jetzt:

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 Codeschleifen, bis die erzeugte Zufallseinheitsvektor innerhalb des definierten Kegels ist. Gibt es einen besseren Weg, das zu tun?

resultierendes Bild von Testcode Balg

Resultierende Frequenzverteilung unter Verwendung Ahmed Fasih Code (in den Kommentaren). Ich frage mich, wie man eine rechteckige oder normale Verteilung bekommt.

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

Testing 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

Hier ist die Lösung. Es basiert auf der wunderbaren Antwort bei https://math.stackexchange.com/a/205589/81266. Ich fand diese Antwort, indem ich "zufällige Punkte auf der Kugelkalotte" googelte, nachdem ich auf Mathworld gelernt habe, dass a spherical cap dieser Schnitt einer 3-Sphäre mit einer Ebene ist.

Hier ist die Funktion:

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 

Dieser Code nur implementiert joriki’s answer on math.stackexchange, in allen Details füllen, die weggelassen joriki.

Hier ist ein Skript, das zeigt, wie man es benutzt.

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

Hier ist ein 3D-Plot von 10'000 Punkten aus dem 30 ° Kalotte um den Vektor [1; 1; 1] zentriert:

hier 120 ° Kalotte:

---

Jetzt, wenn yo Wenn Sie das Histogramm der Winkel zwischen diesen zufälligen Punkten bei coneDir = [1;1;1] betrachten, werden Sie sehen, dass die Verteilung verzerrt ist. Hier ist die Verteilung:

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

Nun, dies Sinn macht:

-Code dies zu generieren! Wenn Sie Punkte aus der 120 ° Kugelkalotte um coneDir, natürlich erzeugen, wird die 1 ° Kappe nur sehr wenige dieser Proben haben, während der Streifen zwischen den 10 ° und 11 ° Kappen viel mehr Punkte haben wird. Also, was wir tun wollen, ist die Anzahl der Punkte in einem bestimmten Winkel durch die Oberfläche der Kugelkalotte bei , dass Winkel zu normalisieren.

Hier ist eine Funktion, die uns die Oberfläche der Kalotte mit Radius R und Winkel in Radiant gibt theta (Gleichung 16 auf Mathworld’s spherical cap Artikel):

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

Dann wir die Histogrammzählwert für jede normalisieren Behälter (n oben) durch die Differenz der spezifischen Oberfläche der Kalotten an den Kanten des bin:

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

der Figur:

Dies sagt uns, „die Anzahl der Zufallsvektoren durch die Oberfläche des sphärischen Segments zwischen Histogrammbin Kanten geteilt“. Das ist einheitlich!

(NB Wenn Sie das normalisierte Histogramm für die von Ihrem ursprünglichen Code erzeugten Vektoren tun, Verwerfungsmethode verwenden, das gleiche gilt:.. Das normalisierte Histogramm einheitlich Es ist nur, dass die Ablehnung Sampling teuer ist im Vergleich zu diesem)

(NB 2: Beachten Sie, dass die naive Art, zufällige Punkte auf einer Kugel zu sammeln - indem zuerst Azimut-/Elevationswinkel erzeugt und dann diese Kugelkoordinaten in kartesische Koordinaten umgewandelt werden - nicht gut ist, da sie Punkte in der Nähe der Pole bündelt: Mathworld, example, example 2 Eine Möglichkeit, Punkte auf der gesamten Kugel zu sammeln, ist die 3D-Normalverteilung: Auf diese Weise werden Sie sich nicht in der Nähe von Stöcken sammeln, also glaube ich, dass Ihre ursprüngliche Technik perfekt ist geeignet, um Ihnen schöne, gleichmäßig verteilte Punkte auf der Kugel zu geben, ohne sich zu verbiegen. Dieser oben beschriebene Algorithmus tut auch das "Richtige" und sollte eine Bündelung vermeiden. Überprüfen Sie sorgfältig alle vorgeschlagenen Algorithmen, um sicherzustellen, dass das Problem der Bündelung in der Nähe der Pole vermieden wird.)

Answer 2

es besser ist, Kugelkoordinaten zu verwenden und es zu kartesischen Koordinaten umwandeln:

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) 

Wenn alle Punkte die Länge 1 haben sollen, setzen Sie r = Einsen (N, 1);

Edit: Da Schnittpunkt von Kegel mit Kugel einen Kreis bildet, erstellen wir zunächst zufällige Punkte innerhalb eines Kreises mit Raduis von (45/2) in Polarkoordinaten und als @Ahmed Fashi kommentierte, um Konzentration von Punkten in der Nähe des Poles zu verhindern soll zuerst diese zufälligen Punkte verwandeln, dann polar konvertieren kartesische 2D x0 und y0

bilden Koordinaten können wir x0 und y0 verwenden als & Theta-Winkel von Kugelkoordinaten phi und coneDirtheta & coneDirphi als Offsets zu diesen Koordinaten hinzuzufügen. dann sphärisch in kartesische 3D-Koordinaten konvertieren