diff --git a/DeepSDFStruct/flexicubes/flexicubes.py b/DeepSDFStruct/flexicubes/flexicubes.py index bdd39ec9..83e1debf 100644 --- a/DeepSDFStruct/flexicubes/flexicubes.py +++ b/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 diff --git a/tests/test_tet_orientation.py b/tests/test_tet_orientation.py new file mode 100644 index 00000000..58658420 --- /dev/null +++ b/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")