diff --git a/gempy_engine/API/dual_contouring/_interpolate_on_edges.py b/gempy_engine/API/dual_contouring/_interpolate_on_edges.py index 897bda3..476ad54 100644 --- a/gempy_engine/API/dual_contouring/_interpolate_on_edges.py +++ b/gempy_engine/API/dual_contouring/_interpolate_on_edges.py @@ -1,4 +1,4 @@ -from typing import List, Tuple, Optional +from typing import List, Optional import numpy as np @@ -24,10 +24,15 @@ def interpolate_on_edges_for_dual_contouring( octree_leaves: OctreeLevel, mask: Optional[np.ndarray] = None ) -> DualContouringData: - - # region define location where we need to interpolate the gradients for dual contouring - output_corners: InterpOutput = octree_leaves.outputs_corners[n_scalar_field] - intersection_xyz, valid_edges = _get_intersection_on_edges(octree_leaves, output_corners, mask) + + output: InterpOutput = octree_leaves.outputs_centers[n_scalar_field] + intersection_xyz, valid_edges = find_intersection_on_edge( + _xyz_corners=octree_leaves.grid_centers.corners_grid.values, + scalar_field_on_corners=output.exported_fields.scalar_field[output.grid.corners_grid_slice], + scalar_at_sp=output.scalar_field_at_sp, + masking=mask + ) + interpolation_input.set_temp_grid(EngineGrid(custom_grid=GenericGrid(values=intersection_xyz))) # endregion @@ -43,22 +48,7 @@ def interpolate_on_edges_for_dual_contouring( xyz_on_centers=octree_leaves.grid_centers.octree_grid.values if mask is None else octree_leaves.grid_centers.octree_grid.values[mask], dxdydz=octree_leaves.grid_centers.octree_dxdydz, exported_fields_on_edges=output_on_edges[n_scalar_field].exported_fields, - n_surfaces_to_export=output_corners.scalar_field_at_sp.shape[0], + n_surfaces_to_export=output.scalar_field_at_sp.shape[0], tree_depth=options.number_octree_levels, ) return dc_data - - -# TODO: These two functions could be moved to the module -def _get_intersection_on_edges(octree_level: OctreeLevel, output_corners: InterpOutput, - mask: Optional[np.ndarray] = None) -> Tuple[np.ndarray, np.ndarray]: - # First find xyz on edges: - intersection_xyz, valid_edges = find_intersection_on_edge( - _xyz_corners=octree_level.grid_corners.values, - scalar_field_on_corners=output_corners.exported_fields.scalar_field, - scalar_at_sp=output_corners.scalar_field_at_sp, - masking=mask - ) - return intersection_xyz, valid_edges - - diff --git a/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py b/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py index 2dd8b58..f7c0e91 100644 --- a/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py +++ b/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py @@ -97,9 +97,11 @@ def dual_contouring_multi_scalar(data_descriptor: InputDataDescriptor, interpola def _mask_generation(octree_leaves, masking_option: MeshExtractionMaskingOptions) -> np.ndarray | None: - all_scalar_fields_outputs: list[InterpOutput] = octree_leaves.outputs_corners + all_scalar_fields_outputs: list[InterpOutput] = octree_leaves.outputs_centers n_scalar_fields = len(all_scalar_fields_outputs) - grid_size = all_scalar_fields_outputs[0].grid_size + outputs_ = all_scalar_fields_outputs[0] + slice_corners = outputs_.grid.corners_grid_slice + grid_size = outputs_.cornersGrid_values.shape[0] mask_matrix = BackendTensor.t.zeros((n_scalar_fields, grid_size // 8), dtype=bool) onlap_chain_counter = 0 @@ -123,15 +125,15 @@ def _mask_generation(octree_leaves, masking_option: MeshExtractionMaskingOptions # raise NotImplementedError("Onlap is not supported yet") # return octree_leaves.outputs_corners[n_scalar_field].squeezed_mask_array.reshape((1, -1, 8)).sum(-1, bool)[0] case MeshExtractionMaskingOptions.INTERSECT, StackRelationType.ERODE: - x = all_scalar_fields_outputs[i + onlap_chain_counter].squeezed_mask_array.reshape((1, -1, 8)) + x = all_scalar_fields_outputs[i + onlap_chain_counter].squeezed_mask_array[slice_corners].reshape((1, -1, 8)) mask_matrix[i] = BackendTensor.t.sum(x, -1, bool)[0] onlap_chain_counter = 0 case MeshExtractionMaskingOptions.INTERSECT, StackRelationType.BASEMENT: - x = all_scalar_fields_outputs[i].squeezed_mask_array.reshape((1, -1, 8)) + x = all_scalar_fields_outputs[i].squeezed_mask_array[slice_corners].reshape((1, -1, 8)) mask_matrix[i] = BackendTensor.t.sum(x, -1, bool)[0] onlap_chain_counter = 0 case MeshExtractionMaskingOptions.INTERSECT, StackRelationType.ONLAP: - x = all_scalar_fields_outputs[i].squeezed_mask_array.reshape((1, -1, 8)) + x = all_scalar_fields_outputs[i].squeezed_mask_array[slice_corners].reshape((1, -1, 8)) mask_matrix[i] = BackendTensor.t.sum(x, -1, bool)[0] onlap_chain_counter += 1 case _, StackRelationType.FAULT: diff --git a/gempy_engine/API/interp_single/_octree_generation.py b/gempy_engine/API/interp_single/_octree_generation.py index c2b2067..87490b3 100644 --- a/gempy_engine/API/interp_single/_octree_generation.py +++ b/gempy_engine/API/interp_single/_octree_generation.py @@ -17,7 +17,7 @@ import numpy as np -def interpolate_on_octree(interpolation_input: InterpolationInput, options: InterpolationOptions, +def interpolate_on_octree_(interpolation_input: InterpolationInput, options: InterpolationOptions, data_shape: InputDataDescriptor) -> OctreeLevel: if BackendTensor.engine_backend is not AvailableBackends.PYTORCH and NOT_MAKE_INPUT_DEEP_COPY is False: temp_interpolation_input = copy.deepcopy(interpolation_input) @@ -30,8 +30,9 @@ def interpolate_on_octree(interpolation_input: InterpolationInput, options: Inte # * Interpolate - corners grid_0_centers: EngineGrid = temp_interpolation_input.grid # ? This could be moved to the next section if options.compute_corners: - grid_0_corners: Optional[EngineGrid] = EngineGrid.from_xyz_coords( - xyz_coords=_generate_corners(regular_grid=grid_0_centers.octree_grid) + xyz_corners = _generate_corners(regular_grid=grid_0_centers.octree_grid) + grid_0_corners: EngineGrid = EngineGrid.from_xyz_coords( + xyz_coords=xyz_corners ) # ! Here we need to swap the grid temporarily but it is @@ -58,6 +59,51 @@ def interpolate_on_octree(interpolation_input: InterpolationInput, options: Inte return next_octree_level +def interpolate_on_octree(interpolation_input: InterpolationInput, options: InterpolationOptions, + data_shape: InputDataDescriptor) -> OctreeLevel: + if BackendTensor.engine_backend is not AvailableBackends.PYTORCH and NOT_MAKE_INPUT_DEEP_COPY is False: + temp_interpolation_input = copy.deepcopy(interpolation_input) + else: + temp_interpolation_input = interpolation_input + + # * Interpolate - corners + if options.compute_corners: + grid_0_centers: EngineGrid = temp_interpolation_input.grid # ? This could be moved to the next section + xyz_corners = _generate_corners(regular_grid=grid_0_centers.octree_grid) + + corner_grid = GenericGrid(values=xyz_corners) + grid_0_centers.corners_grid = corner_grid + output_0_centers: List[InterpOutput] = interpolate_all_fields(temp_interpolation_input, options, data_shape) # interpolate - centers + + # * DEP + grid_0_corners = None + output_0_corners = [] + + # * Create next octree level + next_octree_level = OctreeLevel( + grid_centers=grid_0_centers, + grid_corners=grid_0_corners, + outputs_centers=output_0_centers, + outputs_corners=output_0_corners + ) + else: + grid_0_centers: EngineGrid = temp_interpolation_input.grid # ? This could be moved to the next section + output_0_centers: List[InterpOutput] = interpolate_all_fields(temp_interpolation_input, options, data_shape) # interpolate - centers + + # * DEP + output_0_corners = [] + grid_0_corners = None + + # * Create next octree level + next_octree_level = OctreeLevel( + grid_centers=grid_0_centers, + grid_corners=grid_0_corners, + outputs_centers=output_0_centers, + outputs_corners=output_0_corners + ) + + return next_octree_level + def _generate_corners_DEP(regular_grid: RegularGrid, level=1) -> np.ndarray: if regular_grid is None: raise ValueError("Regular grid is None") diff --git a/gempy_engine/core/data/engine_grid.py b/gempy_engine/core/data/engine_grid.py index fe5bd61..2d03712 100644 --- a/gempy_engine/core/data/engine_grid.py +++ b/gempy_engine/core/data/engine_grid.py @@ -22,6 +22,7 @@ class EngineGrid: topography: Optional[GenericGrid] = None sections: Optional[GenericGrid] = None geophysics_grid: Optional[CenteredGrid] = None # TODO: Not implemented this probably will need something different that the generic grid? + corners_grid: Optional[GenericGrid] = None # TODO: Not implemented this probably will need something different that the generic grid? debug_vals = None @@ -29,13 +30,15 @@ class EngineGrid: def __init__(self, octree_grid: Optional[RegularGrid] = None, dense_grid: Optional[RegularGrid] = None, custom_grid: Optional[GenericGrid] = None, topography: Optional[GenericGrid] = None, - sections: Optional[GenericGrid] = None, geophysics_grid: Optional[CenteredGrid] = None): + sections: Optional[GenericGrid] = None, geophysics_grid: Optional[CenteredGrid] = None, + corners_grid: Optional[GenericGrid] = None): self.octree_grid = octree_grid self.dense_grid = dense_grid self.custom_grid = custom_grid self.topography = topography self.sections = sections self.geophysics_grid = geophysics_grid + self.corners_grid = corners_grid @property def regular_grid(self): @@ -78,6 +81,8 @@ def values(self) -> np.ndarray: values.append(self.sections.values) if self.geophysics_grid is not None: values.append(self.geophysics_grid.values) + if self.corners_grid is not None: + values.append(self.corners_grid.values) values_array = BackendTensor.t.concatenate(values, dtype=BackendTensor.dtype) values_array = BackendTensor.t.array(values_array, dtype=BackendTensor.dtype) @@ -131,6 +136,14 @@ def geophysics_grid_slice(self) -> slice: start + len(self.geophysics_grid) if self.geophysics_grid is not None else start ) + @property + def corners_grid_slice(self) -> slice: + start = self.geophysics_grid_slice.stop + return slice( + start, + start + len(self.corners_grid) if self.corners_grid is not None else start + ) + @property def len_all_grids(self) -> int: return self.values.shape[0] diff --git a/gempy_engine/core/data/interp_output.py b/gempy_engine/core/data/interp_output.py index a5cc86e..e795ca2 100644 --- a/gempy_engine/core/data/interp_output.py +++ b/gempy_engine/core/data/interp_output.py @@ -84,6 +84,14 @@ def custom_grid_values(self): def geophysics_grid_values(self): return self.block[self.grid.geophysics_grid_slice] + @property + def cornersGrid_values(self): + return self.block[self.grid.corners_grid_slice] + + @property + def ids_cornersGrid(self): + return BackendTensor.t.rint(self.block[self.grid.corners_grid_slice]) + @property def ids_geophysics_grid(self): return BackendTensor.t.rint(self.block[self.grid.geophysics_grid_slice]) @@ -132,6 +140,21 @@ def litho_faults_ids(self): # Generate the unique IDs unique_ids = litho_ids + faults_ids * multiplier return unique_ids + + @property + def litho_faults_ids_corners_grid(self): + if self.combined_scalar_field is None: # * This in principle is only used for testing + return self.ids_cornersGrid + + litho_ids = BackendTensor.t.rint(self.block[self.grid.corners_grid_slice]) + faults_ids = BackendTensor.t.rint(self.faults_block[self.grid.corners_grid_slice]) + + # Get the number of unique lithology IDs + multiplier = len(BackendTensor.t.unique(litho_ids)) + + # Generate the unique IDs + unique_ids = litho_ids + faults_ids * multiplier + return unique_ids def get_block_from_value_type(self, value_type: ValueType, slice_: slice): match value_type: diff --git a/gempy_engine/core/data/octree_level.py b/gempy_engine/core/data/octree_level.py index 8103ede..e1c3748 100644 --- a/gempy_engine/core/data/octree_level.py +++ b/gempy_engine/core/data/octree_level.py @@ -48,7 +48,11 @@ def output_corners(self): # * Alias @property def last_output_corners(self): return self.outputs_corners[-1] - + + @property + def litho_faults_ids_corners_grid(self): + return self.outputs_centers[-1].litho_faults_ids_corners_grid + @property def number_of_outputs(self): return len(self.outputs_centers) \ No newline at end of file diff --git a/gempy_engine/core/data/regular_grid.py b/gempy_engine/core/data/regular_grid.py index bea2e48..f6326b9 100644 --- a/gempy_engine/core/data/regular_grid.py +++ b/gempy_engine/core/data/regular_grid.py @@ -20,7 +20,8 @@ class RegularGrid: original_values: np.ndarray = field(default=None, repr=False, init=False) #: When the regular grid is representing a octree level, only active cells are stored in values. This is the original values of the regular grid. def __len__(self): - return self.regular_grid_shape.prod() + # return self.regular_grid_shape.prod() + return self.values.shape[0] def __post_init__(self): self.regular_grid_shape = _check_and_convert_list_to_array(self.regular_grid_shape) diff --git a/gempy_engine/modules/octrees_topology/_octree_internals.py b/gempy_engine/modules/octrees_topology/_octree_internals.py index 7178031..8e6d315 100644 --- a/gempy_engine/modules/octrees_topology/_octree_internals.py +++ b/gempy_engine/modules/octrees_topology/_octree_internals.py @@ -14,7 +14,7 @@ def compute_next_octree_locations(prev_octree: OctreeLevel, evaluation_options: EvaluationOptions, current_octree_level: int) -> EngineGrid: - ids = prev_octree.last_output_corners.litho_faults_ids + ids = prev_octree.litho_faults_ids_corners_grid uv_8 = ids.reshape((-1, 8)) # Old octree diff --git a/tests/test_common/test_integrations/test_multi_fields.py b/tests/test_common/test_integrations/test_multi_fields.py index d31e5b5..5ecd6de 100644 --- a/tests/test_common/test_integrations/test_multi_fields.py +++ b/tests/test_common/test_integrations/test_multi_fields.py @@ -147,9 +147,7 @@ def test_plot_corners(unconformity_complex, n_oct_levels=2): interpolation_input, options, structure = unconformity_complex options.number_octree_levels = n_oct_levels solutions: Solutions = compute_model(interpolation_input, options, structure) - output_corners: InterpOutput = solutions.octrees_output[-1].outputs_corners[-1] - - vertices = output_corners.grid.values + vertices = solutions.octrees_output[-1].grid_centers.corners_grid.values if plot_pyvista or False: helper_functions_pyvista.plot_pyvista(solutions.octrees_output, v_just_points=vertices) diff --git a/tests/test_common/test_modules/test_dual.py b/tests/test_common/test_modules/test_dual.py index eedbfdd..586e254 100644 --- a/tests/test_common/test_modules/test_dual.py +++ b/tests/test_common/test_modules/test_dual.py @@ -1,33 +1,25 @@ import copy -from typing import List - +import numpy as np +import os import pytest +from typing import List -from gempy_engine.API.dual_contouring._interpolate_on_edges import _get_intersection_on_edges -from gempy_engine.API.interp_single._interp_scalar_field import _evaluate_sys_eq -from gempy_engine.API.interp_single._interp_single_feature import input_preprocess +import gempy_engine.API.interp_single.interp_features as interp +from gempy_engine.API.dual_contouring._dual_contouring import compute_dual_contouring from gempy_engine.API.interp_single._multi_scalar_field_manager import interpolate_all_fields +from gempy_engine.API.interp_single.interp_features import interpolate_n_octree_levels, interpolate_and_segment from gempy_engine.API.model.model_api import compute_model from gempy_engine.config import AvailableBackends from gempy_engine.core.backend_tensor import BackendTensor +from gempy_engine.core.data import InterpolationOptions +from gempy_engine.core.data.dual_contouring_data import DualContouringData +from gempy_engine.core.data.dual_contouring_mesh import DualContouringMesh from gempy_engine.core.data.engine_grid import EngineGrid, RegularGrid -import numpy as np -import os - -import gempy_engine.API.interp_single.interp_features as interp - -from gempy_engine.API.dual_contouring._dual_contouring import compute_dual_contouring from gempy_engine.core.data.input_data_descriptor import InputDataDescriptor -from gempy_engine.modules.activator.activator_interface import activate_formation_block -from gempy_engine.core.data.internal_structs import SolverInput - -from gempy_engine.core.data.solutions import Solutions -from gempy_engine.core.data.dual_contouring_mesh import DualContouringMesh -from gempy_engine.core.data.dual_contouring_data import DualContouringData -from gempy_engine.core.data.octree_level import OctreeLevel from gempy_engine.core.data.interp_output import InterpOutput from gempy_engine.core.data.interpolation_input import InterpolationInput -from gempy_engine.API.interp_single.interp_features import interpolate_n_octree_levels, interpolate_and_segment +from gempy_engine.core.data.octree_level import OctreeLevel +from gempy_engine.core.data.solutions import Solutions from gempy_engine.modules.dual_contouring.dual_contouring_interface import QEF, find_intersection_on_edge, triangulate_dual_contouring from gempy_engine.modules.octrees_topology.octrees_topology_interface import get_regular_grid_value_for_level from gempy_engine.plugins.plotting import helper_functions_pyvista @@ -42,6 +34,18 @@ plot_pyvista = False +def _grab_xyz_edges(last_octree_level: OctreeLevel) -> tuple: + corners = last_octree_level.outputs_centers[0] + # First find xyz on edges: + xyz, edges = find_intersection_on_edge( + _xyz_corners=last_octree_level.grid_centers.corners_grid.values, + scalar_field_on_corners=corners.exported_fields.scalar_field[corners.grid.corners_grid_slice], + scalar_at_sp=corners.scalar_field_at_sp, + masking=None + ) + return xyz, edges + + def test_compute_dual_contouring_api(simple_model, simple_grid_3d_octree): # region Test find_intersection_on_edge spi, ori_i, options, data_shape = simple_model @@ -56,19 +60,12 @@ def test_compute_dual_contouring_api(simple_model, simple_grid_3d_octree): last_octree_level: OctreeLevel = octree_list[-1] - intersection_xyz, valid_edges = _get_intersection_on_edges(last_octree_level, last_octree_level.outputs_corners[0]) - interpolation_input.set_temp_grid(EngineGrid.from_xyz_coords(intersection_xyz)) - - - output_on_edges: List[InterpOutput] = interpolate_all_fields(interpolation_input, options, data_shape) - - dc_data = DualContouringData( - xyz_on_edge=intersection_xyz, - valid_edges=valid_edges, - xyz_on_centers=last_octree_level.grid_centers.values, - dxdydz=last_octree_level.grid_centers.octree_dxdydz, - exported_fields_on_edges=output_on_edges[0].exported_fields, - n_surfaces_to_export=data_shape.tensors_structure.n_surfaces + intersection_xyz, valid_edges = _grab_xyz_edges(last_octree_level) + dc_data = _gen_dc_data( + octree_level=last_octree_level, + interpolation_input=interpolation_input, + data_shape=data_shape, + options=options, ) gradients = dc_data.gradients @@ -96,6 +93,25 @@ def test_compute_dual_contouring_api(simple_model, simple_grid_3d_octree): # endregion +def _gen_dc_data(octree_level: OctreeLevel, interpolation_input: InterpolationInput, + options: InterpolationOptions, data_shape: InputDataDescriptor) -> DualContouringData: + intersection_xyz, valid_edges = _grab_xyz_edges(octree_level) + + interpolation_input.set_temp_grid(EngineGrid.from_xyz_coords(intersection_xyz)) + + output_on_edges: List[InterpOutput] = interpolate_all_fields(interpolation_input, options, data_shape) + + dc_data = DualContouringData( + xyz_on_edge=intersection_xyz, + valid_edges=valid_edges, + xyz_on_centers=octree_level.grid_centers.octree_grid.values, + dxdydz=octree_level.grid_centers.octree_dxdydz, + exported_fields_on_edges=output_on_edges[0].exported_fields, + n_surfaces_to_export=data_shape.tensors_structure.n_surfaces + ) + return dc_data + + @pytest.mark.skipif(BackendTensor.engine_backend != AvailableBackends.numpy, reason="Only numpy supported") def test_compute_mesh_extraction_fancy_triangulation(simple_model, simple_grid_3d_octree): from gempy_engine.modules.dual_contouring.fancy_triangulation import get_left_right_array, triangulate @@ -124,23 +140,43 @@ def _simple_grid_3d_octree_regular(): octree_level_for_surface: OctreeLevel = octree_list[options.number_octree_levels_surface - 1] - intersection_xyz, valid_edges = _get_intersection_on_edges( + # corners = octree_level_for_surface.outputs_corners[0] + # # First find xyz on edges: + # xyz, edges = find_intersection_on_edge( + # _xyz_corners=octree_level_for_surface.grid_corners.values, + # scalar_field_on_corners=corners.exported_fields.scalar_field, + # scalar_at_sp=corners.scalar_field_at_sp, + # masking=None + # ) + # result = xyz, edges + # intersection_xyz, valid_edges = result + # + # last_octree_level: OctreeLevel = octree_list[-1] + # intersection_xyz, valid_edges = _grab_xyz_edges(octree_level_for_surface) + # interpolation_input.set_temp_grid(EngineGrid.from_xyz_coords(intersection_xyz)) + # + # output_on_edges: List[InterpOutput] = interpolate_all_fields(interpolation_input, options, data_shape) + # + # dc_data = DualContouringData( + # xyz_on_edge=intersection_xyz, + # valid_edges=valid_edges, + # xyz_on_centers=octree_level_for_surface.grid_centers.values, + # dxdydz=octree_level_for_surface.grid_centers.octree_dxdydz, + # exported_fields_on_edges=output_on_edges[0].exported_fields, + # n_surfaces_to_export=data_shape.tensors_structure.n_surfaces + # ) + + # dc_data = _gen_dc_data(data_shape, interpolation_input, intersection_xyz, last_octree_level, options, valid_edges) + + dc_data = _gen_dc_data( octree_level=octree_level_for_surface, - output_corners=octree_level_for_surface.outputs_corners[0] + interpolation_input=interpolation_input, + data_shape=data_shape, + options=options, ) - - interpolation_input.set_temp_grid(EngineGrid.from_xyz_coords(intersection_xyz)) - output_on_edges: List[InterpOutput] = interpolate_all_fields(interpolation_input, options, data_shape) - - dc_data = DualContouringData( - xyz_on_edge=intersection_xyz, - valid_edges=valid_edges, - xyz_on_centers=octree_level_for_surface.grid_centers.values, - dxdydz=octree_level_for_surface.grid_centers.octree_dxdydz, - exported_fields_on_edges=output_on_edges[0].exported_fields, - n_surfaces_to_export=data_shape.tensors_structure.n_surfaces - ) - + + intersection_xyz, valid_edges = _grab_xyz_edges(octree_level_for_surface) + dc_meshes: List[DualContouringMesh] = compute_dual_contouring(dc_data) dc_data = dc_meshes[0].dc_data valid_voxels = dc_data.valid_voxels @@ -212,7 +248,7 @@ def test_compute_dual_contouring_complex(unconformity_complex_one_layer, n_oct_l dc_data = solutions.dc_meshes[0].dc_data if plot_pyvista or False: - output_corners: InterpOutput = solutions.octrees_output[-1].outputs_corners[-1] + output_corners: InterpOutput = solutions.octrees_output[-1].outputs_centers[-1] vertices = output_corners.grid.values intersection_xyz = dc_data.xyz_on_edge @@ -245,7 +281,16 @@ def test_compute_dual_contouring_several_meshes(simple_model_3_layers, simple_gr last_octree_level: OctreeLevel = octree_list[-1] - intersection_xyz, valid_edges = _get_intersection_on_edges(last_octree_level, last_octree_level.outputs_corners[0]) + corners = last_octree_level.outputs_centers[0] + # First find xyz on edges: + xyz, edges = find_intersection_on_edge( + _xyz_corners=last_octree_level.grid_centers.corners_grid.values, + scalar_field_on_corners=corners.exported_fields.scalar_field[corners.grid.corners_grid_slice], + scalar_at_sp=corners.scalar_field_at_sp, + masking=None + ) + result = xyz, edges + intersection_xyz, valid_edges = result interpolation_input.set_temp_grid(EngineGrid.from_xyz_coords(intersection_xyz)) output_on_edges = interp.interpolate_single_field(interpolation_input, options, data_shape.tensors_structure) @@ -253,7 +298,7 @@ def test_compute_dual_contouring_several_meshes(simple_model_3_layers, simple_gr dc_data = DualContouringData( xyz_on_edge=intersection_xyz, valid_edges=valid_edges, - xyz_on_centers=last_octree_level.grid_centers.values, + xyz_on_centers=last_octree_level.grid_centers.octree_grid.values, dxdydz=last_octree_level.grid_centers.octree_dxdydz, exported_fields_on_edges=output_on_edges.exported_fields, n_surfaces_to_export=data_shape.tensors_structure.n_surfaces @@ -279,24 +324,24 @@ def test_compute_dual_contouring_several_meshes(simple_model_3_layers, simple_gr @pytest.mark.skipif(BackendTensor.engine_backend != AvailableBackends.numpy, reason="Only numpy supported") -def test_find_edges_intersection_step_by_step(simple_model, simple_grid_3d_octree): +def test_find_edges_intersection_pro(simple_model, simple_grid_3d_octree): # region Test find_intersection_on_edge spi, ori_i, options, data_shape = simple_model ids = np.array([1, 2]) - grid_0_centers = copy.deepcopy(simple_grid_3d_octree) + grid_0_centers = simple_grid_3d_octree interpolation_input = InterpolationInput(spi, ori_i, grid_0_centers, ids) options.number_octree_levels = 5 options.number_octree_levels_surface = 5 options.compute_scalar_gradient = True - + octree_list = interpolate_n_octree_levels(interpolation_input, options, data_shape) last_octree_level: OctreeLevel = octree_list[-1] - sfsp = last_octree_level.last_output_corners.scalar_field_at_sp + sfsp = last_octree_level.last_output_center.scalar_field_at_sp - xyz_on_edge, valid_edges = find_intersection_on_edge(last_octree_level.grid_corners.values, last_octree_level.output_corners.exported_fields.scalar_field, sfsp, ) + xyz_on_edge, valid_edges = find_intersection_on_edge(last_octree_level.grid_centers.corners_grid.values, last_octree_level.last_output_center.exported_fields.scalar_field[last_octree_level.last_output_center.grid.corners_grid_slice], sfsp, ) # endregion @@ -395,111 +440,6 @@ def test_find_edges_intersection_step_by_step(simple_model, simple_grid_3d_octre plot_label=False, plot_marching_cubes=False) -@pytest.mark.skipif(BackendTensor.engine_backend != AvailableBackends.numpy, reason="Only numpy supported") -def test_find_edges_intersection_pro(simple_model, simple_grid_3d_octree): - # region Test find_intersection_on_edge - spi, ori_i, options, data_shape = simple_model - ids = np.array([1, 2]) - grid_0_centers = simple_grid_3d_octree - interpolation_input = InterpolationInput(spi, ori_i, grid_0_centers, ids) - - options.compute_scalar_gradient = True - octree_list = interpolate_n_octree_levels(interpolation_input, options, data_shape) - - last_octree_level: OctreeLevel = octree_list[-1] - - sfsp = last_octree_level.output_corners.scalar_field_at_sp - # sfsp = np.append(sfsp, -0.1) - xyz_on_edge, valid_edges = find_intersection_on_edge(last_octree_level.grid_corners.values, last_octree_level.output_corners.exported_fields.scalar_field, sfsp, ) - # endregion - - # region Get Normals - interpolation_input.set_temp_grid(EngineGrid.from_xyz_coords(xyz_on_edge)) - - output_on_edges = interpolate_and_segment(interpolation_input, options, data_shape.tensors_structure) - # stack gradients output_on_edges.exported_fields.gx_field - gradients = np.stack( - (output_on_edges.exported_fields.gx_field, - output_on_edges.exported_fields.gy_field, - output_on_edges.exported_fields.gz_field), axis=0).T - - # endregion - - # region Prepare data for vectorized QEF - - n_edges = valid_edges.shape[0] - - # Coordinates for all posible edges (12) and 3 dummy normals in the center - xyz = np.zeros((n_edges, 15, 3)) - normals = np.zeros((n_edges, 15, 3)) - - xyz[:, :12][valid_edges] = xyz_on_edge - normals[:, :12][valid_edges] = gradients - - BIAS_STRENGTH = 0.1 - - xyz_aux = np.copy(xyz[:, :12]) - - # Numpy zero values to nans - xyz_aux[np.isclose(xyz_aux, 0)] = np.nan - # Mean ignoring nans - mass_points = np.nanmean(xyz_aux, axis=1) - - xyz[:, 12] = mass_points - xyz[:, 13] = mass_points - xyz[:, 14] = mass_points - - normals[:, 12] = np.array([BIAS_STRENGTH, 0, 0]) - normals[:, 13] = np.array([0, BIAS_STRENGTH, 0]) - normals[:, 14] = np.array([0, 0, BIAS_STRENGTH]) - - # Remove unused voxels - bo = valid_edges.sum(axis=1, dtype=bool) - xyz = xyz[bo] - normals = normals[bo] - - # Compute LSTSQS in all voxels at the same time - A1 = normals - b1 = xyz - bb1 = (A1 * b1).sum(axis=2) - s1 = np.einsum("ijk, ilj->ikl", A1, np.transpose(A1, (0, 2, 1))) - s2 = np.linalg.inv(s1) - s3 = np.einsum("ijk,ik->ij", np.transpose(A1, (0, 2, 1)), bb1) - v_pro = np.einsum("ijk, ij->ik", s2, s3) - - # endregion - - # endregion - - # region triangulate - grid_centers = last_octree_level.grid_centers - valid_voxels = valid_edges.sum(axis=1, dtype=bool) - - temp_ids = octree_list[-1].last_output_center.ids_block # ! I need this because setters in python sucks - temp_ids[valid_voxels] = 5 - octree_list[-1].last_output_center.ids_block = temp_ids # paint valid voxels - - dc_data = DualContouringData( - xyz_on_edge=xyz_on_edge, - xyz_on_centers=grid_centers.values, - dxdydz=grid_centers.octree_dxdydz, - valid_edges=valid_edges, - exported_fields_on_edges=None, - n_surfaces_to_export=data_shape.tensors_structure.n_surfaces - ) - - indices = triangulate_dual_contouring(dc_data) - # endregion - - if plot_pyvista or False: - # ! I leave this test for the assert as comparison to the other implementation. The model looks bad - # ! with this level of BIAS - center_mass = xyz[:, 12:].reshape(-1, 3) - normals = normals[:, 12:].reshape(-1, 3) - _plot_pyvista(last_octree_level, octree_list, simple_model, ids, grid_0_centers, - xyz_on_edge, gradients, a=center_mass, b=normals, - v_pro=v_pro, indices=indices, plot_marching_cubes=True - ) @pytest.mark.skipif(BackendTensor.engine_backend != AvailableBackends.numpy, reason="Only numpy supported") @@ -518,8 +458,8 @@ def test_find_edges_intersection_bias_on_center_of_the_cell(simple_model, simple last_octree_level: OctreeLevel = octree_list[-1] - sfsp = last_octree_level.output_corners.scalar_field_at_sp - xyz_on_edge, valid_edges = find_intersection_on_edge(last_octree_level.grid_corners.values, last_octree_level.output_corners.exported_fields.scalar_field, sfsp, ) + sfsp = last_octree_level.last_output_center.scalar_field_at_sp + xyz_on_edge, valid_edges = find_intersection_on_edge(last_octree_level.grid_centers.corners_grid.values, last_octree_level.last_output_center.exported_fields.scalar_field[last_octree_level.last_output_center.grid.corners_grid_slice], sfsp, ) valid_voxels = valid_edges.sum(axis=1, dtype=bool) # endregion @@ -549,7 +489,7 @@ def test_find_edges_intersection_bias_on_center_of_the_cell(simple_model, simple BIAS_STRENGTH = 1 - mass_points = last_octree_level.grid_centers.values + mass_points = last_octree_level.grid_centers.octree_grid.values xyz[:, 12] = mass_points xyz[:, 13] = mass_points @@ -580,13 +520,13 @@ def test_find_edges_intersection_bias_on_center_of_the_cell(simple_model, simple # region triangulate grid_centers = last_octree_level.grid_centers - temp_ids = octree_list[-1].output_centers.ids_block # ! I need this because setters in python sucks + temp_ids = octree_list[-1].last_output_center.ids_block # ! I need this because setters in python sucks temp_ids[valid_voxels] = 5 - octree_list[-1].output_centers.ids_block = temp_ids # paint valid voxels + octree_list[-1].last_output_center.ids_block = temp_ids # paint valid voxels dc_data = DualContouringData( xyz_on_edge=xyz_on_edge, - xyz_on_centers=grid_centers.values, + xyz_on_centers=grid_centers.octree_grid.values, dxdydz=grid_centers.octree_dxdydz, valid_edges=valid_edges, exported_fields_on_edges=None, @@ -614,15 +554,6 @@ def _plot_pyvista(last_octree_level, octree_list, simple_model, ids, grid_0_cent ): p = pv.Plotter() - # Plot Actual mesh (from marching cubes) - if plot_marching_cubes: - output_1_centers = last_octree_level.output_centers - resolution = [20, 20, 20] - mesh = _compute_actual_mesh(simple_model, ids, grid_0_centers, resolution, - output_1_centers.scalar_field_at_sp, output_1_centers.weights) - - p.add_mesh(mesh, opacity=.8, silhouette=True) - # Plot Regular grid Octree regular_grid_values = octree_list[n].grid_centers.octree_grid.values_vtk_format regular_grid_scalar = get_regular_grid_value_for_level(octree_list, n) @@ -666,40 +597,3 @@ def _plot_pyvista(last_octree_level, octree_list, simple_model, ids, grid_0_cent p.add_axes() p.show() - -def _compute_actual_mesh(simple_model, ids, grid, resolution, scalar_at_surface_points, weights): - raise NotImplementedError("This is broken since a while") - interpolation_input = simple_model[0] - options = simple_model[1] - shape: InputDataDescriptor = simple_model[2] - - from gempy_engine.core.data.engine_grid import EngineGrid, RegularGrid - - # region interpolate high res grid - grid_high_res:EngineGrid = EngineGrid.from_regular_grid(RegularGrid([0.25, .75, 0.25, .75, 0.25, .75], resolution)) - interpolation_input.grid = grid_high_res - input1: SolverInput = input_preprocess(shape.tensors_structure, interpolation_input) - exported_fields_high_res = _evaluate_sys_eq(input1, weights, options) - - exported_fields_high_res.set_structure_values( - reference_sp_position=shape.tensors_structure.reference_sp_position, - slice_feature=interpolation_input.slice_feature, - grid_size=interpolation_input.grid.len_all_grids) - - res = activate_formation_block(exported_fields_high_res, ids, sigmoid_slope=50000) - result = res, exported_fields_high_res, grid_high_res.octree_dxdydz - values_block_high_res, scalar_high_res, dxdydz = result - # endregion - - from skimage.measure import marching_cubes - import pyvista as pv - spacing = np.array(dxdydz)/8 - vert, edges, _, _ = marching_cubes( - volume=(scalar_high_res.scalar_field[:-8].reshape(resolution)), - level=scalar_at_surface_points[0], - spacing=spacing - ) - loc_0 = np.array([0.25, .25, .25]) + np.array(spacing) / 2 - vert += np.array(loc_0).reshape(1, 3) - mesh = pv.PolyData(vert, np.insert(edges, 0, 3, axis=1).ravel()) - return mesh diff --git a/tests/test_common/test_modules/test_octrees.py b/tests/test_common/test_modules/test_octrees.py index 2a3e5d6..c120dcf 100644 --- a/tests/test_common/test_modules/test_octrees.py +++ b/tests/test_common/test_modules/test_octrees.py @@ -4,13 +4,9 @@ import pytest import gempy_engine.API.interp_single.interp_features as interp -from gempy_engine.API.interp_single._octree_generation import _generate_corners -from gempy_engine.core.data.engine_grid import EngineGrid from gempy_engine.core.data.interpolation_input import InterpolationInput -from gempy_engine.core.data.octree_level import OctreeLevel from gempy_engine.modules.octrees_topology.octrees_topology_interface import get_next_octree_grid, \ get_regular_grid_value_for_level -from .test_dual import _compute_actual_mesh from ...conftest import plot_pyvista, TEST_SPEED dir_name = os.path.dirname(__file__)