Determining internal coordinates of a rotated ellipsoid

I would like to find all the internal points of a rotated ellipsoid.

• 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 (FEX)

I am unsure why intriangulation generates points outside my ellipsoid (see attached picture). I have attached my function and running script below. I have increased the heavytest parameter to 20, but this doesn't improve the result.

Any suggestions / recommendations are welcome!

% 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(testpoints(:,1),testpoints(:,2),testpoints(:,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

Answer 1

[Matt J]

Why is triangulation part of the strategy if you want actual points inside the ellipsoid? If you have the means to sample the surface of one ellipsoid, you could just create multiple concentric/coaxial smaller ellipsoids to sample the interior.

[wim]

I used triangulation because it was suggested by a Rik in a previous, similar post. That question however did not deal with rotated ellipsoids.

[Matt J]

OK, but I've suggested a viable alternative?

[wim]

Thanks for the suggestion Matt. The function below plots all internal points.

Edit: More internal points can be obtained by increasing N

function internalPointsEllipsoid(r1,r2,r3,Cx,Cy,Cz,theta,psi,phi)
% r1 r2 r3 double length of semi-axes 
% Cx Cy Cz center of ellipsoid 
% theta psi phi euler angles of rotation

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

% 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;
sPoints = [xe(:) ye(:) ze(:)];

%% create smaller ellipsoids to sample the interior of the ellipsoid 
stepsize = 0.1; 
tol = 0.05; % tolerance value 
interiorPoints = struct(); % store points 
count = 0; 
while 1
    r1 = r1-stepsize; 
    r2 = r2-stepsize; 
    r3 = r3-stepsize; 
    if r1 < tol || r2 < tol || r3 < tol
        break 
    else
        count = count + 1;
    end 
    % make ellipsoid 
    [xe,ye,ze] = ellipsoid(0,0,0,r1,r2,r3,N);
    % 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;
    
    % save coordinates in the struct file 
    interiorPoints(count).coords = [xe(:) ye(:) ze(:)];
end 

%% plotting the points 
figure; 
scatter3(sPoints(:,1),sPoints(:,2),sPoints(:,3),...
        'MarkerEdgeColor','k',...
         'MarkerFaceColor','b') 
hold on 
n = length(interiorPoints);

for i = 1:n
    coords = interiorPoints(i).coords;
    scatter3(coords(:,1),coords(:,2),coords(:,3),'.', ...
        'MarkerEdgeColor','r')
end 
hold off
end 

[Matt J]

You're welcome, but please Accept-click the answer to indicate that your question has been resolved.