5

如何在matlab中绘制地球圈?

我的意思是球体上离散化点的方式Geosphere(Geosphere 是,例如在 3Ds Max 中)。

在下图中,显示为Sphere(左侧)和Geosphere(右侧) 在此处输入图像描述

在 Matlab 中有一个函数sphere,它给出了这样的结果: 在此处输入图像描述

我需要得到这样的地圈图像。我有一个矩阵 Nx3,其中包含 Geosphere 的每个点的 xyz 坐标。

更新:

我只在显示(绘制)地球圈时遇到了问题,因为我已经有了数据 - 每个点的 xyz 坐标。这就是为什么Gunther Struyf我的回答对我有帮助,我接受了。

这是我用这种显示方式得到的:(Geosphere of factor = 6, N=362) 在此处输入图像描述

如何获得地圈的 3d 点?我使用免费库SPHERE_GRID来获取 3d 球体的点。(在库中有不同的离散化方式)。

另外,为了计算 Geosphere 的点,感谢@Rody Oldenhuis下面的回答。

4

2 回答 2

5

如果您已经掌握了要点,我认为您可以使用与此处相同的解决方案:

tri = convhull(xyz);
trisurf(tri,xyz(:,1),xyz(:,2),xyz(:,3));
于 2012-10-30T13:22:18.310 回答
5

我在这里有一个函数,它能够生成任意精度的球体三角剖分。这是我基于buildsphereGiaccari Luigi 的一个功能,可悲的是,他已经完全从互联网上消失了(连同这个功能)。

所以,我会把它贴在这里,并在我的文件交换帐户上。一旦获得批准,我将用链接替换此代码。

请注意,我的计算机上有一些预先生成的模型,此功能依赖于这些模型。您必须先删除该部分,然后才能运行它。

要绘制生成的球体:请参阅 Gunther Struyf 的答案。

function [p, t] = TriSphere(N, R)
% TRISPHERE: Returns the triangulated model of a sphere using the
% icosaedron subdivision method.
%
% INPUT:
% N (integer number) indicates the number of subdivisions,
%   it can assumes values between 0-inf. The greater N the better will look
%   the surface but the more time will be spent in surface costruction and
%   more triangles will be put in the output model.
%
% OUTPUT:
% In p (nx3) and t(mx3) we can find points and triangles indexes
% of the model. The sphere is supposed to be of unit radius and centered in
% (0,0,0). To obtain spheres centered in different location, or with
% different radius, is just necessary a traslation and a scaling
% trasformation. These operation are not perfomed by this code beacuse it is
% extrimely convinient, in order of time perfomances, to do this operation
% out of the function avoiding to call the costruction step each time.
%
% NOTE:
% This function is more efficient than the matlab command sphere in
% terms of dimension fo the model/ accuracy of recostruction. This due to
% well traingulated model that requires a minor number of patches for the
% same geometrical recostruction accuracy. Possible improvement are possible
% in time execution and model subdivision flexibilty.
%
% EXAMPLE:
%
%  N=5;
%
%  [p,t] = TriSphere(N);
%
%  figure(1) axis equal hold on trisurf(t,p(:,1),p(:,2),p(:,3)); axis vis3d
%  view(3)

% Author: Giaccari Luigi Created:25/04/2009%
% For info/bugs/questions/suggestions : giaccariluigi@msn.com
% ORIGINAL NAME: BUILDSPHERE
%
% Adjusted by Rody Oldenhuis (speed/readability)

    % error traps
    error(nargchk(1,1,nargin));
    error(nargoutchk(1,2,nargout));
    if ~isscalar(N)
        error('Buildsphere:N_mustbe_scalar',...
            'Input N must be a scalar.');
    end    
    if round(N) ~= N
        error('Buildsphere:N_mustbe_scalar',...
            'Input N must be an integer value.');
    end

    % Coordinates have been taken from Jon Leech' files

    % Twelve vertices of icosahedron on unit sphere
    tau = 0.8506508083520400; % t   = (1+sqrt(5))/2, tau = t/sqrt(1+t^2)
    one = 0.5257311121191336; % one = 1/sqrt(1+t^2)  (unit sphere)    
    p = [
        +tau  +one  +0     % ZA
        -tau  +one  +0     % ZB
        -tau  -one  +0     % ZC
        +tau  -one  +0     % ZD
        +one  +0    +tau   % YA
        +one  +0    -tau   % YB
        -one  +0    -tau   % YC
        -one  +0    +tau   % YD
        +0    +tau  +one   % XA
        +0    -tau  +one   % XB
        +0    -tau  -one   % XC
        +0    +tau  -one]; % XD

    % Structure for unit icosahedron
    t = [  
         5  8  9 
         5 10  8 
         6 12  7 
         6  7 11 
         1  4  5 
         1  6  4 
         3  2  8 
         3  7  2 
         9 12  1 
         9  2 12 
        10  4 11 
        10 11  3 
         9  1  5 
        12  6  1 
         5  4 10 
         6 11  4 
         8  2  9 
         7 12  2 
         8 10  3 
         7  3 11 ];

    % possible quick exit
    if N == 0, return, end

    % load pre-generated trispheres (up to 8 now...)
    if N <= 8
        S = load(['TriSphere', num2str(N), '.mat'],'pts','idx');
        p = S.pts; t = S.idx; 
        if nargin == 2, p = p*R; end
        return
    else
        % if even more is requested (why on Earth would you?!), make sure to START 
        % from the maximum pre-loadable trisphere
        S = load('TriSphere8.mat','pts','idx');
        p = S.pts; t = S.idx; clear S; N0 = 10;
    end

    % how many triangles/vertices do we have? 
    nt = size(t,1); np = size(p,1); totp = np;    
    % calculate the final number of points    
    for ii=N0:N, totp = 4*totp - 6; end    
    % initialize points array
    p = [p; zeros(totp-12, 3)];

    % determine the appropriate class for the triangulation indices
    numbits   = 2^ceil(log(log(totp+1)/log(2))/log(2));
    castToInt = ['uint',num2str(numbits)];

    % issue warning when required
    if numbits > 64
        warning('TriSphere:too_many_notes',...
            ['Given number of iterations would require a %s class to accurately ',...
            'represent the triangulation indices. Using double instead; Expect ',...
            'strange results!']);
        castToInt = @double;
    else
        castToInt = str2func(castToInt);
    end

    % refine icosahedron N times
    for ii = N0:N
        % initialize inner loop
        told  = t;
        t = zeros(nt*4, 3);
        % Use sparse. Yes, its slower in a loop, but for N = 6 the size is
        % already ~10,000x10,000, growing by a factor of 4 with every
        % increasing N; its simply too memory intensive to use zeros().
        peMap = sparse(np,np); 
        ct    = 1;        
        % loop trough all old triangles        
        for j = 1:nt

            % some helper variables
            p1 = told(j,1);
            p2 = told(j,2);
            p3 = told(j,3);
            x1 = p(p1,1); x2 = p(p2,1); x3 = p(p3,1);
            y1 = p(p1,2); y2 = p(p2,2); y3 = p(p3,2);
            z1 = p(p1,3); z2 = p(p2,3); z3 = p(p3,3);

            % first edge
            % -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

            % preserve triangle orientation
            if p1 < p2, p1m = p1; p2m = p2; else p2m = p1; p1m = p2; end

            % If the point does not exist yet, calculate the new point            
            p4 = peMap(p1m,p2m);
            if p4 == 0
                np = np+1;
                p4 = np;
                peMap(p1m,p2m) = np;%#ok
                p(np,1) = (x1+x2)/2;
                p(np,2) = (y1+y2)/2; 
                p(np,3) = (z1+z2)/2;
            end

            % second edge
            % -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

            % preserve triangle orientation
            if p2 < p3; p2m = p2; p3m = p3; else p2m = p3; p3m = p2; end

            % If the point does not exist yet, calculate the new point
            p5 = peMap(p2m,p3m);
            if p5 == 0
                np = np+1;
                p5 = np;
                peMap(p2m,p3m) = np;%#ok
                p(np,1) = (x2+x3)/2; 
                p(np,2) = (y2+y3)/2; 
                p(np,3) = (z2+z3)/2;
            end

            % third edge
            % -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

            % preserve triangle orientation
            if p1 < p3; p1m = p1; p3m = p3; else p3m = p1; p1m = p3; end

            % If the point does not exist yet, calculate the new point            
            p6 = peMap(p1m,p3m);
            if p6 == 0
                np = np+1;
                p6 = np;
                peMap(p1m,p3m) = np;%#ok
                p(np,1) = (x1+x3)/2; 
                p(np,2) = (y1+y3)/2;
                p(np,3) = (z1+z3)/2;                
            end

            % allocate new triangles
            %   refine indexing
            %          p1
            %          /\
            %         /t1\
            %      p6/____\p4
            %       /\    /\
            %      /t4\t2/t3\
            %     /____\/____\
            %    p3    p5     p2            
            t(ct,1) = p1; t(ct,2) = p4; t(ct,3) = p6; ct = ct+1;            
            t(ct,1) = p4; t(ct,2) = p5; t(ct,3) = p6; ct = ct+1;
            t(ct,1) = p4; t(ct,2) = p2; t(ct,3) = p5; ct = ct+1;            
            t(ct,1) = p6; t(ct,2) = p5; t(ct,3) = p3; ct = ct+1;

        end % end subloop
        % update number of triangles
        nt = ct-1;         
    end % end main loop

    % normalize all points to 1 (or R)
    p = bsxfun(@rdivide, p, sqrt(sum(p.^2,2)));
    if (nargin == 2), p = p*R; end
    % convert t to proper integer class
    t = castToInt(t);    

end % funciton TriSphere
于 2012-10-30T13:35:06.687 回答