Skip to content

Commit 96d8df2

Browse files
committed
[ENH] Adding first version of marching cubes
1 parent d4b4277 commit 96d8df2

File tree

3 files changed

+85
-73
lines changed

3 files changed

+85
-73
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
Lines changed: 78 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -1,64 +1,88 @@
11
import numpy as np
2+
from typing import Optional
23
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
36

7+
from gempy.core.data import GeoModel, StructuralElement, StructuralGroup
8+
from gempy.core.data.grid_modules import RegularGrid
49

5-
def compute_marching_cubes(model):
6-
# Empty lists to store vertices and edges
7-
mc_vertices = []
8-
mc_edges = []
9-
# Boolean list of fault groups
10-
faults = model.structural_frame.group_is_fault
11-
# MC for faults, directly on fault block not on scalar field
12-
if faults is not None:
13-
_extract_fault_mesh(mc_edges, mc_vertices, model)
14-
else:
15-
pass
10+
11+
def set_meshes_with_marching_cubes(model: GeoModel) -> None:
12+
"""Extract meshes for all structural elements using the marching cubes algorithm.
1613
17-
# Extract scalar field values for elements
18-
scalar_values = model.solutions.raw_arrays.scalar_field_at_surface_points
14+
Parameters
15+
----------
16+
model : GeoModel
17+
The geological model containing solutions and structural elements.
1918
20-
# Get indices of non fault elements
21-
false_indices = _get_lithology_idx(faults, model)
22-
# Extract marching cubes for non fault elements
23-
for idx in false_indices:
24-
_extract_meshes_for_lithologies(idx, mc_edges, mc_vertices, model, scalar_values)
25-
return mc_edges, mc_vertices
26-
27-
28-
def _extract_meshes_for_lithologies(idx, mc_edges, mc_vertices, model, scalar_values):
29-
# Get correct scalar field for structural group
30-
scalar_field = model.solutions.raw_arrays.scalar_field_matrix[idx].reshape(model.grid.regular_grid.resolution)
31-
# Extract marching cubes for each scalar value for all elements of a group
32-
for i in range(len(scalar_values[idx])):
33-
verts, faces, _, _ = measure.marching_cubes(scalar_field, scalar_values[idx][i],
34-
spacing=(model.grid.regular_grid.dx,
35-
model.grid.regular_grid.dy,
36-
model.grid.regular_grid.dz))
37-
38-
mc_vertices.append(verts + [model.grid.regular_grid.extent[0],
39-
model.grid.regular_grid.extent[2],
40-
model.grid.regular_grid.extent[4]])
41-
mc_edges.append(faces)
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
4236

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+
)
4350

44-
def _get_lithology_idx(faults, model):
45-
if faults is not None:
46-
false_indices = [i for i, fault in enumerate(faults) if not fault]
47-
else:
48-
false_indices = np.arange(len(model.structural_frame.structural_groups))
49-
return false_indices
5051

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]])
5185

52-
def _extract_fault_mesh(mc_edges, mc_vertices, model):
53-
# TODO: This should also use the scalar fields probably
54-
for i in np.unique(model.solutions.raw_arrays.fault_block)[:-1]:
55-
fault_block = model.solutions.raw_arrays.fault_block.reshape(model.grid.regular_grid.resolution)
56-
verts, faces, _, _ = measure.marching_cubes(fault_block,
57-
i,
58-
spacing=(model.grid.regular_grid.dx,
59-
model.grid.regular_grid.dy,
60-
model.grid.regular_grid.dz))
61-
mc_vertices.append(verts + [model.grid.regular_grid.extent[0],
62-
model.grid.regular_grid.extent[2],
63-
model.grid.regular_grid.extent[4]])
64-
mc_edges.append(faces)
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: 5 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -7,9 +7,9 @@
77
from gempy.modules.mesh_extranction import marching_cubes
88
from gempy.optional_dependencies import require_gempy_viewer
99

10-
1110
PLOT = True
1211

12+
1313
def test_marching_cubes_implementation():
1414
model = gp.generate_example_model(ExampleModel.COMBINATION, compute_model=False)
1515

@@ -26,6 +26,7 @@ def test_marching_cubes_implementation():
2626
reset=True
2727
)
2828

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

@@ -36,28 +37,13 @@ def test_marching_cubes_implementation():
3637

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

39-
mc_edges, mc_vertices = marching_cubes.compute_marching_cubes(model)
40+
marching_cubes.set_meshes_with_marching_cubes(model)
4041

4142
if PLOT:
4243
gpv = require_gempy_viewer()
4344
gtv: gpv.GemPyToVista = gpv.plot_3d(
4445
model=model,
4546
show_data=True,
46-
image=False,
47-
show=False
47+
image=True,
48+
show=True
4849
)
49-
import pyvista as pv
50-
51-
# TODO: This opens interactive window as of now
52-
pyvista_plotter: pv.Pltter = gtv.p
53-
54-
# Add the meshes to the plot
55-
for i in range(len(mc_vertices)):
56-
pyvista_plotter.add_mesh(
57-
pv.PolyData(mc_vertices[i],
58-
np.insert(mc_edges[i], 0, 3, axis=1).ravel()),
59-
color='blue')
60-
61-
pyvista_plotter.show()
62-
63-

0 commit comments

Comments
 (0)