[REFACTOR] Generalize dual contouring triangulation for arbitrary grid shapes and remove legacy code

gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py

@@ -62,7 +62,7 @@ def dual_contouring_multi_scalar(
         return all_meshes
 
     # * 1) Triangulation code
-    left_right_codes = get_triangulation_codes(octree_list, options)
+    left_right_codes = get_triangulation_codes(octree_list)
 
     # * 2) Dual contouring mask
     # ? I guess this mask is different that erosion mask

gempy_engine/core/data/options/evaluation_options.py

@@ -1,6 +1,9 @@
 import enum
 
 from dataclasses import dataclass
+from typing import Annotated
+
+from typing_extensions import deprecated
 
 
 class MeshExtractionMaskingOptions(enum.Enum):
@@ -20,7 +23,7 @@ class EvaluationOptions:
 
     mesh_extraction: bool = True
     mesh_extraction_masking_options: MeshExtractionMaskingOptions = MeshExtractionMaskingOptions.INTERSECT
-    mesh_extraction_fancy: bool = True
+    mesh_extraction_fancy: Annotated[bool, deprecated("Old extraction method not in use anymore")] = True
 
     evaluation_chunk_size: int = 500_000
 

gempy_engine/modules/dual_contouring/_parallel_triangulation.py

@@ -5,7 +5,6 @@
 from gempy_engine.config import AvailableBackends
 from ...core.backend_tensor import BackendTensor
 from ...core.data.dual_contouring_data import DualContouringData
-from ._triangulate import triangulate_dual_contouring
 from ...modules.dual_contouring.fancy_triangulation import triangulate
 
 # Multiprocessing imports

gempy_engine/modules/dual_contouring/_sequential_triangulation.py

@@ -1,14 +1,9 @@
 from typing import Any
 
-import numpy as np
-import warnings
-
 from ._aux import _surface_slicer
-from ...config import AvailableBackends
 from ._gen_vertices import _compute_vertices
 from ...core.backend_tensor import BackendTensor
 from ...core.data.dual_contouring_data import DualContouringData
-from ._triangulate import triangulate_dual_contouring
 from ...modules.dual_contouring.fancy_triangulation import triangulate
 
 
@@ -20,7 +15,7 @@ def _sequential_triangulation(dc_data_per_stack: DualContouringData,
                               ) -> tuple[DualContouringData, Any, Any]:
     """Orchestrator function that combines vertex computation and triangulation."""
     dc_data_per_surface, vertices = _compute_vertices(
-        dc_data_per_stack=dc_data_per_stack, 
+        dc_data_per_stack=dc_data_per_stack,
         debug=debug,
         surface_i=i,
         valid_edges_per_surface=valid_edges_per_surface
@@ -35,72 +30,37 @@ def _sequential_triangulation(dc_data_per_stack: DualContouringData,
     edges_normals[valid_edges] = dc_data_per_stack.gradients[slice_object]
 
     indices_numpy = _compute_triangulation(
-        dc_data_per_surface=dc_data_per_surface, 
+        dc_data_per_surface=dc_data_per_surface,
         left_right_codes=left_right_codes,
         edges_normals=edges_normals,
-        vertex= vertices
+        vertex=vertices
     )
 
-
     vertices_numpy = BackendTensor.t.to_numpy(vertices)
     return dc_data_per_surface, indices_numpy, vertices_numpy
 
 
-def _compute_triangulation(dc_data_per_surface: DualContouringData, 
+def _compute_triangulation(dc_data_per_surface: DualContouringData,
                            left_right_codes, edges_normals, vertex) -> Any:
     """Compute triangulation indices for a specific surface."""
 
-    if left_right_codes is None:
-        # * Legacy triangulation
-        indices = triangulate_dual_contouring(dc_data_per_surface)
-    else:
-        # * Fancy triangulation 👗
-        if BackendTensor.engine_backend != AvailableBackends.PYTORCH:
-            with warnings.catch_warnings():
-                warnings.simplefilter("ignore", category=RuntimeWarning)
-                voxel_normal = np.nanmean(edges_normals, axis=1)
-                voxel_normal = voxel_normal[(~np.isnan(voxel_normal).any(axis=1))]  # drop nans
-        else:
-            # Assuming edges_normals is a PyTorch tensor
-            nan_mask = BackendTensor.t.isnan(edges_normals)
-            valid_count = (~nan_mask).sum(dim=1)
-
-            # Replace NaNs with 0 for sum calculation
-            safe_normals = edges_normals.clone()
-            safe_normals[nan_mask] = 0
-
-            # Compute the sum of non-NaN elements
-            safe_normals[:,:,0]
-            sum_normals = BackendTensor.t.sum(safe_normals, 1)
-
-            # Calculate the mean, avoiding division by zero
-            voxel_normal = sum_normals / valid_count.clamp(min=1)
+    # * Fancy triangulation 👗
+    valid_voxels = dc_data_per_surface.valid_voxels
 
-            # Remove rows where all elements were NaN (and hence valid_count is 0)
-            voxel_normal = voxel_normal[valid_count > 0].reshape(-1, 3)
+    left_right_per_surface = left_right_codes[valid_voxels]
+    valid_voxels_per_surface = dc_data_per_surface.valid_edges[valid_voxels]
+    tree_depth_per_surface = dc_data_per_surface.tree_depth
 
-        valid_voxels = dc_data_per_surface.valid_voxels
+    voxels_normals = edges_normals[valid_voxels]
 
-        left_right_per_surface = left_right_codes[valid_voxels]
-        valid_voxels_per_surface = dc_data_per_surface.valid_edges[valid_voxels]
-        voxel_normal_per_surface = voxel_normal
-        tree_depth_per_surface = dc_data_per_surface.tree_depth
-
-        voxels_normals = edges_normals[valid_voxels]
-        # nan_mask = BackendTensor.t.isnan(voxels_normals)
-        # voxels_normals[nan_mask] = 0
-
-        indices = triangulate(
-            left_right_array=left_right_per_surface,
-            valid_edges=valid_voxels_per_surface,
-            tree_depth=tree_depth_per_surface,
-            voxel_normals=voxels_normals,
-            vertex=vertex
-        )
-        # indices = BackendTensor.t.concatenate(indices, axis=0)
+    indices = triangulate(
+        left_right_array=left_right_per_surface,
+        valid_edges=valid_voxels_per_surface,
+        tree_depth=tree_depth_per_surface,
+        voxel_normals=voxels_normals,
+        vertex=vertex
+    )
 
     # @on
     indices_numpy = BackendTensor.t.to_numpy(indices)
     return indices_numpy
-
-

gempy_engine/modules/dual_contouring/dual_contouring_interface.py

@@ -83,32 +83,17 @@ def find_intersection_on_edge(_xyz_corners, scalar_field_on_corners,
 # endregion
 
 # region Triangulation Codes
-def get_triangulation_codes(octree_list: List[OctreeLevel], options: InterpolationOptions) -> np.ndarray | None:
+def get_triangulation_codes(octree_list: List[OctreeLevel]) -> np.ndarray | None:
     """
     Determine the appropriate triangulation codes based on options and octree structure.
 
     Args:
         octree_list: List of octree levels
-        options: Interpolation options
 
     Returns:
         Left-right codes array if fancy triangulation is enabled and supported, None otherwise
     """
-    is_pure_octree = bool(np.all(octree_list[0].grid_centers.octree_grid_shape == 2))
-
-    match (options.evaluation_options.mesh_extraction_fancy, is_pure_octree):
-        case (True, True):
-            return get_left_right_array(octree_list)
-        case (True, False):
-            warnings.warn(
-                "Fancy triangulation only works with regular grid of resolution [2,2,2]. "
-                "Defaulting to regular triangulation"
-            )
-            return None
-        case (False, _):
-            return None
-        case _:
-            raise ValueError("Invalid combination of options")
+    return get_left_right_array(octree_list)
 
 
 def get_masked_codes(left_right_codes: np.ndarray, mask: np.ndarray) -> np.ndarray: