球面上的最佳轨迹

Question

我正在尝试找到一点在单位球面上跳跃到不同位置的参数方程，使得：

1. 每次跳跃很小（pi/4 < d < pi/2），跳跃范围很窄，例如[1.33, 1.34]
2. 该点能够尽可能快速且均匀地访问球面的大部分区域
3. 该点沿着尽可能不同的“方向向量”行进

这是我的尝试：

N = 3600;    % number of points
t = (1:N) * pi / 180;    % parameter
theta_sph = sqrt(2) * t * pi;    % first angle
phi_sph = sqrt(3) * t * pi;    % second angle
rho_sph = 1;    % radius
% Coordinates of a point on the surface of a sphere
x_sph = rho_sph * sin(phi_sph) .* cos(theta_sph);
y_sph = rho_sph * sin(phi_sph) .* sin(theta_sph);
z_sph = rho_sph * cos(phi_sph);

% Check length of jumps (it is intended that this is valid only for small jumps!!!)
aa = [x_sph(1:(N-1)); y_sph(1:(N-1)); z_sph(1:(N-1))];
bb = [x_sph(2:N); y_sph(2:N); z_sph(2:N)];
cc = cross(aa, bb);
d = rho_sph * atan2(arrayfun(@(n) norm(cc(:, n)), 1:size(cc,2)), dot(aa, bb));
figure
plot(d, '.')
figure
plot(diff(d), '.')

% Check trajectory on the surface of the sphere
figure
hh = 1;
h_plot3 = plot3(x_sph(hh), y_sph(hh), z_sph(hh), '-');
hold on
axis square
% axis off
set(gca, 'XLim', [-1 1])
set(gca, 'YLim', [-1 1])
set(gca, 'ZLim', [-1 1])
for hh = 1:N
  h_point3 = plot3(x_sph(hh), y_sph(hh), z_sph(hh), ...
      'o', 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'r');
  drawnow
  delete(h_point3)
  set(h_plot3, 'XData', x_sph(1:hh))
  set(h_plot3, 'YData', y_sph(1:hh))
  set(h_plot3, 'ZData', z_sph(1:hh))
end

编辑--> 有没有人能找到一条更规则的轨迹，可以更快地覆盖球面（即跳数最少）并更加均匀？所谓规则轨迹是指方向应该平滑变化，而不是突然改变。美学美感是一个额外的奖励。这些点应该尽可能均匀地分布在球面上。

- user2875617

---

你能进一步澄清标准吗？特别是第2点和第3点，以及你在最后所说的“常规”和“均匀”的含义是什么？ - aganders3

均匀分布意味着点应该遍布整个球面，彼此之间距离最远（这可以很好地测量）。当然，不正确的理解是：通过规则轨迹，我指的是一条不会太多变化的轨迹，即它应该非常平滑地改变方向，而不是急剧地改变方向。如果它在美学上令人愉悦，那就更好了 :-) - user2875617

你能澄清第三点吗？如果我忽略它，我期望沿着球体螺旋下降（这是您正在使用的路径）是最优的。只需调整极角经过的完整旋转次数即可控制跳跃大小和覆盖范围。如果您控制极速，您可以“交错”您的点，以便它们不会全部落在同一纬线上，因为似乎您更关心“点”的覆盖范围而不是连续轨迹。 - kalhartt

我不喜欢我所描绘的情况的原因是，点一直在访问南极和北极，即在两极的点的密度远高于其他任何地方。螺旋是否可以相应地改变？你将如何调整完整旋转的次数？你将如何控制跳跃大小？你将如何控制极速率？你能展示一个例子吗？ - user2875617

Answer 1

修改代码开头，引入对底层球体的旋转。这会产生一条轨迹，不会频繁返回极点。可能需要调整旋转速度才能看起来“好看”（只绕一个轴旋转可能更好看）。rot_angle1是绕x轴旋转，rot_angle2和rot_angle3是绕y和z轴旋转。也许这能给你一个想法！

N = 3600;    % number of points
t = (1:N) * pi / 180;    % parameter
theta_sph = sqrt(2) * t * pi;    % first angle
phi_sph = sqrt(3) * t * pi;    % second angle
rho_sph = 1;    % radius
rot_angle1 = sqrt(2) * t * pi;
rot_angle2 = sqrt(2.5) * t * pi;
rot_angle3 = sqrt(3) * t * pi;
% Coordinates of a point on the surface of a sphere
x_sph0 = rho_sph * sin(phi_sph) .* cos(theta_sph);
y_sph0 = rho_sph * sin(phi_sph) .* sin(theta_sph);
z_sph0 = rho_sph * cos(phi_sph);

x_sph1 = x_sph0;
y_sph1 = y_sph0.*cos(rot_angle1)-z_sph0.*sin(rot_angle1);
z_sph1 = y_sph0.*sin(rot_angle1)+z_sph0.*cos(rot_angle1);

x_sph2 = x_sph1.*cos(rot_angle2)+z_sph1.*sin(rot_angle2);
y_sph2 = y_sph1;
z_sph2 = -x_sph1.*sin(rot_angle2)+z_sph1.*cos(rot_angle2);

x_sph = x_sph2.*cos(rot_angle3)-y_sph2.*sin(rot_angle3);
y_sph = x_sph2.*sin(rot_angle3)+y_sph2.*cos(rot_angle3);
z_sph = z_sph2;

- David

Answer 2

我手头没有Matlab的副本，但我会发布我对你的曲线所做修改。

仅为明确起见，因为有无限多+1个不同的球面角定义，我将使用以下定义，与您的定义相反，但如果我尝试切换，我可能会犯错误。

- \ phi-从z轴的角度 - \ theta-在x-y平面上的投影角度。

参数化

让t是从0到pi（含）的N个等间距点的离散集。

\phi(t) = t
\theta = 2 * c * t

相当直接和简单，围绕球体的螺旋线在 \phi 和 \theta 中是线性的。 c 是代表在 \theta 中完全旋转的数量的常数，它不必是整数。

相邻点

在您的示例中，您使用 atan2(norm(cross....) 计算向量之间的角度，这很好，但并没有给出任何关于问题的见解。您的问题在球体表面上，利用这个事实。因此，我考虑点之间的距离 使用此公式。

现在找到相邻点，这些点出现在 t +- dt 和 \theta +- 2\pi 中，无论发生什么 t。

在第一种情况t +- dt中，很容易计算cos(gamma) = 1 - 2 c^2 sin^2(t) dt^2。由于sin^2(t)的影响，极点更加密集。理想情况下，您需要选择theta(t)和phi(t)，使得dtheta^2 * sin^2(phi)是恒定且最小的，以满足此情况。

第二种情况稍微复杂一些，并引起了我对“错开”点的评论。如果我们选择一个N，使得dtheta不能均匀地除以2pi，则在theta的完整旋转后，我无法直接位于先前点的正下方。 在这种情况下比较点之间的距离，使用delta t，使得c delta t = 1。然后，您有delta phi = delta t和delta theta = 2 c delta t - 2pi。根据您选择的c，delta phi可能太大或太小，无法使用小角度近似。

最后说明

很显然，c=0 是一个沿球体直线向下的直线。通过增加 c，可以增加“螺旋密度”，获得更多的覆盖面积。不过，这也会增加相邻点之间的距离。您需要为所选的 N 选择一个使上述两个距离公式大致相等的 c。

- kalhartt

---

我已经发布了我的原始代码的修改版本 - 你认为它是否包含了你的提示，还是我有什么地方忽略了？要在Matlab/Octave中进行编译，可以使用http://www.compileonline.com/execute_matlab_online.php。 - user2875617

我还添加了您建议的两点之间距离的公式。 - user2875617

Answer 3

我用C语言写了一个快速版本，对于给定的固定点行为非常良好。你可以在ideone上玩一下。如果你的浏览器支持WebGL（Chrome、Firefox），你可以将这些结果here粘贴到那里进行绘图。由于在推导公式时使用了一些积分近似，所以极点略有偏差，但除此之外很难找出其他缺陷。除了输出所需点数之外，没有其他需要调整的常数。

#include <stdio.h>
#include <math.h>

int main(void) {
    int i, numPoints = 200;
    double slope = sqrt(1.2) / sqrt(numPoints);
    for (i = 0; i < numPoints; i++) {
        double s = asin((double)i / (double)(numPoints - 1) * 2.0 - 1.0);
        double z = sin(s);
        double r = sqrt(1.0 - z * z);
        double ss = (2.0 * s + M_PI) / slope;
        double x = cos(ss) * r;
        double y = sin(ss) * r;
        printf("%lf,%lf,%lf,\"1\"\n", x, y, z);
    }
    return 0;
}

- Kaganar

Answer 4

我在这里编辑，因为代码很长。在David和kalhartt的提示下，我尝试了这个：

N = 3600;    % number of points
t = (1:N) * pi / 180;    % parameter
% theta_sph much faster than phi_sph to avoid overly visiting the poles
theta_sph = sqrt(20.01) * t * pi;    % first angle
phi_sph = sqrt(.02) * t * pi;    % second angle
rho_sph = 1;    % radius
% Coordinates of a point on the surface of a sphere
x_sph0 = rho_sph * sin(phi_sph) .* cos(theta_sph);
y_sph0 = rho_sph * sin(phi_sph) .* sin(theta_sph);
z_sph0 = rho_sph * cos(phi_sph);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use David's hint to rotate axes (but only the first now):
rot_angle1 = t * pi / 10;
rot_angle2 = 0;
rot_angle3 = 0;

x_sph1 = x_sph0;
y_sph1 = y_sph0.*cos(rot_angle1)-z_sph0.*sin(rot_angle1);
z_sph1 = y_sph0.*sin(rot_angle1)+z_sph0.*cos(rot_angle1);

x_sph2 = x_sph1.*cos(rot_angle2)+z_sph1.*sin(rot_angle2);
y_sph2 = y_sph1;
z_sph2 = -x_sph1.*sin(rot_angle2)+z_sph1.*cos(rot_angle2);

x_sph = x_sph2.*cos(rot_angle3)-y_sph2.*sin(rot_angle3);
y_sph = x_sph2.*sin(rot_angle3)+y_sph2.*cos(rot_angle3);
z_sph = z_sph2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Check length of jumps
aa = [x_sph(1:(N-1)); y_sph(1:(N-1)); z_sph(1:(N-1))];
bb = [x_sph(2:N); y_sph(2:N); z_sph(2:N)];
cc = cross(aa, bb);
d = rho_sph * atan2(arrayfun(@(n) norm(cc(:, n)), 1:size(cc,2)), dot(aa, bb));
figure
plot(d, '.')

% Check trajectory on the surface of the sphere
figure
hh = 1;
h_plot3 = plot3(x_sph(hh), y_sph(hh), z_sph(hh), '-');
hold on
axis square
% axis off
set(gca, 'XLim', [-1 1])
set(gca, 'YLim', [-1 1])
set(gca, 'ZLim', [-1 1])
for hh = 1:N
  h_point3 = plot3(x_sph(hh), y_sph(hh), z_sph(hh), ...
      'o', 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'r');
  drawnow
  delete(h_point3)
  set(h_plot3, 'XData', x_sph(1:hh))
  set(h_plot3, 'YData', y_sph(1:hh))
  set(h_plot3, 'ZData', z_sph(1:hh))
end

我认为现在比以前好多了！我发现两件事很重要：1. theta_sph 必须比 phi_sph 快得多，以避免经常访问两个相反的极点；2. 如果 theta_sph 比 phi_sph 快，那么你必须缓慢地旋转 rot_angle1 或 rot_angle2 中的任意一个，以获得不太混乱的轨迹。我仍然愿意接受任何其他改进建议。

- user2875617

Answer 5

clear all
close all

u = pi/2:(-pi/36):-pi/2;
v = 0:pi/36:2*pi;

nv = length(v);
nu = length(u);

f = myfigure(1);
ax = myaxes(f,1,1);

hold on

for aa = 1:nv
    tv = v(aa);
    for bb = 1:nu
        tu = u(bb);
        x = sin(tu)*cos(tv);
        y = cos(tu)*cos(tv);
        z = sin(tv);
        plot3(x,y,z,'*')
    end
end

编辑：myfigure和myaxes是我用于创建图形和坐标轴的函数

- Alex Byasse