Skip to content

Commit 3e4614e

Browse files
LeguarkLeguark
authored
[ENH/WIP] Adding first test to implement marching cubes (#1000)
# Description Implementing marching cubes Relates to <issue> # Checklist - [ ] My code uses type hinting for function and method arguments and return values. - [ ] I have created tests which cover my code. - [ ] 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 258f5ba + 7cba374 commit 3e4614e

File tree

2 files changed

+185
-1
lines changed

2 files changed

+185
-1
lines changed

Diff for: ‎gempy/core/data/grid.py

+5-1
Original file line numberDiff line numberDiff line change
@@ -254,7 +254,11 @@ def _update_values(self):
254254
if self.GridTypes.CENTERED in self.active_grids:
255255
if self.centered_grid is None: raise AttributeError('Centered grid is active but not defined')
256256
values.append(self.centered_grid.values)
257-
257+
258+
# make sure values is not empty
259+
if len(values) == 0:
260+
return self.values
261+
258262
self.values = np.concatenate(values)
259263

260264
return self.values

Diff for: ‎test/test_modules/test_marching_cubes.py

+180
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,180 @@
1+
import numpy as np
2+
from gempy_engine.core.data.raw_arrays_solution import RawArraysSolution
3+
4+
import gempy as gp
5+
from gempy.core.data.enumerators import ExampleModel
6+
from gempy.core.data.grid_modules import RegularGrid
7+
from gempy.optional_dependencies import require_gempy_viewer
8+
9+
from skimage import measure
10+
11+
PLOT = True
12+
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+
36+
37+
def test_marching_cubes_implementation():
38+
model = gp.generate_example_model(ExampleModel.COMBINATION, compute_model=False)
39+
40+
# Change the grid to only be the dense grid
41+
dense_grid: RegularGrid = RegularGrid(
42+
extent=model.grid.extent,
43+
resolution=np.array([20, 20, 20])
44+
)
45+
46+
model.grid.dense_grid = dense_grid
47+
gp.set_active_grid(
48+
grid=model.grid,
49+
grid_type=[model.grid.GridTypes.DENSE],
50+
reset=True
51+
)
52+
53+
model.interpolation_options.evaluation_options.mesh_extraction = False # * Not extracting the mesh with dual contouring
54+
gp.compute_model(model)
55+
56+
# Assert
57+
assert model.solutions.block_solution_type == RawArraysSolution.BlockSolutionType.DENSE_GRID
58+
assert model.solutions.dc_meshes is None
59+
arrays = model.solutions.raw_arrays # * arrays is equivalent to gempy v2 solutions
60+
61+
assert arrays.scalar_field_matrix.shape == (3, 8_000) # * 3 surfaces, 8000 points
62+
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
163+
164+
if PLOT:
165+
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()

0 commit comments

Comments
 (0)