CalcSpatialInertia fails on this Hello Robot obj

[RussTedrake]

What happened?

Following #20444, I've put a simple failure reproduction here: https://github.com/RobotLocomotion/drake/commit/2fe7f248091a93e20fbb471f8a27502ba9caee73

This attempts to compute the spatial inertia from the hello robot model in mujoco menagerie but fails with

Expected: CalcSpatialInertia(unit_scale_obj, kDensity) doesn't throw an exception.
  Actual: it throws std::runtime_error with description "Spatial inertia fails SpatialInertia::IsPhysicallyValid().
 mass = 2.1647959683444456e-05
 Center of mass = [-0.24465954725602354  -3.041901842037344e-06  0.02583678929053868]
 Inertia about point P, I_BP =
[1.90172e-08  -1.2568e-11  1.35882e-07]
[-1.2568e-11  1.30846e-06  1.45236e-12]
[1.35882e-07  1.45236e-12   1.2914e-06]
 Inertia about center of mass, I_BBcm =
[ 4.56633e-09   3.54305e-12  -9.59847e-10]
[ 3.54305e-12  -1.79939e-09  -2.49019e-13]
[-9.59847e-10  -2.49019e-13  -4.40846e-09]
 Principal moments of inertia about Bcm (center of mass) =
[-4.509963383944518e-09  -1.799395202505969e-09  4.667840740414944e-09]

Additionally the CalcSpatialInertia doc strings do not actually declare that they could throw an exception for this reason.

Version

No response

What operating system are you using?

Ubuntu 22.04

What installation option are you using?

compiled from source code using Bazel

Relevant log output

No response

[jwnimmer-tri]

The CalcSpatialInertia(mesh=, ...) API documents that:

The mesh must fully enclose a volume (no cracks, no open manifolds, etc.). All triangles must be "wound" such that their normals point outward (according to the right-hand rule based on vertex winding). If these requirements are not met, a value will be returned, but its value is meaningless.

I'm guessing that the obj file in question might not meet those conditions.

So there are at least three problems:

(1) Instead of returning a "meaningless" value like it promised, sometimes CalcSpatialInertia(mesh=...) throws.

(2) The docs for the related overload CalcSpatialInertia(shape=...) fail to repeat the caveats from the more specific overload that it delegates to.

(3) Our parser(s?) call CalcSpatialInertia(shape=...) without sufficient regard to whether or not the preconditions for mesh sensibility have been meet. Allowing the "meaningless" value to flow through the parse would be fine, but at least the parser needs to warn the user about that case when encountered -- so it can't be completely ignorant of the circumstance.

Or, if this particular obj file actually meets those requirements, then we have even more going wrong beyond those points.

On first impression, the ideal fix would be to update the contract of CalcSpatialInertia(...) to always produce a valid answer -- no more "meaningless". If we can't do that, then instead we need to provide an overload with a structured error channel, and wire the errors up to the parser. This could be as simple as making the return value optional<> and returning null when the result cannot be computed.

In the meantime, possibly we should hotfix the function to at least never throw.

[RussTedrake]

Thanks Jeremy. I completely agree, and would note that (1) checking properly that the mesh satisfies those conditions in the parser directly would be impractical without a geometry helper method. (2) mujoco parses these files ok (and must calculate some inertia). We could look at what they do.

I think propagating the error condition properly up through the parser (without a throw) would be better than nothing. Removing the apparently overly conservative requirements would of course be better. As a fallback, I could imagine CalcSpatialInertia could try to compute the inertia and if it is not physically valid, then it could fall back to a simpler approximation (like taking the volume of the convex hull, or even obb), so that it always returns something valid.

[SeanCurtis-TRI]

Several thoughts:

Model compliance

The model you've referenced does indeed have arbitrary open faces (they are the red faces you can see):

It looks like it's appropriate for visualization, but not for proximity or physics. I'm curious what mujoco actually does with meshes of that type. We could certainly compute a spatial inertia assuming the mass were strictly distributed on the surface of the triangles. In that case, we wouldn't be dependent on closed manifolds or consistent winding.

Calculating spatial inertia

The primary problem with computing spatial inertia from an obj is there are many ways to produce a "meaningless" but apparently otherwise valid spatial inertia. In this case, it just turns out the non-compliance is sufficient to push the resulting tensor out of "meaningless" and into "non-physical".

Today, we can easily detect the non-physical tensors (which is what happens here belatedly), but we can't generally detect meaningless tensors. So, any change to the API and the promises cannot eliminate the fact that if you invalidate assumptions about the obj, you may get something back meaningless. Still, this would be the "sometimes it throws".

Alternatively, we could change to a different computation mechanism that would be less sensitive to those issues (the aforementioned "mass on surface" approach). I don't have any intuition of how that might mess with the physics -- how different that is from an assumption of uniform density for the enclosed space.

[SeanCurtis-TRI]

P.S. Geometric analysis is still currently beyond Drake's capacity. But it's been on my wish list for years.

[jwnimmer-tri]

As a point of reference, one of the ways that SDFormat specifies handling auto inertia for non-convex meshes is to convert the mesh into a voxelized approximation, and sum the inertia of the voxels.

[RussTedrake]

I think we already support non-convex. It's not clear that would work for open faces.

[jwnimmer-tri]

Well, I was trying to talk about what to do with non-convex meshes that we currently mis-process (not non-convex meshes overall), but anyway that was getting off the point. Sean had a better question..

I'm curious what mujoco actually does ...

This is the only question that matters.

We are writing a parser for the mujoco file format. The only question that matters is what is specified to happen, by the file format documentation, when auto-inertia is being used on an open mesh. Our job is to implement that specification.

If the specification doesn't say, they we ask upstream to clarify it, and in the meantime we can guess what the spec would say by reading the mujoco parsing code and seeing what happens in their first-party parser.

My attempts to reason about this bottom-up from CalcSpatialInertia were missing the forest for the trees. Yes, the function is badly-shaped, and yes we're not calling it properly, but we need to work the problem top-down. What math does the specification command, then therefore what should our code look like.

[SeanCurtis-TRI]

To be clear - non-convexity is not an issue here. Non-convexity doesn't cause problems at all. Open meshes do.

Agreed - we need to find out what mujoco is doing. 100% on board with that.

[jwnimmer-tri]

I was treating the Convex shape as a special case because computing its uniform-density inertia should always be easy and correct, even if some triangles are missing or wound the wrong way or etc. One we take the hull, the inertia should be easy.

Probably a moot point for this thread though, since I don't think the drake:declare_convex marker is part of our mjcf parser in the first place, so there are never any Convex shapes anyway.

[SeanCurtis-TRI]

BTW Mujoco explicitly has a parameter called "shellinertia" indicating whether the mass is distributed on the surface. Its default value is false.

Also, there's a compiler option, exactmeshinertia that defaults to false, but when true, emphasizes that the inertia is exact for "any closed mesh geometry" (emphasis mine).

So, having done a quick pass through the documentation, I didn't see enough information vis a vis computing the spatial inertia from a mesh. I'd invite others to see if they can find it, but I suspect we'll have to turn to the code/developers to find out.

[RussTedrake]

The mujoco ComputeVolume() code for a mesh appears to be here: https://github.com/google-deepmind/mujoco/blob/df7ea3ed3350164d0f111c12870e46bc59439a96/src/user/user_mesh.cc#L1214 . I agree that they aren't doing anything fancy there.

[RussTedrake]

My current mental model is that mujoco is probably generating physically invalid inertias for those models, but then not doing the consistency check. But that should be verified.

[RussTedrake]

in #21403, i've implemented the fallback of using the volume of the OBB if the CalcSpatialInertia fails. (with a TODO reminding me that we could do better by taking the volume of the convex hull).