- 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 ]
Two methods to check if a point exists within an ellipsoid
37 次查看(过去 30 天)
显示 更早的评论
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))
Ry=roty(r_angles(2))
Rz=rotz(r_angles(3))
% total rotation matrix
R=Rz*Rx*Ry
% 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
% 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
% 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))
pie2=p0'*R'*D*R*p0
% for method 2 we just switch the locations of R' and R
pie3=p0'*R*D*R'*p0
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?
0 个评论
回答(1 个)
Simar
2024-8-1
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
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)]);
disp(['Method 1: ', num2str(pie2)]);
disp(['Method 2: ', num2str(pie3)]);
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
0 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!