How to transform general Ellipsoid into rotated axis-oriented one?

Question

Recently, I discovered a long-standing bug in the rendering of Disk and Circle primitives after applying GeometricTransformation with a matrix as the second argument. This bug affects rendering of general ellipsoids of the form Ellipsoid[p, Σ], but does not affect axes-oriented ellipsoids of the form Ellipsoid[p, {r1, …}], because the latter are directly translated into the corresponding Disk objects with the same syntax and arguments. Applying Rotate (or GeometricTransformation with RotationTransform as the second argument) to the axis-oriented ellipsoids is a workaround for the bug. But for this we should be able to obtain from the matrix Σ the corresponding semiaxes lengths {r1, …} and the rotation angle Θ (for the 3D case, we also need the rotation axis w).

There is an elegant and efficient built-in way to perform the opposite task. For the 2D case, we have:

TransformedRegion[Ellipsoid[{x, y}, {r1, r2}], RotationTransform[Θ, {x, y}]]

Ellipsoid[{x, y}, 
   {{r1^2 Cos[Θ]^2 + r2^2 Sin[Θ]^2, r1^2 Cos[Θ] Sin[Θ] - r2^2 Cos[Θ] Sin[Θ]}, 
    {r1^2 Cos[Θ] Sin[Θ]] - r2^2 Cos[Θ] Sin[Θ], r2^2 Cos[Θ]^2 + r1^2 Sin[Θ]^2}}]

But I failed to find a way to transform Ellipsoid[p, Σ] into Rotate[Ellipsoid[p, {r1, …}], Θ, p] (or another suitable syntax form of Rotate or RotationTransform). Is it possible to do this in an elegant and efficient way (2D and 3D cases)?

Answer 1

We can employ Eigensystem in order to find the principal axes. While the rotation angle can be easily found in 2D case, the 3D case involves also the detection of the rotation axis; this is done by utilizing NullSpace.

2D case

SeedRandom[666];
Σ = RandomReal[{-1, 1}, {2, 2}];
Σ = Transpose[Σ].Σ;
p = RandomReal[{-1, 1}, {2}]; 

{λ, U} = Eigensystem[Σ];
r = Sqrt[λ];
rot = Det[U] Transpose[U];
θ = ArcTan @@ (rot[[All, 1]]);
Max[Abs[RotationMatrix[θ] - rot]]
Graphics[{
  FaceForm[ColorData[97][2]], Ellipsoid[p, Σ],
  EdgeForm[{Thickness[0.015], ColorData[97][1]}], 
  Rotate[#, θ, p] &@Ellipsoid[p, r]
  }]

0.

3D case

SeedRandom[666];
Σ = RandomReal[{-1, 1}, {3, 3}];
Σ = Transpose[Σ].Σ;
p = RandomReal[{-1, 1}, {3}];    

{λ, U} = Eigensystem[Σ];
r = Sqrt[λ];
rot = Det[U] Transpose[U];
w = NullSpace[rot - IdentityMatrix[3]][[1]];
V = RotationMatrix[{{1, 0, 0}, w}];
θ = -ArcTan @@ (Transpose[V].rot.V)[[2, 2 ;; 3]];
Max[Abs[RotationMatrix[θ, w] - rot]]
Show[
 Graphics3D[{
   FaceForm[ColorData[97][2]], 
   Ellipsoid[p, Σ],
   ColorData[97][1],
   Rotate[#, θ, w, p] &@Ellipsoid[p, r]
   }, Lighting -> "Neutral"]
 ]

7.77156*10^-16

Since RotationMatrix is the bottleneck here, I also link to this post of mine.