MATLAB Answers

Determining internal coordinates of a rotated ellipsoid

조회 수: 32(최근 30일)
wim
wim 16 Apr 2021 13:17
편집: wim 26 Apr 2021 14:32
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

채택된 답변

Matt J
Matt J 16 Apr 2021 13:33
편집: Matt J 16 Apr 2021 13:33
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.
  댓글 수: 4
Matt J
Matt J 19 Apr 2021 12:06
You're welcome, but please Accept-click the answer to indicate that your question has been resolved.

댓글을 달려면 로그인하십시오.

추가 답변(0개)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by