Skip to content

[FEATURE] Add mesh extraction module using marching cubes algorithm | GEN-12031 #1005

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Mar 23, 2025
Merged

[FEATURE] Add mesh extraction module using marching cubes algorithm | GEN-12031 #1005

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
76756ad
Adding mc
Leguark Mar 22, 2025
d4b4277
[ENH/CLN/TEST/WIP] Moving marching cubes code to a module
Leguark Mar 22, 2025
96d8df2
[ENH] Adding first version of marching cubes
Leguark Mar 22, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
[ENH/CLN/TEST/WIP] Moving marching cubes code to a module
  • Loading branch information
@Leguark
Leguark committed Mar 23, 2025
commit d4b427746f1b2535d453fd7088051df19cee7b50
64 changes: 64 additions & 0 deletions gempy/modules/mesh_extranction/marching_cubes.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,64 @@
import numpy as np
from skimage import measure


def compute_marching_cubes(model):
# Empty lists to store vertices and edges
mc_vertices = []
mc_edges = []
# Boolean list of fault groups
faults = model.structural_frame.group_is_fault
# MC for faults, directly on fault block not on scalar field
if faults is not None:
_extract_fault_mesh(mc_edges, mc_vertices, model)
else:
pass

# Extract scalar field values for elements
scalar_values = model.solutions.raw_arrays.scalar_field_at_surface_points

# Get indices of non fault elements
false_indices = _get_lithology_idx(faults, model)
# Extract marching cubes for non fault elements
for idx in false_indices:
_extract_meshes_for_lithologies(idx, mc_edges, mc_vertices, model, scalar_values)
return mc_edges, mc_vertices


def _extract_meshes_for_lithologies(idx, mc_edges, mc_vertices, model, scalar_values):
# Get correct scalar field for structural group
scalar_field = model.solutions.raw_arrays.scalar_field_matrix[idx].reshape(model.grid.regular_grid.resolution)
# Extract marching cubes for each scalar value for all elements of a group
for i in range(len(scalar_values[idx])):
verts, faces, _, _ = measure.marching_cubes(scalar_field, scalar_values[idx][i],
spacing=(model.grid.regular_grid.dx,
model.grid.regular_grid.dy,
model.grid.regular_grid.dz))

mc_vertices.append(verts + [model.grid.regular_grid.extent[0],
model.grid.regular_grid.extent[2],
model.grid.regular_grid.extent[4]])
mc_edges.append(faces)


def _get_lithology_idx(faults, model):
if faults is not None:
false_indices = [i for i, fault in enumerate(faults) if not fault]
else:
false_indices = np.arange(len(model.structural_frame.structural_groups))
return false_indices


def _extract_fault_mesh(mc_edges, mc_vertices, model):
# TODO: This should also use the scalar fields probably
for i in np.unique(model.solutions.raw_arrays.fault_block)[:-1]:
fault_block = model.solutions.raw_arrays.fault_block.reshape(model.grid.regular_grid.resolution)
verts, faces, _, _ = measure.marching_cubes(fault_block,
i,
spacing=(model.grid.regular_grid.dx,
model.grid.regular_grid.dy,
model.grid.regular_grid.dz))
mc_vertices.append(verts + [model.grid.regular_grid.extent[0],
model.grid.regular_grid.extent[2],
model.grid.regular_grid.extent[4]])
mc_edges.append(faces)
139 changes: 11 additions & 128 deletions test/test_modules/test_marching_cubes.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,36 +4,12 @@
import gempy as gp
from gempy.core.data.enumerators import ExampleModel
from gempy.core.data.grid_modules import RegularGrid
from gempy.modules.mesh_extranction import marching_cubes
from gempy.optional_dependencies import require_gempy_viewer

from skimage import measure

PLOT = True

def marching_cubes(block, elements, spacing, extent):
"""
Extract the surface meshes using marching cubes
Args:
block (np.array): The block to extract the surface meshes from.
elements (list): IDs of unique structural elements in model
spacing (tuple): The spacing between grid points in the block.

Returns:
mc_vertices (list): Vertices of the surface meshes.
mc_edges (list): Edges of the surface meshes.
"""

# Extract the surface meshes using marching cubes
mc_vertices = []
mc_edges = []
for i in range(0, len(elements)):
verts, faces, _, _ = measure.marching_cubes(block, i,
spacing=spacing)
mc_vertices.append(verts + [extent[0], extent[2], extent[4]])
mc_edges.append(faces)
return mc_vertices, mc_edges


def test_marching_cubes_implementation():
model = gp.generate_example_model(ExampleModel.COMBINATION, compute_model=False)

Expand All @@ -60,115 +36,20 @@ def test_marching_cubes_implementation():

assert arrays.scalar_field_matrix.shape == (3, 8_000) # * 3 surfaces, 8000 points

# TODO: Maybe to complicated because it includes accounting for faults, multiple elements in groups
# and transformation to real coordinates

# Empty lists to store vertices and edges
mc_vertices = []
mc_edges = []

# Boolean list of fault groups
faults = model.structural_frame.group_is_fault

# MC for faults, directly on fault block not on scalar field
if faults is not None:
# TODO: This should also use the scalar fields probably
for i in np.unique(model.solutions.raw_arrays.fault_block)[:-1]:
fault_block = model.solutions.raw_arrays.fault_block.reshape(model.grid.regular_grid.resolution)
verts, faces, _, _ = measure.marching_cubes(fault_block,
i,
spacing=(model.grid.regular_grid.dx,
model.grid.regular_grid.dy,
model.grid.regular_grid.dz))
mc_vertices.append(verts + [model.grid.regular_grid.extent[0],
model.grid.regular_grid.extent[2],
model.grid.regular_grid.extent[4]])
mc_edges.append(faces)
else:
pass

# Extract scalar field values for elements
scalar_values = model.solutions.raw_arrays.scalar_field_at_surface_points

# Get indices of non fault elements
if faults is not None:
false_indices = [i for i, fault in enumerate(faults) if not fault]
else:
false_indices = np.arange(len(model.structural_frame.structural_groups))

# Extract marching cubes for non fault elements
for idx in false_indices:

# Get correct scalar field for structural group
scalar_field = model.solutions.raw_arrays.scalar_field_matrix[idx].reshape(model.grid.regular_grid.resolution)

# Extract marching cubes for each scalar value for all elements of a group
for i in range(len(scalar_values[idx])):
verts, faces, _, _ = measure.marching_cubes(scalar_field, scalar_values[idx][i],
spacing=(model.grid.regular_grid.dx,
model.grid.regular_grid.dy,
model.grid.regular_grid.dz))

mc_vertices.append(verts + [model.grid.regular_grid.extent[0],
model.grid.regular_grid.extent[2],
model.grid.regular_grid.extent[4]])
mc_edges.append(faces)

# Reorder everything correctly if faults exist
# TODO: All of the following is just complicated code to reorder the elements to match the order of the elements
# in the structural frame, probably unnecessary in gempy strucuture
#
# if faults is not None:
#
# # TODO: This is a very convoluted way to get a boolean list of faults per element
# bool_list = np.zeros(4, dtype=bool)
# for i in range(len(model.structural_frame.structural_groups)):
# print(i)
# if model.structural_frame.group_is_fault[i]:
# for j in range(len(model.structural_frame.structural_groups[i].elements)):
# bool_list[i + j] = True
# if not model.structural_frame.group_is_fault[i]:
# for k in range(len(model.structural_frame.structural_groups[i].elements)):
# bool_list[i + k] = False
#
# true_count = sum(bool_list)
#
# # Split arr_list into two parts
# true_elements_vertices = mc_vertices[:true_count]
# false_elements_vertices = mc_vertices[true_count:]
# true_elements_edges = mc_edges[:true_count]
# false_elements_edges = mc_edges[true_count:]
#
# # Create a new list to store reordered elements
# mc_vertices = []
# mc_edges = []
#
# # Iterator for both true and false elements
# true_idx, false_idx = 0, 0
#
# # Populate reordered_list based on bool_list
# for is_true in bool_list:
# if is_true:
# mc_vertices.append(true_elements_vertices[true_idx] + [model.grid.regular_grid.extent[0],
# model.grid.regular_grid.extent[2],
# model.grid.regular_grid.extent[4]])
# mc_edges.append(true_elements_edges[true_idx])
# true_idx += 1
# else:
# mc_vertices.append(false_elements_vertices[false_idx] + [model.grid.regular_grid.extent[0],
# model.grid.regular_grid.extent[2],
# model.grid.regular_grid.extent[4]])
# mc_edges.append(false_elements_edges[false_idx])
# false_idx += 1
mc_edges, mc_vertices = marching_cubes.compute_marching_cubes(model)

if PLOT:
gpv = require_gempy_viewer()
# gtv: gpv.GemPyToVista = gpv.plot_3d(model, show_data=True, image=True)
gtv: gpv.GemPyToVista = gpv.plot_3d(
model=model,
show_data=True,
image=False,
show=False
)
import pyvista as pv
# pyvista_plotter: pv.Plotter = gtv.p

# TODO: This opens interactive window as of now
pyvista_plotter = pv.Plotter()
pyvista_plotter: pv.Pltter = gtv.p

# Add the meshes to the plot
for i in range(len(mc_vertices)):
Expand All @@ -178,3 +59,5 @@ def test_marching_cubes_implementation():
color='blue')

pyvista_plotter.show()