A triangle-shaped mesh polygon merged from several hull triangles can get a reversed outline, so mesh_polynormal points inward #3701

Intro

Hi!

My setup

• MuJoCo 3.16.0, Python bindings (PyPI wheel). Also reproduced with the nightly wheel 3.16.1.dev1001381281 (includes the fix for Mesh polygon normal computed from three nearly collinear outline vertices; multiccd then returns a deep contact with a sideways normal #3691) and with a C build of main at 651bae1; all three give identical output.
• macOS 27.0.1, arm64 (Apple Silicon), Python 3.12.13.
• Default options.

What's happening? What did you expect?

When the compiler merges several hull triangles into one polygon and the boundary of that polygon happens to be a triangle, the polygon's outline (its boundary vertices in mesh_polyvert) can be stored in clockwise order as seen from outside. The polygon normal is computed from the outline, so mesh_polynormal then points into the mesh. All other polygons are stored counter-clockwise, with an outward normal, and I expected the same here.

Where the order comes from. MakePolygons groups the hull triangles by normal direction and inserts the triangles of each group into a MeshPolygon one at a time. MeshPolygon::InsertFace stores the polygon's boundary as a list of directed edges. Each edge keeps the direction it has in its own triangle, and all hull triangles are wound counter-clockwise seen from outside. When a new triangle shares edges with the boundary, each shared edge appears twice with opposite directions and is removed, and the new triangle's remaining edges are appended to the end of the list. So once all triangles are inserted, every boundary edge points counter-clockwise around the polygon (seen from outside), but the edges are listed in insertion order, not in the order in which they follow each other around the boundary.

MeshPolygon::Paths turns the edge list into an outline by starting at one edge and repeatedly taking the edge that begins where the current one ends. Because every edge already points counter-clockwise, this walk always produces the counter-clockwise order. Polygons with exactly three edges skip the walk. Instead, a shortcut returns the start vertex of each edge in list order. It was added in ae1c3b7 together with the rest of this code, with the comment "shortcut if polygon is just a triangular face". The shortcut is correct when the polygon is a single triangle, because then the three edges are stored in the triangle's own order. But for a triangle-shaped polygon made of several triangles, the list order depends on the insertion order and need not follow the boundary.

For example, take a polygon with corners a, b, c (counter-clockwise from outside) and one interior vertex r, made of the triangles (a, b, r), (b, c, r) and (c, a, r). Inserting them in that order leaves the edges (a, b), (b, c), (c, a), and the shortcut returns (a, b, c). Inserting them in the order (a, b, r), (c, a, r), (b, c, r) leaves (a, b), (c, a), (b, c), and the shortcut returns (a, c, b), which is clockwise. MakePolygonNormals computes the normal from the outline in stored order, so the second case gives an inward normal.

The insertion order follows the order of the hull's triangles, which in turn depends on the order of the input vertices. So for a given shape, whether such a polygon gets an inward normal depends only on the order in which the vertices are listed. In the model below, a triangular prism (a triangle with sides of about 1 m, extruded 0.1 m) has one extra vertex 0.01 mm below the centroid of its bottom face. The hull splits the bottom face into three triangles around that vertex, and MakePolygons merges them back into one triangular polygon, because the three triangles are nearly coplanar. Over 200 random orders of the same 7 vertices, this polygon's normal points inward in 59 and outward in 141.

Effect on contacts with a plane. For a mesh touching a plane, mjc_PlaneConvex first finds the deepest vertex of the mesh and adds it as a contact. It then looks at the polygons that have this vertex as a corner of their outline (listed in mesh_polymap), picks the one whose stored normal is most opposite to the plane normal, and adds that polygon's other corners that are within the contact margin of the plane as further contacts.

In the reproduction script below, the prism starts 0.5 mm into the floor, so that it is in contact from the first step. With either vertex order, the deepest vertex is at first the extra vertex below the bottom face. That vertex lies inside the merged polygon and is a corner of no polygon outline, so mjc_PlaneConvex makes only this one contact, under the prism's centroid. Up to this point the two vertex orders behave the same: the prism balances on this single contact, and small numerical errors tilt it by 0.002° within the first 1.6 s, until a bottom corner becomes the deepest vertex. From then on the stored normal decides what happens:

• With the outward normal, the bottom polygon is the one most opposite to the plane normal, so its three corners are added as contacts, and on the next step the prism is level again on the single centre contact. This happened twice in the first 4 s, and the printed samples show only the centre contact. The prism stays in place, 0.37 mm deep in the floor.
• With the inward normal, the bottom polygon's normal points the same way as the plane normal, so it is the worst candidate and a side face is picked instead. The only corners of a side face that are near the floor are the two ends of its bottom edge. Two contacts on one edge cannot hold the prism level, so it tips about that edge, another corner becomes the deepest vertex, and the next pair of contacts comes from another side face. The prism never comes back to rest: it sinks to about 6.5 mm into the floor and stays there, its contacts keep switching between pairs of bottom corners, its angular speed (sampled every 4 s) reaches 0.14 rad/s, and it slides 0.67 mm horizontally over 20 s.

#3691 fixed how the normal is computed from an outline; this issue is about the order of the outline itself, which is why the nightly wheel with that fix gives the same output.

Suggested fix, and why it is correct beyond this example

Remove the shortcut, so that three-edge polygons go through the same walk as every other polygon. I have a PR ready (6 deleted lines and one test). Four things make me confident that the change is right for all meshes and not only for this model:

1. No new logic. The walk already builds the outline of every polygon with four or more edges. For a three-edge polygon it starts at the first edge (u, v) and then takes the edge that begins at v, so it returns the boundary in counter-clockwise order whatever the list order is.
2. Only reversed outlines change. Up to rotation, three directed edges around a triangle can be listed in only two ways: in the order in which they follow each other around the boundary, or in the reverse of that order, such as (a, b), (c, a), (b, c) in the example above. For a single triangle the list is (v1, v2), (v2, v3), (v3, v1), and the walk returns v1, v2, v3, which is exactly what the shortcut returns; the same holds for any several-triangle polygon whose list is in boundary order. So the only outlines that can change are the reversed ones, and they become counter-clockwise.
3. Downstream users of the outline still work. Paths has one caller, MakePolygons, and its output reaches mesh_polyvert, mesh_polynormal and mesh_polymap. mesh_polymap only records which polygons have each vertex as a corner, which does not depend on the order. For negative mesh scale, the compiler afterwards reverses every outline to account for the mirroring (Negative mesh scale produces inward convex polygon normals #3604); that step expects counter-clockwise input, so it now also gives the right result for these polygons.
4. Measured on about a million polygons. I compiled main before and after the change and compared every polygon:

Across all random meshes there are 4,488 triangle-shaped polygons made of several triangles, and 2,229 of them were reversed before the change. In all three sets the vertices, the hull triangles and the number of polygons are identical before and after, and every changed outline is the same three vertices in the opposite order.

Apart from the new test, the full test suite gives the same results before and after the change, in double and in single precision. The new test fails before the change and passes after it. The 200 vertex orders of the prism all give an outward normal after the change, and with the inward-normal vertex order the prism then stays 0.37 mm deep in the floor without moving, like the other order.

Steps for reproduction

1. Save the script below as repro.py. The model is inlined.
2. Run python repro.py.
3. The bottom polygon's normal points inward for 59 of the 200 vertex orders. For the second vertex order printed, the prism sinks about 6.5 mm into the floor and keeps rocking; for the first one, it stays 0.37 mm deep and still.

Minimal model for reproduction

<mujoco>
  <asset>
    <mesh name="prism" vertex="0 0 0.1  1 0 0.1  0.5 0.3 -1e-05  0.5 0.9 0  0 0 0  0.5 0.9 0.1  1 0 0"/>
  </asset>
  <worldbody>
    <geom name="floor" type="plane" size="2 2 0.1"/>
    <body name="prism" pos="0 0 -0.0005">
      <freejoint/>
      <geom name="prism" type="mesh" mesh="prism"/>
    </body>
  </worldbody>
</mujoco>

Code required for reproduction

"""Triangular mesh polygon merged from several faces: reversed outline, inward mesh_polynormal.

Prism: a 1 m triangle extruded 0.1 m, plus ONE extra vertex 0.01 mm below the centre of the bottom face.
The hull splits the bottom face into three triangles; MakePolygons merges them back into one triangular polygon.
"""
import itertools
import random

import mujoco
import numpy as np

BOTTOM = [(0, 0, 0), (1, 0, 0), (0.5, 0.9, 0)]
BULGE = [(0.5, 0.3, -0.00001)]
TOP = [(0, 0, 0.1), (1, 0, 0.1), (0.5, 0.9, 0.1)]
VERTICES = BOTTOM + BULGE + TOP

XML = """
<mujoco>
  <asset>
    <mesh name="prism" vertex="{vertex}"/>
  </asset>
  <worldbody>
    <geom name="floor" type="plane" size="2 2 0.1"/>
    <body name="prism" pos="0 0 -0.0005">
      <freejoint/>
      <geom name="prism" type="mesh" mesh="prism"/>
    </body>
  </worldbody>
</mujoco>
"""


def load(order):
  vertex = "  ".join(" ".join(f"{c:g}" for c in VERTICES[i]) for i in order)
  return mujoco.MjModel.from_xml_string(XML.format(vertex=vertex))


def bottom_normal(model, data):
  """World normal of the triangular polygon at the bottom of the prism."""
  g = model.geom("prism").id
  R = data.geom_xmat[g].reshape(3, 3)
  vert = model.mesh_vert[model.mesh_vertadr[0]:][:model.mesh_vertnum[0]]
  for p in range(model.mesh_polyadr[0], model.mesh_polyadr[0] + model.mesh_polynum[0]):
    ids = model.mesh_polyvert[model.mesh_polyvertadr[p]:][:model.mesh_polyvertnum[p]]
    if len(ids) == 3 and (R @ vert[ids].mean(0) + data.geom_xpos[g])[2] < 0.05:
      return R @ model.mesh_polynormal[p]


print("MuJoCo", mujoco.__version__)

# 1. The outward normal of the bottom polygon is (0, 0, -1) whatever the vertex order.
random.seed(0)
orders = random.sample(list(itertools.permutations(range(7))), 200)
inward = []
for order in orders:
  model = load(order)
  data = mujoco.MjData(model)
  mujoco.mj_kinematics(model, data)
  if bottom_normal(model, data)[2] > 0:
    inward.append(order)
print(f"\n200 random vertex orders: bottom polygon normal points inward in {len(inward)}, "
      f"outward in {200 - len(inward)}")

# 2. The prism resting on the floor, for one order of each kind.
for label, order in (("vertex order as listed", tuple(range(7))), ("vertex order " + str(inward[0]), inward[0])):
  model = load(order)
  data = mujoco.MjData(model)
  mujoco.mj_forward(model, data)
  print(f"\n{label}: bottom polygon world normal {np.round(bottom_normal(model, data), 3) + 0.0}")
  start = data.qpos[:2].copy()
  for step in range(1, 10001):
    mujoco.mj_step(model, data)
    if step % 2000 == 0:
      points = " ".join(f"({c.pos[0]:.3f}, {c.pos[1]:.3f})" for c in data.contact[:data.ncon])
      depth = -1000 * min(c.dist for c in data.contact[:data.ncon])
      print(f"  t={step * model.opt.timestep:4.1f} s: penetration {depth:5.2f} mm, "
            f"drift {1000 * np.linalg.norm(data.qpos[:2] - start):4.2f} mm, "
            f"angular speed {np.linalg.norm(data.qvel[3:]):.4f} rad/s, contacts {points}")

MuJoCo 3.16.0

200 random vertex orders: bottom polygon normal points inward in 59, outward in 141

vertex order as listed: bottom polygon world normal [ 0.  0. -1.]
  t= 4.0 s: penetration  0.37 mm, drift 0.00 mm, angular speed 0.0000 rad/s, contacts (0.500, 0.300)
  t= 8.0 s: penetration  0.37 mm, drift 0.00 mm, angular speed 0.0000 rad/s, contacts (0.500, 0.300)
  t=12.0 s: penetration  0.37 mm, drift 0.00 mm, angular speed 0.0001 rad/s, contacts (0.500, 0.300)
  t=16.0 s: penetration  0.37 mm, drift 0.00 mm, angular speed 0.0000 rad/s, contacts (0.500, 0.300)
  t=20.0 s: penetration  0.37 mm, drift 0.00 mm, angular speed 0.0000 rad/s, contacts (0.500, 0.300)

vertex order (4, 5, 3, 2, 0, 6, 1): bottom polygon world normal [0. 0. 1.]
  t= 4.0 s: penetration  6.52 mm, drift 0.09 mm, angular speed 0.1360 rad/s, contacts (1.000, -0.000) (0.500, 0.900)
  t= 8.0 s: penetration  6.54 mm, drift 0.23 mm, angular speed 0.0215 rad/s, contacts (-0.000, -0.000) (0.500, 0.900)
  t=12.0 s: penetration  6.48 mm, drift 0.38 mm, angular speed 0.1152 rad/s, contacts (1.000, -0.000) (-0.000, -0.000)
  t=16.0 s: penetration  6.52 mm, drift 0.53 mm, angular speed 0.1360 rad/s, contacts (0.999, -0.000) (0.500, 0.900)
  t=20.0 s: penetration  6.54 mm, drift 0.67 mm, angular speed 0.0215 rad/s, contacts (-0.001, -0.000) (0.499, 0.900)

Confirmations

• I searched the latest documentation thoroughly before posting.
• I searched previous Issues and Discussions, I am certain this has not been raised before.