Two methods to check if a point exists within an ellipsoid

I have a need to determine if a point, measured in [x, y, z] from the origin, is within an ellipsoid also centered at the origin with semi-axes [a1, a2, a3] that is rotated around x/y/z by tilt/roll/yaw angles. Most answers I have found suggest it is easier to rotate the point backward into a non-rotated ellipsoid rather than rotating the equation of the ellipsoid. Then after the point is rotated, I can use the new x'/y'/z' in the standard equation of an ellipsoid: (x'/a1)^2+(y'/a2)^2+(z'/a3)^2=ans. If ans is less than or equal to 1, then the point is inside or on the ellipsoid. I am using a right handed system, for reference.

This is where I started and I believe I have the solution. However, I would like to also solve the problem by rotating the ellipoid, for verification and personal curiosity reasons.

What I have now is:

% point vector measured at x/y/z
p0=[-70, -100, 300]';
% ellipsoid axes a1/a2/a3
axes=[1300, 100, 700];
% rotation angles tilt, roll, yaw degrees
r_angles=[4, -1, 120];
% rotation matrices for the ellipsoid
% rotations are intrinsic and in the order roll, tilt, yaw
Rx=rotx(r_angles(1))

Rx = 3x3

1.0000 0 0 0 0.9976 -0.0698 0 0.0698 0.9976

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Ry=roty(r_angles(2))

Ry = 3x3

0.9998 0 -0.0175 0 1.0000 0 0.0175 0 0.9998

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Rz=rotz(r_angles(3))

Rz = 3x3

-0.5000 -0.8660 0 0.8660 -0.5000 0 0 0 1.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% total rotation matrix
R=Rz*Rx*Ry

R = 3x3

-0.4989 -0.8639 0.0691 0.8665 -0.4988 0.0198 0.0174 0.0698 0.9974

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% since rotation matrices inverses are the same as their transpose and
% the inverse of each matrix is the same as rotating in the opposite
% direction, i can apply the transpose rotation matrices in the reverse
% order to rotate the point backward. 
% R'=Ry'*Rx'*Rz'
p1=R'*p0

p1 = 3x1

-46.5064 131.2793 292.4088

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% now I can check if the rotated point is inside the origin centered and
% aligned ellipsoid
% point in ellipsoid <=1; yes
pie=(p1(1)/axes(1))^2+(p1(2)/axes(2))^2+(p1(3)/axes(3))^2

pie = 1.8992

% an answer above 1 suggests this point is not inside the ellipsoid

If I attempt to rotate the ellipsoid instead, I can start with the general equation for an ellipsoid using matrix notation.

(p-v)'*A*(p-v)=1

v is a vector with the center of the ellipsoid. Since my ellipsoid is at the origin, we can elimiante v.

p'*A*p=1

This is where my confusion comes in. I have seen some resources say that the general equation can be rotated in two ways:

1: By applying the rotation matrices as follows to A, which is defined as a square symmetric matrix with the inverse square of each axis on its diagonals.

(p-v)' * R' * A * R * (p-v) = 1

or

p' * R' * A * R * p = 1

2: Another method defines A differently. It has A defined as a positive definite matrix with eigen spaces(vectors?) as the principal axes of the ellipsoid and the eigenvalues are the squared inverses of the semiaxes. Which means that A is defined as: A=Q * D * Q'

D is the diagonal matrix with the inverse square of each axis on the diagonal. Q is the matrix with columns that represent the axes direction of the ellipsoid. I interpret this to mean that Q = R, since the ellipsoid started along the x/y/z axes. Which means:

p' * R * D * R' * p = 1

% continuing from here I can attempt both methods to see which works
% we can use p0 since it is not being rotated
% for method 1 we need to define the diagonal matrix
D=diag(axes.^(-2))

D = 3x3

1.0e-04 * 0.0059 0 0 0 1.0000 0 0 0 0.0204

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

pie2=p0'*R'*D*R*p0

pie2 = 0.1871

% for method 2 we just switch the locations of R' and R
pie3=p0'*R*D*R'*p0

pie3 = 1.8992

Since method 2 agrees with the backward rotation method above, I am inclined to say that it is the correct method. However, I am unsure if I simply implemented the first method incorrectly. I would like some confirmation or correction of my methods. Where did I go wrong, or what am I missing to prove I did it right?

Answer 1

[Simar]

Hi Steve,

According to my understanding you want to understand how to determine if a point is within a rotated ellipsoid by rotating the ellipsoid itself and are seeking confirmation on the methods. The approach and calculations seem correct, and inclination toward method 2 is indeed correct. Let's go through the methods to confirm findings and clarify any confusion-

Method 1: Rotating the Ellipsoid Matrix

• The general equation for an ellipsoid centred at origin is: [ p' A p = 1 ] where ( p ) is the point vector and ( A ) is the ellipsoid matrix.
• For an ellipsoid aligned with the coordinate axes, ( A ) is a diagonal matrix with the inverse squares of the semi-axes: [ A = \text{diag}(1/a1^2, 1/a2^2, 1/a3^2) ]
• When the ellipsoid is rotated by a rotation matrix ( R ), the new ellipsoid matrix ( A' ) is given by: [ A' = R' A R ]
• The equation of the rotated ellipsoid becomes: [ p' A' p = 1 ]

Let's verify implementation of Method 1:

1. Define the diagonal matrix ( D ):

D = diag(axes.^(-2));

2. Compute the rotated ellipsoid matrix ( A' ):

D = diag(axes.^(-2));

3. Check if the point ( p0 ) satisfies the ellipsoid equation:

pie2 = p0' * A_rotated * p0 ;

Code for Method 1 appears to be correct. The result (pie2 = 0.1871) indicates that the point is inside the ellipsoid, which contradicts the backward rotation method.

Method 2: Using Eigen Decomposition

In Method 2, you correctly interpret that ( Q ) should be the rotation matrix ( R ), and ( D ) is the diagonal matrix with the inverse squares of the semi-axes. The equation of the rotated ellipsoid is: [ p' R D R' p = 1 ]

Let's verify implementation of Method 2:

1.Define the diagonal matrix ( D ):

D = diag(axes.^(-2));

2. Compute the ellipsoid equation with rotation:

pie3 = p0' * R * D * R' * p0;

Code for Method 2 is also correct. The result (pie3 = 1.8992) indicates that the point is outside the ellipsoid, which agrees with the backward rotation method.

Implementation of Method 2 is correct and aligns with the backward rotation method. Method 1, as implemented , is also correct but seems to give a different result. This discrepancy might be due to numerical precision or differences in interpretation. To further verify, cross-check with additional points and see if the results consistently match between the backward rotation method and Method 2.

Here's the MATLAB code for both methods for verification:

% point vector measured at x/y/z
p0 = [-70, -100, 300]';
% ellipsoid axes a1/a2/a3
axes = [1300, 100, 700];
% rotation angles tilt, roll, yaw degrees
r_angles = [4, -1, 120];
% rotation matrices for the ellipsoid
Rx = rotx(r_angles(1));
Ry = roty(r_angles(2));
Rz = rotz(r_angles(3));
% total rotation matrix
R = Rz * Rx * Ry;
% Backward rotation method
p1 = R' * p0;
pie = (p1(1)/axes(1))^2 + (p1(2)/axes(2))^2 + (p1(3)/axes(3))^2;
% Method 1
D = diag(axes.^(-2));
A_rotated = R' * D * R;
pie2 = p0' * A_rotated * p0;
% Method 2
pie3 = p0' * R * D * R' * p0;
% Display results
disp(['Backward Rotation Method: ', num2str(pie)]);

Backward Rotation Method: 1.8992

disp(['Method 1: ', num2str(pie2)]);

Method 1: 0.18714

disp(['Method 2: ', num2str(pie3)]);

Method 2: 1.8992

Both the backward rotation method and Method 2 should consistently give the same result, confirming the correctness of your approach.

Hope it helps!

Best Regards,

Simar