Skip to content

Commit d4b4277

Browse files
LeguarkLeguark
committed
[ENH/CLN/TEST/WIP] Moving marching cubes code to a module
1 parent 76756ad commit d4b4277

File tree

2 files changed

+75
-128
lines changed

2 files changed

+75
-128
lines changed

‎gempy/modules/mesh_extranction/marching_cubes.py

Lines changed: 64 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,64 @@
1+
import numpy as np
2+
from skimage import measure
3+
4+
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
16+
17+
# Extract scalar field values for elements
18+
scalar_values = model.solutions.raw_arrays.scalar_field_at_surface_points
19+
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)
42+
43+
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
50+
51+
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)

‎test/test_modules/test_marching_cubes.py

Lines changed: 11 additions & 128 deletions
Original file line numberDiff line numberDiff line change
@@ -4,36 +4,12 @@
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
1010

1111
PLOT = True
1212

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-
3713
def test_marching_cubes_implementation():
3814
model = gp.generate_example_model(ExampleModel.COMBINATION, compute_model=False)
3915

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

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

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
39+
mc_edges, mc_vertices = marching_cubes.compute_marching_cubes(model)
16340

16441
if PLOT:
16542
gpv = require_gempy_viewer()
166-
# gtv: gpv.GemPyToVista = gpv.plot_3d(model, show_data=True, image=True)
43+
gtv: gpv.GemPyToVista = gpv.plot_3d(
44+
model=model,
45+
show_data=True,
46+
image=False,
47+
show=False
48+
)
16749
import pyvista as pv
168-
# pyvista_plotter: pv.Plotter = gtv.p
16950

17051
# TODO: This opens interactive window as of now
171-
pyvista_plotter = pv.Plotter()
52+
pyvista_plotter: pv.Pltter = gtv.p
17253

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

18061
pyvista_plotter.show()
62+
63+

0 commit comments

Comments
 (0)