Fix tetrahedron orientation in FlexiCubes output

DeepSDFStruct/flexicubes/flexicubes.py

@@ -28,10 +28,14 @@
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 # See the License for the specific language governing permissions and
 # limitations under the License.
+import logging
+
 import torch
 import warnings
 from DeepSDFStruct.flexicubes.tables import *
 
+logger = logging.getLogger(__name__)
+
 __all__ = ["FlexiCubes"]
 
 
@@ -1002,4 +1006,50 @@ def _tetrahedralize(
 
         tets = torch.cat([tets_surface, tets_inside])
         vertices = torch.cat([vertices, inside_verts, inside_cubes_center])
+        # The surface (pyramid) and interior sub-procedures above emit tets
+        # with inconsistent winding, so a large fraction come out inverted
+        # (negative signed volume). FEA solvers require a positive signed
+        # volume / Jacobian, so normalize every tet to positive orientation
+        # and drop degenerate (zero-volume) elements.
+        tets = self._orient_tets(vertices, tets)
         return vertices, tets
+
+    @staticmethod
+    def _orient_tets(vertices, tets):
+        """Return ``tets`` with a consistent, strictly positive signed volume.
+
+        For a tetrahedron with vertices ``(v0, v1, v2, v3)`` the signed volume
+        is proportional to ``det([v1 - v0, v2 - v0, v3 - v0])``. Elements with
+        a negative determinant are inverted; swapping the last two vertices
+        flips the orientation so the signed volume becomes positive without
+        changing the element's geometry. The vertex indices are merely
+        reordered, so gradients to ``vertices`` are unaffected.
+
+        Degenerate elements (four coplanar vertices, zero volume) cannot be
+        repaired by reordering and are removed instead. Vertices are left
+        untouched, so indices of the remaining tets stay valid.
+        """
+        if tets.shape[0] == 0:
+            return tets
+        # Only integer index masks are derived here, so skip autograd
+        # recording even when ``vertices`` requires gradients.
+        with torch.no_grad():
+            v0 = vertices[tets[:, 0]]
+            v1 = vertices[tets[:, 1]]
+            v2 = vertices[tets[:, 2]]
+            v3 = vertices[tets[:, 3]]
+            e1, e2, e3 = v1 - v0, v2 - v0, v3 - v0
+            signed_vol = torch.einsum("ij,ij->i", e1, torch.linalg.cross(e2, e3, dim=1))
+            inverted = signed_vol < 0
+            # Coplanar tets only evaluate to exactly zero up to floating-point
+            # rounding (which depends on the association order of the triple
+            # product), so compare against a tolerance relative to the Hadamard
+            # bound |e1||e2||e3| of the determinant instead of zero itself.
+            scale = e1.norm(dim=1) * e2.norm(dim=1) * e3.norm(dim=1)
+            degenerate = signed_vol.abs() <= 1e-5 * scale
+        tets = tets.clone()
+        tets[inverted] = tets[inverted][:, [0, 1, 3, 2]]
+        if degenerate.any():
+            logger.info(f"removed {int(degenerate.sum())} elements with 0 volume")
+            tets = tets[~degenerate]
+        return tets

tests/test_tet_orientation.py

@@ -0,0 +1,88 @@
+"""Regression tests for tetrahedron orientation in FlexiCubes output.
+
+FlexiCubes' ``_tetrahedralize`` builds tets from two sub-procedures (surface
+pyramids and interior edges) whose vertex orderings do not share a consistent
+winding, so a large fraction of the elements used to come out inverted
+(negative signed volume). That breaks FEA solvers, which require a positive
+signed volume / Jacobian on every element. The interior sub-procedure can
+additionally emit degenerate (exactly coplanar, zero-volume) tets, which are
+removed during extraction. These tests pin both fixes: every tet returned by
+FlexiCubes must have a strictly positive signed volume.
+"""
+
+import torch
+
+from DeepSDFStruct.flexicubes.flexicubes import FlexiCubes
+from DeepSDFStruct.optimization import tet_signed_vol
+
+
+def _sphere_volume_mesh(res=16, radius=0.7, center=(0.0, 0.0, 0.0)):
+    fc = FlexiCubes(device="cpu")
+    x_nx3, cube_fx8 = fc.construct_voxel_grid(
+        res, bounds=[[-1.0, -1.0, -1.0], [1.0, 1.0, 1.0]]
+    )
+    c = torch.tensor(center, dtype=x_nx3.dtype)
+    s_n = torch.linalg.norm(x_nx3 - c, dim=1) - radius
+    verts, tets, _ = fc(x_nx3, s_n, cube_fx8, res, output_tetmesh=True)
+    return verts, tets
+
+
+def test_orient_tets_flips_inverted_element():
+    """A single inverted tet is flipped to positive signed volume."""
+    verts = torch.tensor(
+        [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]
+    )
+    # [0, 1, 2, 3] is positively oriented; [0, 1, 3, 2] is its inversion.
+    tets = torch.tensor([[0, 1, 2, 3], [0, 1, 3, 2]])
+
+    assert (
+        tet_signed_vol(verts, tets) < 0
+    ).any(), "fixture should contain an inversion"
+
+    oriented = FlexiCubes._orient_tets(verts, tets)
+    vols = tet_signed_vol(verts, oriented)
+    assert (vols > 0).all(), f"expected all positive volumes, got {vols.tolist()}"
+    # Orientation is fixed by reordering indices only: |volume| is preserved.
+    assert torch.allclose(vols.abs(), tet_signed_vol(verts, tets).abs())
+
+
+def test_orient_tets_handles_empty():
+    verts = torch.zeros((0, 3))
+    tets = torch.zeros((0, 4), dtype=torch.long)
+    assert FlexiCubes._orient_tets(verts, tets).shape == (0, 4)
+
+
+def test_orient_tets_removes_degenerate_elements():
+    """Zero-volume (coplanar) tets are dropped."""
+    verts = torch.tensor(
+        [
+            [0.0, 0.0, 0.0],
+            [1.0, 0.0, 0.0],
+            [0.0, 1.0, 0.0],
+            [0.0, 0.0, 1.0],
+            [1.0, 1.0, 0.0],  # coplanar with the first three (z = 0 plane)
+        ]
+    )
+    tets = torch.tensor([[0, 1, 2, 3], [0, 1, 2, 4]])
+
+    oriented = FlexiCubes._orient_tets(verts, tets)
+
+    assert oriented.shape == (1, 4)
+    assert (tet_signed_vol(verts, oriented) > 0).all()
+
+
+def test_flexicubes_volume_mesh_has_only_positive_tets():
+    """End-to-end: extracted volume mesh has strictly positive signed volumes."""
+    for center in [(0.0, 0.0, 0.0), (0.13, 0.07, 0.21)]:
+        verts, tets = _sphere_volume_mesh(center=center)
+        assert tets.shape[0] > 0
+        vols = tet_signed_vol(verts, tets)
+        n_bad = int((vols <= 0).sum())
+        assert n_bad == 0, f"{n_bad} non-positive tets for center {center}"
+
+
+if __name__ == "__main__":
+    test_orient_tets_flips_inverted_element()
+    test_orient_tets_handles_empty()
+    test_flexicubes_volume_mesh_has_only_positive_tets()
+    print("ok")