在一个定义的锥形区域内创建随机单位向量

Question

在一个定义的锥形区域内创建随机单位向量

9

我希望找到一种简单的方法来创建一个在锥形区域内受限的随机单位向量，起点始终为 [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

- Pedro77

2

你提供的代码并没有真正说明你的问题。所以请告诉我是否正确描述了你的问题：你有一个**确定**的圆锥区域（已知：_起点_、_方向_和_发散角度_），你想要生成一个随机向量（同样的_起点_），其方向包含在圆锥范围内？ - Hoki

余弦角的分布是否有任何要求？我可能认为角度应该是均匀的，但您的代码生成的向量相对于“coneDir”的角度非常倾斜。 - Ahmed Fasih

“angles with respect to coneDir are very skewed” 是什么意思？修正了角度计算，子空间不考虑向量方向。现在没问题了。 - Pedro77

1

能否更新一下问题描述，说明您希望结果具有什么统计特性？选择很多，这个问题非常有趣。 - Ahmed Fasih

1

如果你想要一个结果角度的正常或均匀分布，请使用这个作为起点。生成n个随机数，表示你想要的角度分布。每个数字（角度）将定义围绕锥形中心线的圆。然后在每个圆上生成一个随机点（或更多，但是对于每个圆使用相同数量的点）。所有获得的点将成为您的向量（原点为0），并且应该尊重您的分布。 - Hoki

显示剩余6条评论

2个回答

0

最好使用球坐标并将其转换为笛卡尔坐标：

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 = ones(N,1)。

编辑：由于圆锥与球体的交点形成一个圆，因此我们首先在极坐标下创建半径为（45/2）的圆内随机点，如@Ahmed Fasih所评论的那样，为了防止点集聚集在极点附近，我们应该先转换这些随机点，然后将极坐标转换为笛卡尔二维坐标以形成x0和y0。

我们可以使用x0和y0作为球面坐标的phi和theta角，并将coneDirtheta和coneDirphi作为偏移量添加到这些坐标中。然后将球面坐标转换为笛卡尔三维坐标。

- rahnema1

你能解释一下你的解决方案吗？生成的向量不形成一个圆锥，它更像是一个“矩形锥”约束区域。 - Pedro77

@Pedro77的答案已更新为“圆锥形”，我添加了一些解释。 - rahnema1

1

我认为这个解决方案存在与尝试通过均匀生成方位角和俯仰角（球面坐标的角度）来在球体上生成随机点的问题相同：请参见http://mathworld.wolfram.com/SpherePointPicking.html的第一段 - 如果您使用此技术，您的解决方案将会聚集在极点附近。 - Ahmed Fasih

@Ahmed Fasih 很有趣！谢谢，我编辑了答案。 - rahnema1

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Ahmed Fasih · Accepted Answer

以下是解决方案。它基于这个不错的回答：https://math.stackexchange.com/a/205589/81266。我通过在Google上搜索“球冠上的随机点”，找到了这个答案。当时我从Mathworld上学到，球冠是由平面削减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在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')

这是一个以向量[1; 1; 1]为中心的30度球冠周围的一万个点的3D图：

这是一个120度球冠：

现在，如果您查看在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))

嗯，这很有道理！如果您从围绕 coneDir 的 120° 球冠中生成点，则当然 1° 球冠将具有非常少的这些样本，而 10° 和 11° 球冠之间的条带将具有更多的点。因此，我们想做的是通过球冠表面积标准化给定角度处的点数。

这是一个函数，它给出了半径为R且弧度角为theta的球冠的表面积（参见Mathworld's spherical cap文章中的方程式16）：

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

然后，我们可以通过归一化每个区间（上面的n）的直方图计数来得到每个区间在其边缘球冠表面积差异的影响:

接着，我们可以对于每个区间（记为n）将其直方图计数进行标准化，方法是用相邻两个区间边缘球冠表面积之差将该区间的直方图计数除以:

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

该图表：

告诉我们“随机向量数量除以直方图条块之间球形部分的表面积”。这是均匀的！

（注意：如果您对原始代码生成的向量执行此归一化直方图，使用拒绝抽样也同样成立：归一化直方图是均匀的。只是相比于拒绝抽样，拒绝抽样昂贵。）

（注意2：请注意，通过先生成方位/俯仰角，然后将这些球面坐标转换为笛卡尔坐标系的方式选择球面上的随机点是不好的，因为它会使点靠近极点：Mathworld，示例，示例2。选择在整个球体上的点的一种方法是从三维正态分布中采样：这样你就不会在极点附近聚集了。所以我相信您的原始技术是完全适当的，可以为您提供漂亮、均匀分布的球面点，而不会出现任何聚集。上面描述的这个算法也做到了“正确的事情”，应该避免聚集问题。请仔细评估任何拟议的算法，以确保避免极点附近的聚集问题。）