Skip to content

Commit 07f6981

Browse files
javohajavoha
committed
Added MC code.
1 parent a02f174 commit 07f6981

File tree

1 file changed

+140
-13
lines changed

1 file changed

+140
-13
lines changed

‎test/test_modules/test_marching_cubes.py‎

Lines changed: 140 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -6,8 +6,33 @@
66
from gempy.core.data.grid_modules import RegularGrid
77
from gempy.optional_dependencies import require_gempy_viewer
88

9+
from skimage import measure
10+
911
PLOT = True
1012

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+
1136

1237
def test_marching_cubes_implementation():
1338
model = gp.generate_example_model(ExampleModel.COMBINATION, compute_model=False)
@@ -25,7 +50,7 @@ def test_marching_cubes_implementation():
2550
reset=True
2651
)
2752

28-
model.interpolation_options.evaluation_options.mesh_extraction = False # * Not extracting the mesh with dual contouring
53+
model.interpolation_options.evaluation_options.mesh_extraction = False # * Not extracting the mesh with dual contouring
2954
gp.compute_model(model)
3055

3156
# Assert
@@ -35,19 +60,121 @@ def test_marching_cubes_implementation():
3560

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

38-
# TODO: Code to extract the mesh with marching cubes
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
39163

40164
if PLOT:
41165
gpv = require_gempy_viewer()
42-
gtv: gpv.GemPyToVista = gpv.plot_3d(model, show_data=True, image=True)
166+
# gtv: gpv.GemPyToVista = gpv.plot_3d(model, show_data=True, image=True)
43167
import pyvista as pv
44-
pyvista_plotter: pv.Plotter = gtv.p
45-
46-
foo = np.empty(0) # TODO: use the marching cubes algorithm to add mesh to the gempy plot
47-
# Add the mesh to the plot
48-
pyvista_plotter.add_mesh(
49-
foo,
50-
show_edges=True,
51-
color='white',
52-
opacity=0.5,
53-
)
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)