Given an observer looking to the center of an ellipsoid, how can I determine the distance through the ellipsoid if none of the axes are aligned with this look angle?

I have a problem where I need to find out what the distance through an ellipsoid is from a viewer's perspective (also need to know what the two perpendicular distances, vertical and horizontal, through the ellipsoid are as well).

I currently have a observer position, an ellipsoid center position, the semi-axes of the ellipsoid, and the angles the ellipsoid is rotated from a starting point of East = axis 1, North = axis 2, Up = axis 3, tilt, roll, and yaw respectiely.

I started by converting the geodetic coordinates of each party into ECEF and finding the unit vector of the look angle. Then once the ellipsoid is rotated the dot product of each axis unit vector is compared to the unit vector of the look angle to see if any align (dot product result is very clsoe to 1 or -1).

obs=[lat, lon, alt]
ell_c=[lat, lon, alt]
ell_axes=[ax1,ax2,ax3]
ell_ang=[tilt, roll, yaw]
[obsx, obsy, obsz]=geodetic2ecef(obs(1),obs(2),obs(3),wgs84Ellipsoid)
[ellx, elly, ellz]=geodetic2ecef(ell_c(1),ell(2)_c,ell(3)_c,wgs84Ellipsoid)
look_ang=[ellx; elly; ellz]-[obsx; obsy; obsz]
u_la=look_ang/norm(look_ang)
E=blkdiag(ell_axes(1),ell_axes(2),ell_axes(3))
rE=rotz(ell_ang(3))*rotx(ell_ang(1))*roty(ell_ang(2))*E % rotated ellipsoid, convention being used is to rotate in order roll, tilt, yaw
rE_u1=rE(:,1)/norm(rE(:,1)) % unit vector of the first axis after rotations
rE_u2=rE(:,2)/norm(rE(:,2)) % unit vector of the second axis after rotations
rE_u3=rE(:,3)/norm(rE(:,3)) % unit vector of the third axis after rotations
dots=[dot(u_la,rE_u1);dot(u_la,rE_u2),dot(u_la,rE_u3)] % dot products of the look angle against each axis

So this is where I am code wise. I can confirm that none of the axes line up with the observer angle. How do I take the next step and get the dimension through the ellipsoid along the look axis?

Also, please critque my method to get this far already. I will happily accept an error on my part as an explanation why none of the axes line up.

[Cris LaPierre]

Full disclosure, I ran your questions through AI and got the following:

You’re very close conceptually — the main reasons your snippet won’t give the “right” answer are a mix of:

1. Calling geodetic2ecef with the wrong argument order (correct is [X,Y,Z] = geodetic2ecef(spheroid,lat,lon,h)
2. A typo in the ellipsoid-center indexing (ell(2)_c,ell(3)_c should be ell_c(2),ell_c(3))
3. A bug building dots (you accidentally create a 2×2 array instead of 3×1): dots = [dot(u_la,rE_u1); dot(u_la,rE_u2); dot(u_la,rE_u3)];
4. Rotation conventions / degrees-vs-radians / frame assumptions (often the biggest hidden issue).

As for if this completely fixes the issue, I don't know as we don't have enough information to run your code. Consider sharing a runnable example, along with what the correct output should be.

[Steve K]

Thank you, but this isn't my actual code as that lives elsewhere. I am not super concerned with the typos since I have those worked out in the code. I wrote this from memory. I will try to port a working example for context.

Answer 1

[Cris LaPierre]

Again, using AI to draft code that computes the distance through the rotated ellipsoid along the observer → ellipsoid-center look direction (even when no axis is aligned with the look direction). It uses the code you shared, fixing the issues mentioned above, and expands upon it to compute the distance.

Key idea: treat the rotated ellipsoid as a quadratic surface in ECEF and solve the ray–ellipsoid intersection for the ray

Then the distance through the ellipsoid is because u_la is a unit direction vector.

obs=[lat, lon, alt];
ell_c=[lat, lon, alt];
ell_axes=[ax1,ax2,ax3];
ell_ang=[tilt, roll, yaw];

wgs84 = wgs84Ellipsoid("meter");

% geodetic2ecef signature: geodetic2ecef(ellipsoid, lat, lon, h)
[obsx, obsy, obsz] = geodetic2ecef(wgs84, obs(1),   obs(2),   obs(3));
[ellx, elly, ellz] = geodetic2ecef(wgs84, ell_c(1), ell_c(2), ell_c(3));

look_ang = [ellx; elly; ellz] - [obsx; obsy; obsz];
u_la     = look_ang / norm(look_ang);

E  = blkdiag(ell_axes(1), ell_axes(2), ell_axes(3));

% rotated ellipsoid axes (your convention/order as written)
rE = rotz(ell_ang(3)) * rotx(ell_ang(1)) * roty(ell_ang(2)) * E;

rE_u1 = rE(:,1) / norm(rE(:,1));
rE_u2 = rE(:,2) / norm(rE(:,2));
rE_u3 = rE(:,3) / norm(rE(:,3));

% dot products of the look angle against each rotated axis direction
dots = [dot(u_la, rE_u1);
        dot(u_la, rE_u2);
        dot(u_la, rE_u3)];

% ---- Distance through ellipsoid along the look ray ----
% Build the ellipsoid quadratic form in ECEF:
% (x - center)' * A * (x - center) = 1

R = rotz(ell_ang(3)) * rotx(ell_ang(1)) * roty(ell_ang(2));   % rotation only
D = diag(1 ./ (ell_axes(:).^2));                              % body-frame quadric
A = R * D * R.';                                              % world/ECEF quadric

obs_ecef = [obsx; obsy; obsz];
ell_ecef = [ellx; elly; ellz];

% Ray: x(t) = obs_ecef + t*u_la
p = obs_ecef - ell_ecef;

qa = u_la.' * A * u_la;
qb = 2 * (p.' * A * u_la);
qc = p.' * A * p - 1;

disc = qb^2 - 4*qa*qc;

if disc < 0
    dist_through = 0;
    t_enter = NaN;
    t_exit  = NaN;
    x_enter = [NaN; NaN; NaN];
    x_exit  = [NaN; NaN; NaN];
else
    sqrtDisc = sqrt(disc);
    t1 = (-qb - sqrtDisc) / (2*qa);
    t2 = (-qb + sqrtDisc) / (2*qa);

    t_enter = min(t1,t2);
    t_exit  = max(t1,t2);

    x_enter = obs_ecef + t_enter*u_la;
    x_exit  = obs_ecef + t_exit*u_la;

    dist_through = abs(t_exit - t_enter);  % meters (u_la is unit length)
end

[Steve K]

I'll give this a shot. Did the AI cite any sources where this method could be academically verified to be correct?

[Cris LaPierre]

No citations, though it did provide the equations it will be using. I assumed you already knew the approach to use and were just looking for help implementing it in MATLAB.

For this forum, we aim to help with the MATLAB implementation, but leave it to you to determine what the correct approach is.

Answer 2

[Matt J]

Given an observer looking to the center of an ellipsoid

After converting to ECEF coordinates,

 xc=[xc1;xc2;xc3]               ECEF of ellipsoid center
 xobs =[xobs1;xobs2;xobs3]      ECEF of observer point

obtain the equation for the ellipsoid as,

(x-xc)*Q*(x-xc)'=1;

then the ray-ellispsoid intersection length (from an observer through the center of the ellipsoid) would be,

 intersectionLength=2/sqrt((xobs-xc)*Q*(xobs-xc))* norm(xobs-xc);
                                                    

The mathematical derivation is attached.

[Matt J]

obtain the equation for the ellipsoid as, (x-xc)*Q*(x-xc)'=1

There are tools on the File Exchange for obtaining the ellipsoid equation in this form, assuming again you've converted everything to ECEF,

https://www.mathworks.com/matlabcentral/fileexchange/87584-object-oriented-tools-to-fit-plot-conics-and-quadrics?s_tid=prof_contriblnk

gtEllip=ellipsoidalFit.groundtruth([],xc,[ax1,ax2,ax3], [yaw,pitch,roll]);

R=gtEllipse.R';

Q=R'*diag(1./[ax1,ax2,ax3]).^2*R;