Skip to content

Commit 901dbed

Browse files
LeguarkLeguark
authored
[FEATURE] Add mesh extraction module using marching cubes algorithm (#1005)
TODO: Needs new version of gempy viewer # Description Added a new module for mesh extraction using the marching cubes algorithm. The implementation provides a clean way to extract meshes for all structural elements in a geological model. The new `mesh_extraction` module contains functions to: - Extract meshes for all structural elements in a model - Process individual elements with proper coordinate transformations - Handle the conversion from scalar fields to 3D meshes Refactored the test_marching_cubes.py to use the new module, significantly simplifying the test code by removing the manual implementation that was previously used. Relates to #mesh-extraction-enhancement # Checklist - [x] My code uses type hinting for function and method arguments and return values. - [x] I have created tests which cover my code. - [x] The test code either 1. demonstrates at least one valuable use case (e.g. integration tests) or 2. verifies that outputs are as expected for given inputs (e.g. unit tests). - [ ] New tests pass locally with my changes.
2 parents 3e4614e + 96d8df2 commit 901dbed

File tree

4 files changed

+99
-140
lines changed

4 files changed

+99
-140
lines changed

‎gempy/core/data/geo_model.py

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -147,6 +147,8 @@ def solutions(self) -> Solutions:
147147

148148
@solutions.setter
149149
def solutions(self, value):
150+
# * This is set from the gempy engine
151+
150152
self._solutions = value
151153

152154
# * Set solutions per group

‎gempy/modules/mesh_extranction/__init__.py

Whitespace-only changes.

‎gempy/modules/mesh_extranction/marching_cubes.py

Lines changed: 88 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,88 @@
1+
import numpy as np
2+
from typing import Optional
3+
from skimage import measure
4+
from gempy_engine.core.data.interp_output import InterpOutput
5+
from gempy_engine.core.data.raw_arrays_solution import RawArraysSolution
6+
7+
from gempy.core.data import GeoModel, StructuralElement, StructuralGroup
8+
from gempy.core.data.grid_modules import RegularGrid
9+
10+
11+
def set_meshes_with_marching_cubes(model: GeoModel) -> None:
12+
"""Extract meshes for all structural elements using the marching cubes algorithm.
13+
14+
Parameters
15+
----------
16+
model : GeoModel
17+
The geological model containing solutions and structural elements.
18+
19+
Raises
20+
------
21+
ValueError
22+
If the model solutions do not contain dense grid data.
23+
"""
24+
# Verify that solutions contain dense grid data
25+
if (model.solutions is None or
26+
model.solutions.block_solution_type != RawArraysSolution.BlockSolutionType.DENSE_GRID):
27+
raise ValueError("Model solutions must contain dense grid data for mesh extraction.")
28+
29+
regular_grid: RegularGrid = model.grid.regular_grid
30+
structural_groups: list[StructuralGroup] = model.structural_frame.structural_groups
31+
32+
if not model.solutions.octrees_output or not model.solutions.octrees_output[0].outputs_centers:
33+
raise ValueError("No interpolation outputs available for mesh extraction.")
34+
35+
output_lvl0: list[InterpOutput] = model.solutions.octrees_output[0].outputs_centers
36+
37+
for e, structural_group in enumerate(structural_groups):
38+
if e >= len(output_lvl0):
39+
continue
40+
41+
output_group: InterpOutput = output_lvl0[e]
42+
scalar_field_matrix = output_group.exported_fields_dense_grid.scalar_field
43+
44+
for element in structural_group.elements:
45+
extract_mesh_for_element(
46+
structural_element=element,
47+
regular_grid=regular_grid,
48+
scalar_field=scalar_field_matrix
49+
)
50+
51+
52+
def extract_mesh_for_element(structural_element: StructuralElement,
53+
regular_grid: RegularGrid,
54+
scalar_field: np.ndarray,
55+
mask: Optional[np.ndarray] = None) -> None:
56+
"""Extract a mesh for a single structural element using marching cubes.
57+
58+
Parameters
59+
----------
60+
structural_element : StructuralElement
61+
The structural element for which to extract a mesh.
62+
regular_grid : RegularGrid
63+
The regular grid defining the spatial discretization.
64+
scalar_field : np.ndarray
65+
The scalar field used for isosurface extraction.
66+
mask : np.ndarray, optional
67+
Optional mask to restrict the mesh extraction to specific regions.
68+
"""
69+
# Apply mask if provided
70+
volume = scalar_field.reshape(regular_grid.resolution)
71+
if mask is not None:
72+
volume = volume * mask
73+
74+
# Extract mesh using marching cubes
75+
verts, faces, _, _ = measure.marching_cubes(
76+
volume=volume,
77+
level=structural_element.scalar_field,
78+
spacing=(regular_grid.dx, regular_grid.dy, regular_grid.dz)
79+
)
80+
81+
# Adjust vertices to correct coordinates in the model's extent
82+
verts = (verts + [regular_grid.extent[0],
83+
regular_grid.extent[2],
84+
regular_grid.extent[4]])
85+
86+
# Store mesh in the structural element
87+
structural_element.vertices = verts
88+
structural_element.edges = faces

‎test/test_modules/test_marching_cubes.py

Lines changed: 9 additions & 140 deletions
Original file line numberDiff line numberDiff line change
@@ -4,35 +4,11 @@
44
import gempy as gp
55
from gempy.core.data.enumerators import ExampleModel
66
from gempy.core.data.grid_modules import RegularGrid
7+
from gempy.modules.mesh_extranction import marching_cubes
78
from gempy.optional_dependencies import require_gempy_viewer
89

9-
from skimage import measure
10-
1110
PLOT = True
1211

13-
def marching_cubes(block, elements, spacing, extent):
14-
"""
15-
Extract the surface meshes using marching cubes
16-
Args:
17-
block (np.array): The block to extract the surface meshes from.
18-
elements (list): IDs of unique structural elements in model
19-
spacing (tuple): The spacing between grid points in the block.
20-
21-
Returns:
22-
mc_vertices (list): Vertices of the surface meshes.
23-
mc_edges (list): Edges of the surface meshes.
24-
"""
25-
26-
# Extract the surface meshes using marching cubes
27-
mc_vertices = []
28-
mc_edges = []
29-
for i in range(0, len(elements)):
30-
verts, faces, _, _ = measure.marching_cubes(block, i,
31-
spacing=spacing)
32-
mc_vertices.append(verts + [extent[0], extent[2], extent[4]])
33-
mc_edges.append(faces)
34-
return mc_vertices, mc_edges
35-
3612

3713
def test_marching_cubes_implementation():
3814
model = gp.generate_example_model(ExampleModel.COMBINATION, compute_model=False)
@@ -50,6 +26,7 @@ def test_marching_cubes_implementation():
5026
reset=True
5127
)
5228

29+
model.interpolation_options.evaluation_options.number_octree_levels = 1
5330
model.interpolation_options.evaluation_options.mesh_extraction = False # * Not extracting the mesh with dual contouring
5431
gp.compute_model(model)
5532

@@ -60,121 +37,13 @@ def test_marching_cubes_implementation():
6037

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

63-
# TODO: Maybe to complicated because it includes accounting for faults, multiple elements in groups
64-
# and transformation to real coordinates
65-
66-
# Empty lists to store vertices and edges
67-
mc_vertices = []
68-
mc_edges = []
69-
70-
# Boolean list of fault groups
71-
faults = model.structural_frame.group_is_fault
72-
73-
# MC for faults, directly on fault block not on scalar field
74-
if faults is not None:
75-
# TODO: This should also use the scalar fields probably
76-
for i in np.unique(model.solutions.raw_arrays.fault_block)[:-1]:
77-
fault_block = model.solutions.raw_arrays.fault_block.reshape(model.grid.regular_grid.resolution)
78-
verts, faces, _, _ = measure.marching_cubes(fault_block,
79-
i,
80-
spacing=(model.grid.regular_grid.dx,
81-
model.grid.regular_grid.dy,
82-
model.grid.regular_grid.dz))
83-
mc_vertices.append(verts + [model.grid.regular_grid.extent[0],
84-
model.grid.regular_grid.extent[2],
85-
model.grid.regular_grid.extent[4]])
86-
mc_edges.append(faces)
87-
else:
88-
pass
89-
90-
# Extract scalar field values for elements
91-
scalar_values = model.solutions.raw_arrays.scalar_field_at_surface_points
92-
93-
# Get indices of non fault elements
94-
if faults is not None:
95-
false_indices = [i for i, fault in enumerate(faults) if not fault]
96-
else:
97-
false_indices = np.arange(len(model.structural_frame.structural_groups))
98-
99-
# Extract marching cubes for non fault elements
100-
for idx in false_indices:
101-
102-
# Get correct scalar field for structural group
103-
scalar_field = model.solutions.raw_arrays.scalar_field_matrix[idx].reshape(model.grid.regular_grid.resolution)
104-
105-
# Extract marching cubes for each scalar value for all elements of a group
106-
for i in range(len(scalar_values[idx])):
107-
verts, faces, _, _ = measure.marching_cubes(scalar_field, scalar_values[idx][i],
108-
spacing=(model.grid.regular_grid.dx,
109-
model.grid.regular_grid.dy,
110-
model.grid.regular_grid.dz))
111-
112-
mc_vertices.append(verts + [model.grid.regular_grid.extent[0],
113-
model.grid.regular_grid.extent[2],
114-
model.grid.regular_grid.extent[4]])
115-
mc_edges.append(faces)
116-
117-
# Reorder everything correctly if faults exist
118-
# TODO: All of the following is just complicated code to reorder the elements to match the order of the elements
119-
# in the structural frame, probably unnecessary in gempy strucuture
120-
#
121-
# if faults is not None:
122-
#
123-
# # TODO: This is a very convoluted way to get a boolean list of faults per element
124-
# bool_list = np.zeros(4, dtype=bool)
125-
# for i in range(len(model.structural_frame.structural_groups)):
126-
# print(i)
127-
# if model.structural_frame.group_is_fault[i]:
128-
# for j in range(len(model.structural_frame.structural_groups[i].elements)):
129-
# bool_list[i + j] = True
130-
# if not model.structural_frame.group_is_fault[i]:
131-
# for k in range(len(model.structural_frame.structural_groups[i].elements)):
132-
# bool_list[i + k] = False
133-
#
134-
# true_count = sum(bool_list)
135-
#
136-
# # Split arr_list into two parts
137-
# true_elements_vertices = mc_vertices[:true_count]
138-
# false_elements_vertices = mc_vertices[true_count:]
139-
# true_elements_edges = mc_edges[:true_count]
140-
# false_elements_edges = mc_edges[true_count:]
141-
#
142-
# # Create a new list to store reordered elements
143-
# mc_vertices = []
144-
# mc_edges = []
145-
#
146-
# # Iterator for both true and false elements
147-
# true_idx, false_idx = 0, 0
148-
#
149-
# # Populate reordered_list based on bool_list
150-
# for is_true in bool_list:
151-
# if is_true:
152-
# mc_vertices.append(true_elements_vertices[true_idx] + [model.grid.regular_grid.extent[0],
153-
# model.grid.regular_grid.extent[2],
154-
# model.grid.regular_grid.extent[4]])
155-
# mc_edges.append(true_elements_edges[true_idx])
156-
# true_idx += 1
157-
# else:
158-
# mc_vertices.append(false_elements_vertices[false_idx] + [model.grid.regular_grid.extent[0],
159-
# model.grid.regular_grid.extent[2],
160-
# model.grid.regular_grid.extent[4]])
161-
# mc_edges.append(false_elements_edges[false_idx])
162-
# false_idx += 1
40+
marching_cubes.set_meshes_with_marching_cubes(model)
16341

16442
if PLOT:
16543
gpv = require_gempy_viewer()
166-
# gtv: gpv.GemPyToVista = gpv.plot_3d(model, show_data=True, image=True)
167-
import pyvista as pv
168-
# pyvista_plotter: pv.Plotter = gtv.p
169-
170-
# TODO: This opens interactive window as of now
171-
pyvista_plotter = pv.Plotter()
172-
173-
# Add the meshes to the plot
174-
for i in range(len(mc_vertices)):
175-
pyvista_plotter.add_mesh(
176-
pv.PolyData(mc_vertices[i],
177-
np.insert(mc_edges[i], 0, 3, axis=1).ravel()),
178-
color='blue')
179-
180-
pyvista_plotter.show()
44+
gtv: gpv.GemPyToVista = gpv.plot_3d(
45+
model=model,
46+
show_data=True,
47+
image=True,
48+
show=True
49+
)

0 commit comments

Comments
 (0)