[ENH] Refactor marching cubes mesh extraction logic

Simplified and clarified the mesh extraction code by removing outdated comments, unused methods, and redundant handling of masks. Added proper mask reshaping and updated marching_cubes function parameters for improved reliability and efficiency. Adjusted test cases accordingly.

Author: Leguark

gempy/modules/mesh_extranction/marching_cubes.py

@@ -34,92 +34,27 @@ def set_meshes_with_marching_cubes(model: GeoModel) -> None:
 
     output_lvl0: list[InterpOutput] = model.solutions.octrees_output[0].outputs_centers
 
-    # TODO: How to get this properly in gempy
-    # get a list of indices of the lithological groups
-    lith_group_indices = []
-    fault_group_indices = []
-    index = 0
-    for i in model.structural_frame.structural_groups:
-        if i.is_fault:
-            fault_group_indices.append(index)
-        else:
-            lith_group_indices.append(index)
-        index += 1
-
-    # extract scalar field values at surface points
-    scalar_values = model.solutions.raw_arrays.scalar_field_at_surface_points
-
-    # TODO: Here I just get my own masks, cause the gempy masks dont work as expected
-    masks = _get_masking_arrays(lith_group_indices, model, scalar_values)
-
-    # TODO: Attribute of element.scalar_field was None, changed it to scalar field value of that element
-    #  This should probably be done somewhere else and maybe renamed to scalar_field_value?
-    #  This is just the most basic solution to be clear what I did
-    # _set_scalar_field_to_element(model, output_lvl0, structural_groups)
-
-    # Trying to use the exiting gempy masks
-    # masks = []
-    # masks.append(
-    #     np.ones_like(model.solutions.raw_arrays.scalar_field_matrix[0].reshape(model.grid.regular_grid.resolution),
-    #                  dtype=bool))
-    # for idx in lith_group_indices:
-    #     output_group: InterpOutput = output_lvl0[idx]
-    #     masks.append(output_group.mask_components[8:].reshape(model.grid.regular_grid.resolution))
-
-    non_fault_counter = 0
     for e, structural_group in enumerate(structural_groups):
         if e >= len(output_lvl0):
             continue
 
-        # Outdated?
-        # output_group: InterpOutput = output_lvl0[e]
-        # scalar_field_matrix = output_group.exported_fields_dense_grid.scalar_field
-
-        # Specify the correct scalar field, can be removed in the future
-        scalar_field = model.solutions.raw_arrays.scalar_field_matrix[e].reshape(model.grid.regular_grid.resolution)
-
-        # pick mask depending on whether the structural group is a fault or not
-        if structural_group.is_fault:
-            mask = np.ones_like(scalar_field, dtype=bool)
+        output_group: InterpOutput = output_lvl0[e]
+        scalar_field_matrix = output_group.exported_fields_dense_grid.scalar_field
+        if structural_group.is_fault is False:
+            slice_: slice = output_group.grid.dense_grid_slice
+            mask = output_group.combined_scalar_field.squeezed_mask_array[slice_]
         else:
-            mask = masks[non_fault_counter]  # TODO: I need the entry without faults here
-            non_fault_counter += 1
+            mask = np.ones_like(scalar_field_matrix, dtype=bool)
 
         for element in structural_group.elements:
             extract_mesh_for_element(
                 structural_element=element,
                 regular_grid=regular_grid,
-                scalar_field=scalar_field,
+                scalar_field=scalar_field_matrix,
                 mask=mask
             )
 
 
-# TODO: This should be set somewhere else
-def _set_scalar_field_to_element(model, output_lvl0, structural_groups):
-    element: StructuralElement
-    counter = 0
-    for e, structural_group in enumerate(structural_groups):
-        if e >= len(output_lvl0):
-            continue
-
-        for element in structural_group.elements:
-            element.scalar_field_at_interface = model.solutions.scalar_field_at_surface_points[counter]
-            counter += 1
-
-
-# TODO: This should be set somewhere else
-def _get_masking_arrays(lith_group_indices, model, scalar_values):
-    masks = []
-    masks.append(np.ones_like(model.solutions.raw_arrays.scalar_field_matrix[0].reshape(model.grid.regular_grid.resolution),
-                              dtype=bool))
-    for idx in lith_group_indices:
-        mask = model.solutions.raw_arrays.scalar_field_matrix[idx].reshape(model.grid.regular_grid.resolution) <= \
-               scalar_values[idx][-1]
-
-        masks.append(mask)
-    return masks
-
-
 def extract_mesh_for_element(structural_element: StructuralElement,
                              regular_grid: RegularGrid,
                              scalar_field: np.ndarray,
@@ -139,10 +74,12 @@ def extract_mesh_for_element(structural_element: StructuralElement,
     """
     # Extract mesh using marching cubes
     verts, faces, _, _ = measure.marching_cubes(
-        volume=scalar_field,
+        volume=scalar_field.reshape(regular_grid.resolution),
         level=structural_element.scalar_field_at_interface,
         spacing=(regular_grid.dx, regular_grid.dy, regular_grid.dz),
-        mask=mask
+        mask=mask.reshape(regular_grid.resolution) if mask is not None else None,
+        allow_degenerate=False,
+        method="lewiner"
     )
 
     # Adjust vertices to correct coordinates in the model's extent

test/test_modules/test_marching_cubes.py

@@ -52,6 +52,6 @@ def test_marching_cubes_implementation():
         gtv: gpv.GemPyToVista = gpv.plot_3d(
             model=model,
             show_data=True,
-            image=True,
+            image=False,
             show=True
         )