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
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.
Comments
Post a Comment