2012-10-30 19 views
5

Wie zeichne ich eine Geosphäre in Matlab?Wie zeichne Geosphäre in Matlab?

Durch Geosphere meine ich die Art der Diskretisierung Punkte auf einer Kugel (Geosphäre ist, zum Beispiel in 3Ds Max).

Auf dem Bild unten ist es Sphere gezeigt (links) und Geosphere (rechts) enter image description here

In Matlab gibt es eine Funktion sphere, die ein solches Ergebnis gibt: enter image description here

Ich muss ein solches Bild einer Geosphäre bekommen. Ich habe eine Matrix Nx3 mit xyz Koordinaten von jedem Punkt der Geosphäre.

UPDATE:

Ich hatte Problem nur mit der Anzeige (Zeichnung) Geosphäre, weil ich bereits Daten - xyz Koordinaten von jedem Punkt. Deshalb hat Gunther Struyf 's Antwort mir geholfen und ich akzeptierte es.

Das habe ich mit einer solchen Art und Weise bekam der Anzeige: (Geosphäre des Faktors = 6, N = 362) enter image description here

Wie 3D-Punkte der Geosphäre zu bekommen? Ich benutzte freie Bibliothek SPHERE_GRID, um Punkte der Kugel 3d zu erhalten. (In der Bibliothek gibt es verschiedene Arten der Diskretisierungssphäre).

Auch für die Berechnung von Punkten der Geosphäre, danke @Rody Oldenhuis für die Antwort unten.

+0

jetzt bin ich daran interessiert, wie Sie diesen Punkt Array gebaut: p –

+0

@GuntherStruy: siehe meine Antwort :) –

Antwort

5

Wenn Sie bereits die Punkte haben, ich glaube, Sie die gleiche Lösung wie here verwenden:

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

Ja, es hat funktioniert. Vielen Dank –

5

Ich habe hier eine Funktion liegen, die eine Triangulation der Kugel bis zu beliebiger Genauigkeit zu erzeugen in der Lage, . Es ist eine Funktion, die ich auf buildsphere von Giaccari Luigi basiert, der leider komplett aus dem Internet verschwunden ist (zusammen mit dieser Funktion).

Also, ich werde es hier und auf meinem File Exchange-Konto posten. Sobald dies genehmigt wird, ersetze ich diesen Code durch einen Link.

Beachten Sie, dass ich einige vorgenerierte Modelle auf meinem Computer habe, von denen diese Funktion abhängt. Sie müssen diesen Abschnitt entfernen, bevor Sie ihn ausführen können.

Zum Zeichnen der erzeugten Kugel: siehe Gunther Struyfs Antwort.

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 
Verwandte Themen