[CLN] Moving vertex_overlap code to the module

Author: Leguark

gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py

@@ -19,7 +19,7 @@
 from ...core.data.octree_level import OctreeLevel
 from ...core.utils import gempy_profiler_decorator
 from ...modules.dual_contouring.dual_contouring_interface import (find_intersection_on_edge, get_triangulation_codes,
-                                                                  get_masked_codes, mask_generation)
+                                                                  get_masked_codes, mask_generation,apply_faults_vertex_overlap)
 
 
 @gempy_profiler_decorator
@@ -100,7 +100,7 @@ def dual_contouring_multi_scalar(
     # endregion
 
     # region Vertex gen and triangulation
-    foo = []
+    left_right_per_mesh = []
     # Generate meshes for each scalar field
     for n_scalar_field in range(data_descriptor.stack_structure.n_stacks):
         output: InterpOutput = octree_leaves.outputs_centers[n_scalar_field]
@@ -126,84 +126,22 @@ def dual_contouring_multi_scalar(
         )
 
         for m in meshes:
-            foo.append(m.left_right)
+            left_right_per_mesh.append(m.left_right)
 
         # TODO: If the order of the meshes does not match the order of scalar_field_at_surface points, reorder them here
         if meshes is not None:
             all_meshes.extend(meshes)
 
     # endregion
-    # Check for repeated voxels across stacks
-    if (options.debug or len(all_left_right_codes) > 1) and False:
-        voxel_overlaps = find_repeated_voxels_across_stacks(foo)
-        if voxel_overlaps and options.debug:
-            print(f"Found voxel overlaps between stacks: {voxel_overlaps}")
-            _f(all_meshes, 1, 0, voxel_overlaps)
-            _f(all_meshes, 2, 0, voxel_overlaps)
-            _f(all_meshes, 3, 0, voxel_overlaps)
-            _f(all_meshes, 4, 0, voxel_overlaps)
-            _f(all_meshes, 5, 0, voxel_overlaps)
+    if (options.debug or len(all_left_right_codes) > 1) and True:
+        apply_faults_vertex_overlap(all_meshes, data_descriptor, left_right_per_mesh)
 
     return all_meshes
 
+    # ... existing code ...
 
-def _f(all_meshes: list[DualContouringMesh], destination: int, origin: int, voxel_overlaps: dict):
-    key = f"stack_{origin}_vs_stack_{destination}"
-    all_meshes[destination].vertices[voxel_overlaps[key]["indices_in_stack_j"]] = all_meshes[origin].vertices[voxel_overlaps[key]["indices_in_stack_i"]]
 
 
-def find_repeated_voxels_across_stacks(all_left_right_codes: List[np.ndarray]) -> dict:
-    """
-    Find repeated voxels using NumPy operations - better for very large arrays.
-
-    Args:
-        all_left_right_codes: List of left_right_codes arrays, one per stack
-
-    Returns:
-        Dictionary with detailed overlap analysis
-    """
-
-    if not all_left_right_codes:
-        return {}
-
-    # Generate voxel codes for each stack
-
-    from gempy_engine.modules.dual_contouring.fancy_triangulation import _StaticTriangulationData
-    stack_codes = []
-    for left_right_codes in all_left_right_codes:
-        if left_right_codes.size > 0:
-            voxel_codes = (left_right_codes * _StaticTriangulationData.get_pack_directions_into_bits()).sum(axis=1)
-            stack_codes.append(voxel_codes)
-        else:
-            stack_codes.append(np.array([]))
-
-    overlaps = {}
-
-    # Check each pair of stacks
-    for i in range(len(stack_codes)):
-        for j in range(i + 1, len(stack_codes)):
-            if stack_codes[i].size == 0 or stack_codes[j].size == 0:
-                continue
-
-            # Find common voxel codes using numpy
-            common_codes = np.intersect1d(stack_codes[i], stack_codes[j])
-
-            if len(common_codes) > 0:
-                # Get indices of common voxels in each stack
-                indices_i = np.isin(stack_codes[i], common_codes)
-                indices_j = np.isin(stack_codes[j], common_codes)
-
-                overlaps[f"stack_{i}_vs_stack_{j}"] = {
-                        'common_voxel_codes'   : common_codes,
-                        'count'                : len(common_codes),
-                        'indices_in_stack_i'   : np.where(indices_i)[0],
-                        'indices_in_stack_j'   : np.where(indices_j)[0],
-                        'common_binary_codes_i': all_left_right_codes[i][indices_i],
-                        'common_binary_codes_j': all_left_right_codes[j][indices_j]
-                }
-
-    return overlaps
-
 
 def _validate_stack_relations(data_descriptor: InputDataDescriptor, n_scalar_field: int) -> None:
     """

gempy_engine/modules/dual_contouring/_vertex_overlap.py

@@ -0,0 +1,164 @@
+from typing import List
+
+import numpy as np
+
+from ...core.data.dual_contouring_mesh import DualContouringMesh
+
+
+
+def _apply_fault_relations_to_overlaps(
+        all_meshes: List[DualContouringMesh],
+        faults_relations: np.ndarray,
+        voxel_overlaps: dict,
+        n_stacks: int
+) -> None:
+    """
+    Apply fault relations to voxel overlaps by updating mesh vertices.
+
+    Args:
+        all_meshes: List of dual contouring meshes
+        faults_relations: Boolean matrix indicating fault relationships between stacks
+        voxel_overlaps: Dictionary containing overlap information between stacks
+        n_stacks: Total number of stacks
+    """
+    if faults_relations is None:
+        return
+
+    # Calculate mesh indices offset for each stack
+    mesh_indices_offset = _calculate_mesh_indices_offset(all_meshes, n_stacks)
+
+    # Iterate through fault relations matrix
+    for origin_stack in range(n_stacks):
+        for destination_stack in range(n_stacks):
+            # If there's a fault relation from origin to destination
+            if faults_relations[origin_stack, destination_stack]:
+                overlap_key = f"stack_{origin_stack}_vs_stack_{destination_stack}"
+
+                # Check if there are actual overlaps between these stacks
+                if overlap_key in voxel_overlaps:
+                    _apply_vertex_sharing(
+                        all_meshes,
+                        origin_stack,
+                        destination_stack,
+                        voxel_overlaps[overlap_key],
+                        mesh_indices_offset
+                    )
+
+
+def _calculate_mesh_indices_offset(all_meshes: List[DualContouringMesh], n_stacks: int) -> List[int]:
+    """
+    Calculate the starting mesh index for each stack.
+
+    Args:
+        all_meshes: List of all dual contouring meshes
+        n_stacks: Total number of stacks
+
+    Returns:
+        List of starting mesh indices for each stack
+    """
+    # For now, assume each stack has one mesh (this may need adjustment based on actual structure)
+    # This is a simplified approach - you may need to adjust based on how meshes are organized
+    mesh_indices_offset = list(range(n_stacks))
+    return mesh_indices_offset
+
+
+def _apply_vertex_sharing(
+        all_meshes: List[DualContouringMesh],
+        origin_stack: int,
+        destination_stack: int,
+        overlap_data: dict,
+        mesh_indices_offset: List[int]
+) -> None:
+    """
+    Apply vertex sharing between origin and destination meshes based on overlap data.
+
+    Args:
+        all_meshes: List of dual contouring meshes
+        origin_stack: Stack index that serves as the source of vertices
+        destination_stack: Stack index that receives vertices from origin
+        overlap_data: Dictionary containing indices and overlap information
+        mesh_indices_offset: Starting mesh index for each stack
+    """
+    origin_mesh_idx = mesh_indices_offset[origin_stack]
+    destination_mesh_idx = mesh_indices_offset[destination_stack]
+
+    # Ensure mesh indices are valid
+    if (origin_mesh_idx >= len(all_meshes) or
+            destination_mesh_idx >= len(all_meshes)):
+        return
+
+    # Apply the vertex sharing (same logic as original _f function)
+    origin_mesh = all_meshes[origin_mesh_idx]
+    destination_mesh = all_meshes[destination_mesh_idx]
+
+    indices_in_origin = overlap_data["indices_in_stack_i"]
+    indices_in_destination = overlap_data["indices_in_stack_j"]
+
+    destination_mesh.vertices[indices_in_destination] = origin_mesh.vertices[indices_in_origin]
+
+
+def _f(all_meshes: list[DualContouringMesh], destination: int, origin: int, voxel_overlaps: dict):
+    """
+    Legacy function - kept for backward compatibility.
+    Consider using _apply_fault_relations_to_overlaps for new implementations.
+    """
+    key = f"stack_{origin}_vs_stack_{destination}"
+    if key in voxel_overlaps:
+        all_meshes[destination].vertices[voxel_overlaps[key]["indices_in_stack_j"]] = all_meshes[origin].vertices[voxel_overlaps[key]["indices_in_stack_i"]]
+
+# def _f(all_meshes: list[DualContouringMesh], destination: int, origin: int, voxel_overlaps: dict):
+#     key = f"stack_{origin}_vs_stack_{destination}"
+#     all_meshes[destination].vertices[voxel_overlaps[key]["indices_in_stack_j"]] = all_meshes[origin].vertices[voxel_overlaps[key]["indices_in_stack_i"]]
+
+
+def find_repeated_voxels_across_stacks(all_left_right_codes: List[np.ndarray]) -> dict:
+    """
+    Find repeated voxels using NumPy operations - better for very large arrays.
+
+    Args:
+        all_left_right_codes: List of left_right_codes arrays, one per stack
+
+    Returns:
+        Dictionary with detailed overlap analysis
+    """
+
+    if not all_left_right_codes:
+        return {}
+
+    # Generate voxel codes for each stack
+
+    from gempy_engine.modules.dual_contouring.fancy_triangulation import _StaticTriangulationData
+    stack_codes = []
+    for left_right_codes in all_left_right_codes:
+        if left_right_codes.size > 0:
+            voxel_codes = (left_right_codes * _StaticTriangulationData.get_pack_directions_into_bits()).sum(axis=1)
+            stack_codes.append(voxel_codes)
+        else:
+            stack_codes.append(np.array([]))
+
+    overlaps = {}
+
+    # Check each pair of stacks
+    for i in range(len(stack_codes)):
+        for j in range(i + 1, len(stack_codes)):
+            if stack_codes[i].size == 0 or stack_codes[j].size == 0:
+                continue
+
+            # Find common voxel codes using numpy
+            common_codes = np.intersect1d(stack_codes[i], stack_codes[j])
+
+            if len(common_codes) > 0:
+                # Get indices of common voxels in each stack
+                indices_i = np.isin(stack_codes[i], common_codes)
+                indices_j = np.isin(stack_codes[j], common_codes)
+
+                overlaps[f"stack_{i}_vs_stack_{j}"] = {
+                        'common_voxel_codes'   : common_codes,
+                        'count'                : len(common_codes),
+                        'indices_in_stack_i'   : np.where(indices_i)[0],
+                        'indices_in_stack_j'   : np.where(indices_j)[0],
+                        'common_binary_codes_i': all_left_right_codes[i][indices_i],
+                        'common_binary_codes_j': all_left_right_codes[j][indices_j]
+                }
+
+    return overlaps

gempy_engine/modules/dual_contouring/dual_contouring_interface.py

@@ -3,9 +3,11 @@
 
 import numpy as np
 
+from ._vertex_overlap import find_repeated_voxels_across_stacks, _apply_fault_relations_to_overlaps
 from .fancy_triangulation import get_left_right_array
 from ...core.backend_tensor import BackendTensor
 from ...core.data import InterpolationOptions
+from ...core.data.dual_contouring_mesh import DualContouringMesh
 from ...core.data.input_data_descriptor import InputDataDescriptor
 from ...core.data.interp_output import InterpOutput
 from ...core.data.octree_level import OctreeLevel
@@ -109,8 +111,7 @@ def get_triangulation_codes(octree_list: List[OctreeLevel], options: Interpolati
             raise ValueError("Invalid combination of options")
 
 
-
-def get_masked_codes(left_right_codes: np.ndarray | None, mask: np.ndarray | None) -> np.ndarray | None:
+def get_masked_codes(left_right_codes: np.ndarray, mask: np.ndarray) -> np.ndarray:
     """
     Apply mask to left-right codes if both are available.
     
@@ -133,7 +134,7 @@ def get_masked_codes(left_right_codes: np.ndarray | None, mask: np.ndarray | Non
 def mask_generation(
         octree_leaves: OctreeLevel,
         masking_option: MeshExtractionMaskingOptions
-) -> np.ndarray | None:
+) -> np.ndarray:
     """
     Generate masks for mesh extraction based on masking options and stack relations.
     
@@ -196,4 +197,13 @@ def mask_generation(
     return mask_matrix
 
 
-# endregion
+# endregion
+def apply_faults_vertex_overlap(all_meshes: list[DualContouringMesh],
+                                data_descriptor: InputDataDescriptor, 
+                                left_right_per_mesh: list[np.ndarray]):
+    faults_relations = data_descriptor.stack_structure.faults_relations
+    voxel_overlaps = find_repeated_voxels_across_stacks(left_right_per_mesh)
+    
+    if voxel_overlaps:
+        print(f"Found voxel overlaps between stacks: {voxel_overlaps}")
+        _apply_fault_relations_to_overlaps(all_meshes, faults_relations, voxel_overlaps, data_descriptor.stack_structure.n_stacks)