2012-10-30 61 views
5

Làm thế nào để vẽ một Geosphere trong MATLAB?Làm thế nào để vẽ Geosphere trong MATLAB?

Bằng Geosphere Tôi có nghĩa là cách điểm cụ thể hóa trên hình cầu (Geosphere là, ví dụ trong 3Ds Max).

Trên hình ảnh dưới đây, nó được hiển thị Sphere (ở bên trái) và Geosphere (ở bên phải) enter image description here

Trong Matlab có một chức năng sphere, mà đưa ra một kết quả như vậy: enter image description here

Tôi cần có một hình ảnh về Địa quyển. Tôi có một ma trận Nx3, với tọa độ xyz của mọi điểm của Địa quyển.

UPDATE:

Tôi có vấn đề duy nhất với hiển thị (vẽ) địa quyển, bởi vì tôi đã có dữ liệu - tọa độ xyz của mỗi điểm. Đó là lý do tại sao câu trả lời của Gunther Struyf đã giúp tôi và tôi chấp nhận nó.

này tôi nhận với cách hiển thị: (địa quyển của yếu tố = 6, N = 362) enter image description here

Làm thế nào để nhận được điểm 3d của địa quyển? Tôi đã sử dụng thư viện miễn phí SPHERE_GRID để nhận điểm 3d sphere. (Trong thư viện có những cách khác nhau của quả cầu giải phóng).

Ngoài ra, để tính điểm Địa quyển, cảm ơn @Rody Oldenhuis để có câu trả lời bên dưới.

+0

bây giờ tôi quan tâm về cách bạn xây dựng điểm đó mảng: p –

+0

@GuntherStruy: xem câu trả lời của tôi :) –

Trả lời

5

Nếu bạn đã có các điểm, tôi nghĩ bạn có thể sử dụng giải pháp tương tự như here:

tri = convhull(xyz); 
trisurf(tri,xyz(:,1),xyz(:,2),xyz(:,3)); 
+0

Vâng, nó làm việc. Cảm ơn bạn –

5

Tôi có một chức năng nằm xung quanh đây có nghĩa là có thể tạo ra một tam giác của hình cầu lên đến độ chính xác tùy ý . Đó là một chức năng tôi đã dựa trên buildsphere bởi Giaccari Luigi, người thật đáng buồn đã biến mất khỏi Internet hoàn toàn (cùng với chức năng này).

Vì vậy, tôi sẽ đăng ở đây và trên tài khoản Trao đổi tệp của tôi. Khi được chấp thuận, tôi sẽ thay thế mã này bằng một liên kết.

Lưu ý rằng tôi có một số mẫu được tạo trước trên máy tính của mình, chức năng này phụ thuộc vào chức năng nào. Bạn sẽ phải xóa phần đó trước khi có thể chạy.

Để vẽ hình cầu được tạo: xem câu trả lời của 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 : [email protected] 
% 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 
Các vấn đề liên quan