1
$\begingroup$

I have a list of points of a convex polyhedron for which I need to compute the eigenvectors of the Moment of Inertia matrix. But, it seems like the built-in MomentOfInertia behave very differently for ConvexHullMesh and ConvexHullRegion

pts3 = {{0.9000000000000000222`8., 0.4500000000000000111`8., 0.5`8.}, {0.6750000000000000444`8., 0.8397114317029974462`8., 0.5`8.}, {0.2250000000000000056`8., 0.8397114317029974462`8., 0.5`8.}, {0, 0.4500000000000000111`8., 0.5`8.}, {0.2250000000000000056`8., 0.060288568297002243`8., 0.5`8.}, {0.6750000000000000444`8., 0.060288568297002243`8., 0.5`8.}, {0.9000000000000000222`8., 0.4500000000000000111`8., 0.1`8.}, {0.6750000000000000444`8., 0.060288568297002243`8., 0.1`8.}, {0.2250000000000000056`8., 0.060288568297002243`8., 0.1`8.}, {0, 0.4500000000000000111`8.,0.1`8.}, {0.2250000000000000056`8., 0.8397114317029974462`8., 0.1`8.}, {0.6750000000000000444`8., 0.8397114317029974462`8., 0.1`8.}};

(*plotting function*)
plt[region_, eigvec_, cent_] := Show[region,Graphics3D[{{Green, Thick, InfiniteLine[cent, #] &/@ eigvec}, 
{Opacity[0.4], Red, Hyperplane[#, cent] & /@ eigvec}}, ImageSize -> Small]]


(*using convexhullmesh *)
reg1 = ConvexHullMesh[pts3];
cent1 = RegionCentroid@reg1;
{eigval1, eigvec1} = Eigensystem@MomentOfInertia@reg1;

plt[reg1, eigvec1, cent1]

enter image description here

(*using convexhullregion*)
reg2 = ConvexHullRegion[pts3];
cent2 = RegionCentroid@reg2;
{eigval2, eigvec2} = Eigensystem@MomentOfInertia@reg2;

plt[Region@reg2, eigvec2, cent2]

enter image description here

Observe the orientation of the cutting planes; they look very different for the two cases. I think I can use DiscretizeRegion@ConvexHullMesh[pts3] to make the cutting planes in the first case to look like the latter. However, I am confused as to why MomentOfInertia is behaving differently in the two cases? Is this a bug?

Improve this question
$\endgroup$
11
  • $\begingroup$ MMA version 12.3.0 I get twice the same symmetric picture. $\endgroup$
    – Daniel Huber
    Commented Dec 7, 2021 at 14:43
  • $\begingroup$ @DanielHuber I am using Mathematica 12.3.1 for Microsoft Window (64-Bit) (June 19,2021) $\endgroup$
    – Ali Hashmi
    Commented Dec 7, 2021 at 14:44
  • 1
    $\begingroup$ @AliHashmi Yes ConvexHullRegion was introduced with v12.2. I only recognized that Show[reg2] doesn't work, I changed it to Show[Region[reg2]] : $\endgroup$
    – Ulrich Neumann
    Commented Dec 7, 2021 at 16:19
  • 2
    $\begingroup$ It's not an artefact I think: The body has more than 3 symmtry axes, concerning the inertia that means complete rotational symmetry. Every pair of eigenvectors might be arbitrarily rotated! $\endgroup$
    – Ulrich Neumann
    Commented Dec 7, 2021 at 21:23
  • 1
    $\begingroup$ @UlrichNeumann's answer seems like the correct one. But if you are looking for a more 'natural' looking eigen system, you can turn your almost regular extruded hexagon into a regular and exact hexagon through RootApproximant[N[pts3]]. $\endgroup$
    – Greg Hurst
    Commented Dec 8, 2021 at 12:48

2 Answers 2

Reset to default
4
$\begingroup$

It's not an artefact, Mathematica result is ok for both cases!

The hexagon plane (perpendicular to {0,0,1} ) has six axes of symmetry. That means, concerning the inertia , the body is rotationally symmetrical!

J1 = MomentOfInertia@reg1;
{eigval1, eigvec1} = Eigensystem@J1;

J2 = MomentOfInertia@reg2;
{eigval2, eigvec2} = Eigensystem@J2;

The second and third eigenvalues are identical, which confirms mentioned rotational symmetry of the two regions

That's why rotation of the eigenvectors around {0,0,1} gives a new set of equivalent eigenvectors

rot1 = RotationMatrix[\[CurlyPhi], eigvec1[[1]]] // Chop; 
ev1rot = eigvec1 . rm // Chop; (* new eigenvectors*)
(J1 . Transpose[#] - Transpose[#] . DiagonalMatrix[eigval1]) &[ev1rot] // Chop; 
(*{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}*)

Same is true for reg2.

That's why the different appearing eigenvectors eigvec1, eigvec2 fullfill the eigensystem of both regions:

(J2 . Transpose[#] - Transpose[#] . DiagonalMatrix[eigval1]) &[eigvec1] // Chop
(*{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}*)

(J1 . Transpose[#] - Transpose[#] . DiagonalMatrix[eigval1]) &[eigvec2] // Chop
(*{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}*)
Improve this answer
$\endgroup$
2
  • $\begingroup$ I understand your point. However, should not the eigensystem output the same result for convexhullmesh and convexhullregion, unless of course the two represent very different underlying objects? $\endgroup$
    – Ali Hashmi
    Commented Dec 8, 2021 at 12:50
  • $\begingroup$ I don't think so. convexhullmeshand convexhullregion don't know from each other. Both return a correct eigensystem, unfortunately the eigenvectors are different in this case $\endgroup$
    – Ulrich Neumann
    Commented Dec 8, 2021 at 13:22
2
$\begingroup$

A workaround is to threshold values near zero in the inertia matrix.

principleVectors[reg_] := Eigensystem[Threshold[MomentOfInertia[reg]]]

Max[Abs[principleVectors[reg1] - principleVectors[reg2]]]

3.81639*10^-17
Improve this answer
$\endgroup$
2
  • $\begingroup$ thanks. makes sense now. I think it is the small values of the off-diagonal entries that is making Eigensystem return different values. Btw, if you use reg1 = DiscretizeRegion@ConvexHullMesh[pts3] and reg2 = DiscretizeRegion@BoundaryMeshRegion@ConvexHullRegion[pts3] the off-diagonal entries of the matrix are significantly large and therefore Eigensystem returns approx the same eigenvectors $\endgroup$
    – Ali Hashmi
    Commented Dec 8, 2021 at 9:16
  • $\begingroup$ upvoting your answer ! $\endgroup$
    – Ali Hashmi
    Commented Dec 8, 2021 at 9:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.