[REFACTOR] Move dual contouring code to modules directory and improve mesh consistency and water tight prot (#28)

# Refactoring Dual Contouring Implementation

This PR reorganizes the dual contouring code to improve maintainability and fix mesh generation issues:

- Moved `_dual_contouring.py` from API to modules directory to better reflect its role
- Refactored the dual contouring implementation into smaller, more focused files:
  - `_gen_vertices.py` for vertex generation
  - `_triangulate.py` for triangulation logic
  - `_sequential_triangulation.py` and `_parallel_triangulation.py` for processing surfaces
- Added functionality to detect and fix overlapping voxels across stacks
- Added `left_right` property to `DualContouringMesh` to track triangulation codes
- Improved vertex visualization in PyVista helper by increasing point size
- Extracted common functionality into reusable functions like `_surface_slicer`
- Moved mask generation and triangulation code functions to the interface module

These changes improve code organization and lay groundwork for fixing mesh generation issues at stack boundaries.

Author: Leguark

gempy_engine/API/dual_contouring/_experimental_water_tight_DC_1.py

@@ -2,7 +2,7 @@
 
 import numpy as np
 
-from gempy_engine.API.dual_contouring._dual_contouring import compute_dual_contouring
+from gempy_engine.modules.dual_contouring._dual_contouring import compute_dual_contouring
 from gempy_engine.API.dual_contouring._interpolate_on_edges import interpolate_on_edges_for_dual_contouring
 from gempy_engine.core.data.dual_contouring_data import DualContouringData
 from gempy_engine.core.data.dual_contouring_mesh import DualContouringMesh

gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py

@@ -10,6 +10,7 @@ class DualContouringMesh:
     vertices: np.ndarray
     edges: np.ndarray
     dc_data: Optional[DualContouringData] = None  # * In principle we need this just for testing
+    left_right: Optional[np.ndarray] = None 
 
     def __repr__(self):
         return f"DualContouringMesh({self.vertices.shape[0]} vertices, {self.edges.shape[0]} edges)"

gempy_engine/core/data/dual_contouring_mesh.py

@@ -6,33 +6,33 @@
 from ...core.data.dual_contouring_data import DualContouringData
 from ...core.data.dual_contouring_mesh import DualContouringMesh
 from ...core.utils import gempy_profiler_decorator
-from ...modules.dual_contouring._parallel_triangulation import _should_use_parallel_processing, _process_surface_batch, _init_worker
-from ...modules.dual_contouring._sequential_triangulation import _sequential_triangulation
+from ._parallel_triangulation import _should_use_parallel_processing, _process_surface_batch, _init_worker
+from ._sequential_triangulation import _sequential_triangulation
+from ._gen_vertices import _compute_vertices
 
 # Multiprocessing imports
 try:
     import torch.multiprocessing as mp
+
     MULTIPROCESSING_AVAILABLE = True
 except ImportError:
     import multiprocessing as mp
-    MULTIPROCESSING_AVAILABLE = False
-
-
 
+    MULTIPROCESSING_AVAILABLE = False
 
 
 @gempy_profiler_decorator
-def compute_dual_contouring(dc_data_per_stack: DualContouringData, left_right_codes=None, debug: bool = False) -> List[DualContouringMesh]:
+def compute_dual_contouring(dc_data_per_stack: DualContouringData,
+                            left_right_codes=None, debug: bool = False) -> List[DualContouringMesh]:
     valid_edges_per_surface = dc_data_per_stack.valid_edges.reshape((dc_data_per_stack.n_surfaces_to_export, -1, 12))
 
     # Check if we should use parallel processing
     use_parallel = _should_use_parallel_processing(dc_data_per_stack.n_surfaces_to_export, BackendTensor.engine_backend)
     parallel_results = None
-    
-    if use_parallel and False: # ! (Miguel Sep 25) I do not see a speedup
+
+    if use_parallel and False:  # ! (Miguel Sep 25) I do not see a speedup
         print(f"Using parallel processing for {dc_data_per_stack.n_surfaces_to_export} surfaces")
         parallel_results = _parallel_process_surfaces(dc_data_per_stack, left_right_codes, debug)
-        
 
     # Fall back to sequential processing
     print(f"Using sequential processing for {dc_data_per_stack.n_surfaces_to_export} surfaces")
@@ -41,24 +41,27 @@ def compute_dual_contouring(dc_data_per_stack: DualContouringData, left_right_co
     for i in range(dc_data_per_stack.n_surfaces_to_export):
         # @off
         if parallel_results is not None:
-            _, vertices_numpy = _sequential_triangulation(
-                dc_data_per_stack,
-                debug, 
-                i, 
-                left_right_codes, 
-                valid_edges_per_surface,
-                compute_indices=False 
+            _, vertices_numpy = _compute_vertices(
+                dc_data_per_stack=dc_data_per_stack,
+                debug=debug,
+                surface_i=i,
+                valid_edges_per_surface=valid_edges_per_surface
             )
             indices_numpy = parallel_results[i]
         else:
-            indices_numpy, vertices_numpy = _sequential_triangulation(
-                dc_data_per_stack,
-                debug, 
-                i, 
-                left_right_codes, 
-                valid_edges_per_surface,
-                compute_indices=True
+            dc_data_per_surface, indices_numpy, vertices_numpy = _sequential_triangulation(
+                dc_data_per_stack=dc_data_per_stack,
+                debug=debug,
+                i=i,
+                valid_edges_per_surface=valid_edges_per_surface,
+                left_right_codes=left_right_codes
             )
+            
+            # Handle None left_right_codes case
+            if left_right_codes is not None:
+                valid_left_right_codes = left_right_codes[dc_data_per_surface.valid_voxels]
+            else:
+                valid_left_right_codes = None
 
         if TRIMESH_LAST_PASS := True:
             vertices_numpy, indices_numpy = _last_pass(vertices_numpy, indices_numpy)
@@ -67,14 +70,13 @@ def compute_dual_contouring(dc_data_per_stack: DualContouringData, left_right_co
             DualContouringMesh(
                 vertices_numpy,
                 indices_numpy,
-                dc_data_per_stack
+                dc_data_per_stack,
+                left_right=valid_left_right_codes
             )
         )
     return stack_meshes
 
 
-
-
 def _parallel_process_surfaces(dc_data_per_stack, left_right_codes, debug, num_workers=None, chunk_size=2):
     """Process surfaces in parallel using multiprocessing."""
     if num_workers is None:
@@ -96,6 +98,9 @@ def _parallel_process_surfaces(dc_data_per_stack, left_right_codes, debug, num_w
     surface_indices = list(range(dc_data_per_stack.n_surfaces_to_export))
     chunks = [surface_indices[i:i + chunk_size] for i in range(0, len(surface_indices), chunk_size)]
 
+    # Handle None left_right_codes case - ensure we pass a serializable value
+    serializable_left_right_codes = left_right_codes
+
     try:
         # Use spawn context for better PyTorch compatibility
         ctx = mp.get_context("spawn") if MULTIPROCESSING_AVAILABLE else mp
@@ -106,7 +111,7 @@ def _parallel_process_surfaces(dc_data_per_stack, left_right_codes, debug, num_w
             for chunk in chunks:
                 result = pool.apply_async(
                     _process_surface_batch,
-                    (chunk, dc_data_dict, left_right_codes, debug)
+                    (chunk, dc_data_dict, serializable_left_right_codes, debug)
                 )
                 async_results.append(result)
 
@@ -122,6 +127,7 @@ def _parallel_process_surfaces(dc_data_per_stack, left_right_codes, debug, num_w
         print(f"Parallel processing failed: {e}. Falling back to sequential processing.")
         return None
 
+
 def _last_pass(vertices, indices):
     # Check if trimesh is available
     try:
@@ -130,4 +136,4 @@ def _last_pass(vertices, indices):
         mesh.fill_holes()
         return mesh.vertices, mesh.faces
     except ImportError:
-        return vertices, indices
+        return vertices, indices

gempy_engine/API/dual_contouring/_dual_contouring.py renamed to gempy_engine/modules/dual_contouring/_dual_contouring.py

@@ -0,0 +1,125 @@
+from typing import Any
+
+import numpy as np
+
+from ...config import AvailableBackends
+from ...core.backend_tensor import BackendTensor
+from ...core.data.dual_contouring_data import DualContouringData
+
+
+def _compute_vertices(dc_data_per_stack: DualContouringData,
+                      debug: bool,
+                      surface_i: int,
+                      valid_edges_per_surface) -> tuple[DualContouringData, Any]:
+    """Compute vertices for a specific surface."""
+    valid_edges: np.ndarray = valid_edges_per_surface[surface_i]
+    next_surface_edge_idx: int = valid_edges_per_surface[:surface_i + 1].sum()
+    if surface_i == 0:
+        last_surface_edge_idx = 0
+    else:
+        last_surface_edge_idx: int = valid_edges_per_surface[:surface_i].sum()
+    slice_object: slice = slice(last_surface_edge_idx, next_surface_edge_idx)
+
+    dc_data_per_surface = DualContouringData(
+        xyz_on_edge=dc_data_per_stack.xyz_on_edge,
+        valid_edges=valid_edges,
+        xyz_on_centers=dc_data_per_stack.xyz_on_centers,
+        dxdydz=dc_data_per_stack.dxdydz,
+        exported_fields_on_edges=dc_data_per_stack.exported_fields_on_edges,
+        n_surfaces_to_export=dc_data_per_stack.n_surfaces_to_export,
+        tree_depth=dc_data_per_stack.tree_depth
+    )
+
+    vertices_numpy = _generate_vertices(dc_data_per_surface, debug, slice_object)
+    return dc_data_per_surface, vertices_numpy
+
+
+def _generate_vertices(dc_data_per_surface: DualContouringData, debug: bool, slice_object: slice) -> Any:
+    vertices: np.ndarray = generate_dual_contouring_vertices(
+        dc_data_per_stack=dc_data_per_surface,
+        slice_surface=slice_object,
+        debug=debug
+    )
+    vertices_numpy = BackendTensor.t.to_numpy(vertices)
+    return vertices_numpy
+
+
+def generate_dual_contouring_vertices(dc_data_per_stack: DualContouringData, slice_surface: slice, debug: bool = False):
+    # @off
+    n_edges = dc_data_per_stack.n_edges
+    valid_edges = dc_data_per_stack.valid_edges
+    valid_voxels = dc_data_per_stack.valid_voxels
+    xyz_on_edge = dc_data_per_stack.xyz_on_edge[slice_surface]
+    gradients = dc_data_per_stack.gradients[slice_surface]
+    # @on
+
+    # * Coordinates for all posible edges (12) and 3 dummy edges_normals in the center
+    edges_xyz = BackendTensor.tfnp.zeros((n_edges, 15, 3), dtype=BackendTensor.dtype_obj)
+    valid_edges = valid_edges > 0
+    edges_xyz[:, :12][valid_edges] = xyz_on_edge
+
+    # Normals
+    edges_normals = BackendTensor.tfnp.zeros((n_edges, 15, 3), dtype=BackendTensor.dtype_obj)
+    edges_normals[:, :12][valid_edges] = gradients
+
+    if OLD_METHOD := False:
+        # ! Moureze model does not seems to work with the new method
+        # ! This branch is all nans at least with ch1_1 model
+        bias_xyz = BackendTensor.tfnp.copy(edges_xyz[:, :12])
+        isclose = BackendTensor.tfnp.isclose(bias_xyz, 0)
+        bias_xyz[isclose] = BackendTensor.tfnp.nan  # zero values to nans
+        mass_points = BackendTensor.tfnp.nanmean(bias_xyz, axis=1)  # Mean ignoring nans
+    else:  # ? This is actually doing something
+        bias_xyz = BackendTensor.tfnp.copy(edges_xyz[:, :12])
+        if BackendTensor.engine_backend == AvailableBackends.PYTORCH:
+            # PyTorch doesn't have masked arrays, so we'll use a different approach
+            mask = bias_xyz == 0
+            # Replace zeros with NaN for mean calculation
+            bias_xyz_masked = BackendTensor.tfnp.where(mask, float('nan'), bias_xyz)
+            mass_points = BackendTensor.tfnp.nanmean(bias_xyz_masked, axis=1)
+        else:
+            # NumPy approach with masked arrays
+            bias_xyz = BackendTensor.tfnp.to_numpy(bias_xyz)
+            import numpy as np
+            mask = bias_xyz == 0
+            masked_arr = np.ma.masked_array(bias_xyz, mask)
+            mass_points = masked_arr.mean(axis=1)
+            mass_points = BackendTensor.tfnp.array(mass_points)
+
+    edges_xyz[:, 12] = mass_points
+    edges_xyz[:, 13] = mass_points
+    edges_xyz[:, 14] = mass_points
+
+    BIAS_STRENGTH = 1
+
+    bias_x = BackendTensor.tfnp.array([BIAS_STRENGTH, 0, 0], dtype=BackendTensor.dtype_obj)
+    bias_y = BackendTensor.tfnp.array([0, BIAS_STRENGTH, 0], dtype=BackendTensor.dtype_obj)
+    bias_z = BackendTensor.tfnp.array([0, 0, BIAS_STRENGTH], dtype=BackendTensor.dtype_obj)
+
+    edges_normals[:, 12] = bias_x
+    edges_normals[:, 13] = bias_y
+    edges_normals[:, 14] = bias_z
+
+    # Remove unused voxels
+    edges_xyz = edges_xyz[valid_voxels]
+    edges_normals = edges_normals[valid_voxels]
+
+    # Compute LSTSQS in all voxels at the same time
+    A = edges_normals
+    b = (A * edges_xyz).sum(axis=2)
+
+    if BackendTensor.engine_backend == AvailableBackends.PYTORCH:
+        transpose_shape = (2, 1, 0)  # For PyTorch: (batch, dim2, dim1)
+    else:
+        transpose_shape = (0, 2, 1)  # For NumPy: (batch, dim2, dim1)
+
+    term1 = BackendTensor.tfnp.einsum("ijk, ilj->ikl", A, BackendTensor.tfnp.transpose(A, transpose_shape))
+    term2 = BackendTensor.tfnp.linalg.inv(term1)
+    term3 = BackendTensor.tfnp.einsum("ijk,ik->ij", BackendTensor.tfnp.transpose(A, transpose_shape), b)
+    vertices = BackendTensor.tfnp.einsum("ijk, ij->ik", term2, term3)
+
+    if debug:
+        dc_data_per_stack.bias_center_mass = edges_xyz[:, 12:].reshape(-1, 3)
+        dc_data_per_stack.bias_normals = edges_normals[:, 12:].reshape(-1, 3)
+
+    return vertices

gempy_engine/modules/dual_contouring/_gen_vertices.py

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