matlab - 如何在 MATLAB 中绘制球体上的平滑连接图？

Question

单位球体上有某些点，我想沿着球体平滑地连接它们，我如何在 MATLAB 中做到这一点，因为当我使用 Matlab 中的 3dplot 函数执行此操作时，它只是使用直线连接点。

例如，第一象限有一个点，第八象限有第二个点，然后用直线将它们连接起来。不遵循弯曲的路径。

这些是 theta 的值:

theta = [ 80.0000   73.2995   65.7601   95.5007  100.4861   97.8834   94.0849     52.5174 74.4710  104.6674   52.7177   97.0538   75.7018   83.2817   97.5423   85.1797 84.2677  126.2296   81.1814   66.1376   91.6953  167.7085   46.5980   87.8220 113.4588  180.0000   80.7624   95.8623  115.0538   76.5773   61.9858  141.0402  109.9872   76.1273   84.4166   75.2734  110.4489   82.2434   96.8303  100.0815 73.2454   82.0755   64.6457   76.3510   87.7863  133.2706   86.1305   76.8670  86.3225   96.8016   49.2653  107.2900  145.9905   59.2158  107.7546  180.0000 93.9687   87.5474  103.1400  180.0000  136.8251  180.0000  106.2629  109.0069 ];

以及 phi 的值是:

phi = [ -90.0000  -78.5230  -51.6764   84.6854   58.1182  -75.9705   78.0541  -60.0560  88.8935  -84.6539  -44.1415  -86.7643   61.7764  -87.4767  -86.9440  -80.2459  -76.8752   88.9510   64.7297  -51.1245  -83.1606  -88.7280  -32.7110   81.0951  86.8393   -0.0000   52.6243  -88.7833  -75.4600   84.1374   79.8300  -86.7258  -65.8055   80.9829  -89.3172   57.1802  -80.6346   72.5277  -87.4452   74.2778  -86.1069   76.6124  -80.4604   89.2202   85.0649   89.2164  -79.0290   84.9961  -88.2301  -87.5064   50.4016   83.0830   82.4863  -50.8481   87.0335   -0.0000  88.4613   79.7583  -80.6474   -0.0000   80.0771   -0.0000   89.2428  -82.769 ];

这些可以很容易地处理吗

Answer 1

如果您希望 MATLAB 像这样沿着单位球体绘图，则需要指定其间的所有点，因为 MATLAB 只会用直线连接点。

为此，我们可以调整 Roger Stafford's great solution在 MATLAB Central 上绘制任意两个连续点之间的最短大圆路径。

使用以下函数我们可以做到这一点。我们将确定两个连续点之间的最短大圆路径，然后在它们之间进行插值以在单位圆上绘制直线

function plotOnSphere(x,y,z,varargin)

    %// Vectors representing each point
    xyz = [x(:), y(:), z(:)].';  %'

    %// One vector of the "first" points and one of the "next" points
    v1 = xyz(:, 1:end-1);
    v2 = xyz(:, 2:end);

    %// Cross product between the vectors of one point and the next
    cv1v2 = cross(v1, v2);

    %// Compute unit vector in the plane defined by v1 and v2
    v3 = normc(cross(cv1v2, v1));

    %// Figure out the range of the inner angle between v1 and v2
    nc = sqrt(sum(cv1v2.^2, 1));
    t = atan2(nc, dot(v1, v2, 1));

    %// Number of points to sample between any two points on the sphere
    nPoints = 100;

    %// Compute the interpolant
    V = zeros([nPoints, fliplr(size(v1))]);
    for k = 1:numel(t)
        T = linspace(0, t(k), 100);
        V(:,k,:) = (v1(:,k) * cos(T) + v3(:,k) * sin(T)).';    %'
    end

    %// Break the result out into x,y,z parts
    xx = V(:,:,1);
    yy = V(:,:,2);
    zz = V(:,:,3);

    %// Plot the lines
    h = plot3(xx(:), yy(:), zz(:), varargin{:});
    hold on

    %// Plot the original data points
    plot3(x,y,z, 'o', ...
        'Color', get(h, 'Color'), ...
        'Parent', get(h, 'Parent'), varargin{:});
end

如果我们应用您提供的输入数据。

theta = [ 80.0000   73.2995   65.7601   95.5007  100.4861   97.8834   94.0849     52.5174 74.4710  104.6674   52.7177   97.0538   75.7018   83.2817   97.5423   85.1797 84.2677  126.2296   81.1814   66.1376   91.6953  167.7085   46.5980   87.8220 113.4588  180.0000   80.7624   95.8623  115.0538   76.5773   61.9858  141.0402  109.9872   76.1273   84.4166   75.2734  110.4489   82.2434   96.8303  100.0815 73.2454   82.0755   64.6457   76.3510   87.7863  133.2706   86.1305   76.8670  86.3225   96.8016   49.2653  107.2900  145.9905   59.2158  107.7546  180.0000 93.9687   87.5474  103.1400  180.0000  136.8251  180.0000  106.2629  109.0069 ];
phi = [ -90.0000  -78.5230  -51.6764   84.6854   58.1182  -75.9705   78.0541  -60.0560  88.8935  -84.6539  -44.1415  -86.7643   61.7764  -87.4767  -86.9440  -80.2459  -76.8752   88.9510   64.7297  -51.1245  -83.1606  -88.7280  -32.7110   81.0951  86.8393   -0.0000   52.6243  -88.7833  -75.4600   84.1374   79.8300  -86.7258  -65.8055   80.9829  -89.3172   57.1802  -80.6346   72.5277  -87.4452   74.2778  -86.1069   76.6124  -80.4604   89.2202   85.0649   89.2164  -79.0290   84.9961  -88.2301  -87.5064   50.4016   83.0830   82.4863  -50.8481   87.0335   -0.0000  88.4613   79.7583  -80.6474   -0.0000   80.0771   -0.0000   89.2428  -82.769 ];

%// Convert to cartesian coordinates
[x,y,z] = sph2cart(deg2rad(theta), deg2rad(phi), 1);

figure;
plotOnSurface(x,y,z);

%// Plot a unit sphere for reference
sphere()
s = findall(gca, 'type', 'surf');
set(s, 'FaceColor', 'k', 'FaceAlpha', 0.01, 'EdgeAlpha', 0.1)