How to plot a solid ellipsoid and get the xyz coordinates of all its points?

Hello, I would like to plot a 3D solid ellipsoid. The ellipsoid(...) function only allows me to plot the surface of the ellipsoid, but I would like to have the inside volume filled too. How would this be possible? Thank you in advance!

Answer 1

[Rik]

You can use the coordinates to convert it to a patch with faces and vertices, after which you should be able to modify intriangulation (FEX) or mesh voxelisation (FEX).

I have written a comment to keep in mind if you use the second:

I came across a weird thing, which I want to let you know about. I have delineations in .wrl format that are read in to a triangulated mesh. The z-coordinates of the vertices are all discrete multiples apart, as they follow the slice thickness from the CT. When I then use this script, it may happen that slices are returned empty when the values in gridZ happen to be equal to one of the z-values. This is not always the case. The fix I used is simply using a tiny offset in the z-direction, but the artifact took me by surprise.

If you plan on using intriangulation, check if you need to use heavytest, because for me that reduced the error from 5% to 0.

[zhen]

Hi, I think there may be a misunderstanding. I currently only have the dimensions (semi-axes lengths) of the ellipsoid, and (1) I would like to have a 3D plot of full/solid ellipsoid and not just a surface plot. From the plot, (2) I would like to extract its xyz coordinates.

[Rik]

That makes it even easier (as long as you're not picky about rotation).

That is even easier if you don't care about rotation.

%set semi-axis values
r=[5 5 3];
r2=r.^2;
Nsteps=20;
%generate coordinate grid
[X,Y,Z]=meshgrid(...
    linspace(-r(1),r(1),Nsteps),...
    linspace(-r(2),r(2),Nsteps),...
    linspace(-r(3),r(3),Nsteps));
%apply equation (https://en.wikipedia.org/wiki/Ellipsoid)
[X2,Y2,Z2]=deal(X.^2,Y.^2,Z.^2);
bin= (X2/r2(1) + Y2/r2(2) + Z2/r2(3))<=1;

[wim]

How can I visualize the points in bin?

[Rik]

It is a logical volume, so most volume viewing tools would work. You could also treat it as point cloud, although you would need to extract the coordinates with something like findND.

[Walter Roberson]

findND does accept logical.

if ~(isnumeric(X) || islogical(X)) || numel(X)==0
    error('HJW:findND:FirstInput',...
        'Expected first input (X) to be a non-empty numeric or logical array.')
end

[Rik]

The code Walter quoted clearly shows findND does accept logicals as input, so you will have to show the error message. Did you use the Add-On Manager to install it? If not, did you make sure you put the function file in a folder on your path?

There is an even easier way to extract the point cload coordinates: just use bin on the coordinates.

[X2,Y2,Z2]=deal(X.^2,Y.^2,Z.^2);
bin= (X2/r2(1) + Y2/r2(2) + Z2/r2(3))<=1;
Xe=X2(bin);
Ye=Y2(bin);
Ze=Z2(bin);

[wim]

Sorry my mistake. I have deleted my previous comment!

[wim]

A follow up question on the intriangulation function.

• My ellipsoid with prinicipal axes r1,r2,r3, rotated by the euler-angles theta,psi,phi, and it's center translated by Cx,Cy,Cz.
• I generate surface coordinates of the ellipsoid using ellipsoid.
• Then I used surf2patch to obtain the faces and the vertices of my ellipsoid.
• I generate a testpoints on a meshgrid using ndgrid and test these using intriangulation

I am unsure why intriangulation generates points outside my ellipsoid (see attached picture). I have attached my function and running script below.

% running script 
clc; clear all; close all 
% principal lengths 
r1 = 10; 
r2 = 10; 
r3 = 10; 
% center of ellipsoid 
Cx = 10; 
Cy = 20; 
Cz = 30; 
% euler angles 
theta = 2*pi*rand(1);
psi = 2*pi*rand(1); 
phi = 2*pi*rand(1); 
% calling function 
ExceedEllipsoid(r1,r2,r3,Cx,Cy,Cz,theta,psi,phi)

function ExceedEllipsoid(r1,r2,r3,Cx,Cy,Cz,theta,psi,phi)

% make ellipsoid 
N = 20 ;
[xe,ye,ze] = ellipsoid(0,0,0,r1,r2,r3,20); 

% define rotation matrix 
Dr=[cos(psi)*cos(phi)-cos(theta)*sin(psi)*sin(phi)  sin(psi)*cos(phi)+cos(theta)*cos(psi)*sin(phi) sin(theta)*sin(phi)  ...
    ; -cos(psi)*sin(phi)-cos(theta)*sin(psi)*cos(phi) -sin(psi)*sin(phi)+cos(theta)*cos(psi)*cos(phi) sin(theta)*cos(phi)  ...
    ;                  sin(theta)*sin(psi)                         -sin(theta)*cos(psi)                         cos(theta)     ];

% rotate the ellipsoid by Dr 
for j=1:1:(N+1)
    for k=1:1:(N+1)       
        V=Dr'*[xe(j,k) ; ye(j,k) ; ze(j,k)];      
        xe(j,k)=V(1);
        ye(j,k)=V(2);
        ze(j,k)=V(3);        
    end
end

% translate coordinates of ellipsoid 
xe=xe+Cx;
ye=ye+Cy;
ze=ze+Cz;

% obtain patch object of ellipse 
myellipsoid_patch = surf2patch(xe,ye,ze,ze);
vertices = myellipsoid_patch.vertices;
faces = myellipsoid_patch.faces;
faces(:,4) = []; % remove the fourth column of faces.  

% make meshgrid of points which are separated by 1-unit in x-, y- and z-axes 
[xm,ym,zm] = ndgrid(floor(min(xe(:))):1:ceil(max(xe(:))), ...
                    floor(min(ye(:))):1:ceil(max(ye(:))),...
                    floor(min(ze(:))):1:ceil(max(ze(:)))) ;
testpoints = [xm(:) ym(:) zm(:)];

figure;
% uncomment to make scatter plot of surface of ellipsoid 
% scatter3(xe(:),ye(:),ze(:),...
%         'MarkerEdgeColor','k',...
%          'MarkerFaceColor','b') % 
% hold on

% uncomment to see the meshgrid of testpoints 
% scatter3(testp(:,1),testp(:,2),testp(:,3),'k.');
% hold on 

% using intriangulation to find testpoints within myellipsoid 
heavytest = 5;
in = intriangulation(vertices,faces,testpoints,heavytest);

% plotting the surfaces and the vertices 
h = trisurf(faces,vertices(:,1),vertices(:,2),vertices(:,3));
set(h,'FaceColor','black','FaceAlpha',1/3,'EdgeColor','none');
hold on;

% plotting the testpoints determined to be inside the ellipsoid 
plot3(testpoints(in==1,1),testpoints(in==1,2),testpoints(in==1,3),'ro');

end

[Rik]

I think your question is becoming a bit too big te be answered in a comment. Consider posting it as a separate question.