diff --git a/docs/extent_capping.md b/docs/extent_capping.md new file mode 100644 index 00000000..5a7d9860 --- /dev/null +++ b/docs/extent_capping.md @@ -0,0 +1,114 @@ +# Best-effort extent capping + +Extent capping is opt-in and operates independently on each exported isosurface: + +```python +from gempy_engine.core.data import MeshExtentCapping + +options.evaluation_options.mesh_extraction_extent_capping = ( + MeshExtentCapping.SCALAR_LESS_EQUAL +) +``` + +The equivalent string is `"scalar_less_equal"`; the default is `"none"`. +The inside convention is `scalar <= isovalue`. Nested surfaces therefore produce +overlapping solids, not a partition into geological unit volumes. + +## Geometry + +The capper constructs three connected pieces: + +1. The ordinary interior dual-contouring triangles. +2. Transition triangles from boundary-cell QEF vertices to the extent contour, + including primal-edge triangles between adjacent boundary cells. +3. Flat, outward-facing triangles covering inside regions of the six box faces. + +Boundary squares use a consistent diagonal and piecewise-linear triangle +clipping. Lattice points and edge intersections have shared integer identities, +including intersections on face diagonals and on box edges and corners. + +Only cap components connected to a contour with an available, contained QEF +vertex are retained. An enclosed isosurface receives no disconnected box shell. +If no isosurface intersects the box, no box shell is manufactured, even when +the whole box satisfies the inside condition. This is closure of existing +isosurfaces, not unconditional extraction of the complete clipped sublevel set. + +`RegularGrid.orthogonal_extent` is the single requested extent used for grid +sampling, refinement, and capping. All octree levels retain that exact extent, +whether capping is enabled or disabled. There is no conditional private root-grid +resampling or extent translation. Caps are generated in engine coordinates, +before downstream output transforms. Original vertex indices are preserved and +cap vertices are appended. + +Enabled extraction and cap clipping share strict scalar crossing rules. An +endpoint equal to the isovalue is inside; an entire iso-valued edge is not a +crossing. Interpolation parameters are clamped and snapped within `1e-12` of an +endpoint. Nonfinite scalar inputs are rejected with a diagnostic. QEF mass +points use edge-validity masks globally, including genuine zero coordinates, +regardless of capping mode. Enabled QEF solving avoids PyTorch's additional +origin-centered regularization. + +With capping disabled, no boundary evaluation or cap construction is performed. +Exact grid sampling and zero-coordinate handling apply in both modes; uncapped +mesh arrays are not guaranteed to match legacy output. Enabled interior geometry +can differ because its crossing and QEF-solving conventions are stricter. + +## Best-Effort Contract + +Inspect `mesh.capping_report["closure_success"]` before treating a result as a +closed solid. Unsupported geometry emits `RuntimeWarning` and returns the +available mesh; it is not silently declared successful. + +- Fault stacks, fault-affected stacks, and stacks with removed extraction-mask + cells currently skip capping. Scalar samples alone do not encode enough cut + ownership to safely prevent bridging intentional fault or relation boundaries. + These meshes retain their extracted geometry and report `skipped_reason`. +- Missing refinement support is never repaired by extending a cap into the + interior. Missing boundary QEFs are reported and transitions are omitted. +- Escaped QEF vertices are reported and excluded from transition construction; + this implementation does not introduce a bounded QEF solver. +- Exact ties and multiple contour components in one cell remain best effort. + Edge-incidence, directed-winding, and vertex-link audits detect many failures, + but there is no geometric self-intersection test. +- New duplicate or zero-area triangles are removed. Invalid original geometry + is preserved and reported unsuccessful rather than silently repaired. +- `GEMPY_SKIP_TRIANGULATION` also skips capping and boundary evaluation. + +The full issue's general fault/mask-aware capping and all ambiguous topology +cases are not implemented by this first version. + +## Reports and Metadata + +Each mesh exposes `stack_index`, `surface_index`, `exported_surface_index`, +`isovalue`, and, when enabled, `inside_convention`. + +The capping report includes before/after open-edge counts, inferred non-extent +openings, added cap and transition triangle counts, removed triangle counts, +missing/escaped QEF counts, manifoldness and winding checks, maximum cap-plane +error, evaluated boundary-point count, and the original vertex count. + +`watertight` describes topology; `closure_success` additionally checks finite, +nondegenerate geometry and unresolved construction diagnostics. Neither proves +absence of self-intersections. Physical-opening classification for original +QEF edges uses owning-cell membership, not the interior QEF positions; it is +therefore a cell-level inference rather than proof of the opening's cause. +Cap-plane checks use `64 * float64_epsilon * max(1, abs(extent))` tolerance. + +`vertices_tensor` remains the pre-overlap differentiable QEF snapshot. Final +triangle indices reference `mesh.vertices`, **not** `vertices_tensor`. Appended +cap geometry and topology are NumPy postprocessing and are not differentiable. + +## Cost + +One unique finest-level six-face lattice is shared across surfaces. Boundary +scalar values are evaluated once per stack per batch and reused across that +stack's isosurfaces. The API module `API/dual_contouring/extent_capping.py` +orchestrates capping and restores the grid and stack cursor after boundary +evaluation, including failures. `evaluation_chunk_size` bounds the boundary point batch; +the existing evaluator also applies its own kernel-workload chunking. + +Boundary storage scales as `O(nx*ny + nx*nz + ny*nz)`, but Python dictionaries, +triangle connectivity, and topology audits have substantial additional memory +cost. This is a correctness-first implementation: there is no adaptive cap +simplification, cross-call cache, or reuse of existing corner samples. +Large-depth GPU and production-scale memory benchmarks remain outstanding. diff --git a/gempy_engine/API/dual_contouring/extent_capping.py b/gempy_engine/API/dual_contouring/extent_capping.py new file mode 100644 index 00000000..de205e65 --- /dev/null +++ b/gempy_engine/API/dual_contouring/extent_capping.py @@ -0,0 +1,76 @@ +"""Orchestrate extent capping and batched boundary scalar evaluation.""" + +import copy +import os + +import numpy as np + +from ..interp_single.interp_features import interpolate_all_fields_no_octree +from ...core.backend_tensor import BackendTensor +from ...core.data.engine_grid import EngineGrid +from ...core.data.generic_grid import GenericGrid +from ...core.data.stack_relation_type import StackRelationType +from ...modules.dual_contouring._extent_capping import boundary_lattice, cap_mesh + + +def cap_meshes_at_extent(all_meshes, dc_data_per_surface, all_mask_arrays, base_number, + orthogonal_extent, interpolation_input, options, data_descriptor): + """Cap extracted meshes in place, sharing boundary samples across surfaces.""" + skip_triangles = os.getenv("GEMPY_SKIP_TRIANGULATION", "0").lower() in ("true", "1", "t", "y", "yes") + if not all_meshes or skip_triangles: + return + + extent = BackendTensor.t.to_numpy(orthogonal_extent) + skip_reasons = [] + stack_relations = data_descriptor.stack_structure.masking_descriptor + fault_relations = data_descriptor.stack_structure.faults_relations + for mesh in all_meshes: + stack = mesh.stack_index + mask = all_mask_arrays[stack] + faulted = (stack_relations[stack] is StackRelationType.FAULT + or (fault_relations is not None and np.any(fault_relations[:, stack]))) + masked = mask is not None and not bool(mask.all()) + skip_reasons.append( + "Capping skipped: fault or extraction-mask boundary ownership is not supported" + if faulted or masked else None + ) + geometry = None + coordinates = points = np.empty((0, 3)) + boundary_scalars = None + if any(reason is None for reason in skip_reasons): + geometry = boundary_lattice(base_number, extent) + coordinates, points, _ = geometry + boundary_scalars = _interp_on_boundary(points, interpolation_input, options, data_descriptor) + for mesh, dc_data, reason in zip(all_meshes, dc_data_per_surface, skip_reasons): + cell_coordinates = BackendTensor.t.to_numpy(dc_data.left_right_codes[dc_data.valid_voxels]) + cap_mesh(mesh, cell_coordinates, base_number, extent, coordinates, + boundary_scalars[mesh.stack_index] if reason is None else np.empty(0), + mesh.isovalue, boundary_geometry=geometry, skip_reason=reason) + mesh.capping_report["boundary_scalar_points"] = len(points) if reason is None else 0 + mesh.capping_report["original_vertex_count"] = len(cell_coordinates) + + +def _interp_on_boundary(points, interpolation_input, options, data_descriptor): + """Evaluate each stack once per boundary batch, reusing it for its surfaces.""" + boundary_options = copy.deepcopy(options) + boundary_options.evaluation_options.compute_scalar = True + boundary_options.evaluation_options.compute_scalar_gradient = False + saved_grid = interpolation_input.grid + saved_stack = data_descriptor.stack_structure.stack_number + scalars = [np.empty(len(points)) for _ in range(data_descriptor.stack_structure.n_stacks)] + batch_size = max(1, int(options.evaluation_options.evaluation_chunk_size)) + try: + for start in range(0, len(points), batch_size): + stop = min(start + batch_size, len(points)) + interpolation_input.set_temp_grid(EngineGrid(custom_grid=GenericGrid( + values=BackendTensor.t.array(points[start:stop], dtype=BackendTensor.dtype) + ))) + outputs = interpolate_all_fields_no_octree(interpolation_input, boundary_options, data_descriptor) + for stack_index, output in enumerate(outputs): + scalars[stack_index][start:stop] = BackendTensor.t.to_numpy( + output.exported_fields.scalar_field[output.grid.custom_grid_slice] + ) + finally: + interpolation_input.set_temp_grid(saved_grid) + data_descriptor.stack_structure.stack_number = saved_stack + return scalars 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 3a357871..c5f0b12c 100644 --- a/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py +++ b/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py @@ -24,7 +24,8 @@ get_masked_codes, mask_generation) from ...modules.dual_contouring.overlapping import average_overlapping_vertices, remove_fault_overlap_triangles from ...modules.dual_contouring._support_report import mesh_support_report -from ...core.data.options.evaluation_options import OctreeRefinementMode +from ...core.data.options.evaluation_options import OctreeRefinementMode, MeshExtentCapping +from .extent_capping import cap_meshes_at_extent @gempy_profiler_decorator @@ -52,6 +53,7 @@ def dual_contouring_multi_scalar( octree_leaves = octree_list[-1] all_meshes: List[DualContouringMesh] = [] + cap_enabled = MeshExtentCapping(options.evaluation_options.mesh_extraction_extent_capping) != MeshExtentCapping.NONE dual_contouring_options = copy.deepcopy(options) dual_contouring_options.evaluation_options.compute_scalar_gradient = True @@ -86,7 +88,8 @@ def dual_contouring_multi_scalar( _xyz_corners=octree_leaves.grid.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 + masking=mask, + strict_crossings=cap_enabled ) all_surfaces_intersection.append(intersection_xyz) @@ -110,6 +113,7 @@ def dual_contouring_multi_scalar( # Generate meshes for each scalar field dc_data_per_surface_all = [] support_reports = [] + surface_metadata = [] stack_relations = data_descriptor.stack_structure.masking_descriptor for n_scalar_field in range(data_descriptor.stack_structure.n_stacks): if stack_relations[n_scalar_field] is StackRelationType.NULL_SPACE: @@ -125,7 +129,8 @@ def dual_contouring_multi_scalar( output.exported_fields.scalar_field[output.grid.corners_grid_slice], output.scalar_field_at_sp[surface_i], base_number, mask, surface_index=surface_i, - ancestor_coordinates=[level.grid.octree_grid.integer_coordinates for level in octree_list[:-1]] + ancestor_coordinates=[level.grid.octree_grid.integer_coordinates for level in octree_list[:-1]], + strict_crossings=cap_enabled ) report['stack_index'] = n_scalar_field if report['internal_refinement_boundary_edge_count']: @@ -154,11 +159,13 @@ def dual_contouring_multi_scalar( tree_depth=options.number_octree_levels, base_number=base_number, triangulation_method=options.evaluation_options.triangulation_method, - generated_cell_coordinates=left_right_codes + generated_cell_coordinates=left_right_codes, + strict_crossings=cap_enabled ) dc_data_per_surface_all.append(dc_data_per_surface) surface_to_stack.append(n_scalar_field) + surface_metadata.append((n_scalar_field, surface_i, float(output.scalar_field_at_sp[surface_i]))) if compute_overlap: left_right_per_mesh.append(all_left_right_codes[n_scalar_field][dc_data_per_surface.valid_voxels]) @@ -213,6 +220,20 @@ def dual_contouring_multi_scalar( mesh.vertices = BackendTensor.t.to_numpy(mesh.vertices) mesh.edges = BackendTensor.t.to_numpy(mesh.edges) + for index, (mesh, (stack_index, surface_index, isovalue)) in enumerate(zip(all_meshes, surface_metadata)): + mesh.stack_index = stack_index + mesh.surface_index = surface_index + mesh.exported_surface_index = index + mesh.isovalue = isovalue + mesh.inside_convention = "scalar <= isovalue" if cap_enabled else None + + if cap_enabled: + cap_meshes_at_extent( + all_meshes, dc_data_per_surface_all, all_mask_arrays, base_number, + octree_list[0].grid.octree_grid.orthogonal_extent, + interpolation_input, options, data_descriptor + ) + return all_meshes diff --git a/gempy_engine/core/data/__init__.py b/gempy_engine/core/data/__init__.py index 2715672e..acc87081 100644 --- a/gempy_engine/core/data/__init__.py +++ b/gempy_engine/core/data/__init__.py @@ -3,5 +3,5 @@ from .kernel_classes.orientations import Orientations, OrientationsInternals from .kernel_classes.surface_points import SurfacePoints, SurfacePointsInternals from .options.interpolation_options import InterpolationOptions -from .options.evaluation_options import OctreeRefinementMode +from .options.evaluation_options import OctreeRefinementMode, MeshExtentCapping from .solutions import Solutions diff --git a/gempy_engine/core/data/dual_contouring_data.py b/gempy_engine/core/data/dual_contouring_data.py index b43f22cd..b794e422 100644 --- a/gempy_engine/core/data/dual_contouring_data.py +++ b/gempy_engine/core/data/dual_contouring_data.py @@ -32,6 +32,7 @@ class DualContouringData: triangulation_method: TriangulationMethod = TriangulationMethod.LEGACY generated_cell_coordinates: Optional[np.ndarray] = None # Before geological masking. triangulation_report: dict = field(default_factory=dict) # Quad support before overlap/fault triangle removal. + strict_crossings: bool = False @property def valid_voxels(self): diff --git a/gempy_engine/core/data/dual_contouring_mesh.py b/gempy_engine/core/data/dual_contouring_mesh.py index 692ad997..25a6108e 100644 --- a/gempy_engine/core/data/dual_contouring_mesh.py +++ b/gempy_engine/core/data/dual_contouring_mesh.py @@ -11,6 +11,12 @@ class DualContouringMesh: edges: np.ndarray dc_data: Optional[DualContouringData] = None # * In principle we need this just for testing support_report: Optional[dict] = None + capping_report: Optional[dict] = None + stack_index: Optional[int] = None + surface_index: Optional[int] = None + exported_surface_index: Optional[int] = None + isovalue: Optional[float] = None + inside_convention: Optional[str] = None def __repr__(self): return f"DualContouringMesh({self.vertices.shape[0]} vertices, {self.edges.shape[0]} edges)" diff --git a/gempy_engine/core/data/options/evaluation_options.py b/gempy_engine/core/data/options/evaluation_options.py index c9aad945..6ed35f60 100644 --- a/gempy_engine/core/data/options/evaluation_options.py +++ b/gempy_engine/core/data/options/evaluation_options.py @@ -15,6 +15,11 @@ class OctreeRefinementMode(str, enum.Enum): PRECISE = "precise" +class MeshExtentCapping(str, enum.Enum): + NONE = "none" + SCALAR_LESS_EQUAL = "scalar_less_equal" + + class MeshExtractionMaskingOptions(enum.Enum): NOTHING = enum.auto() # * This is only for testing DISJOINT = enum.auto() @@ -39,6 +44,7 @@ class EvaluationOptions: triangulation_method: TriangulationMethod = TriangulationMethod.LEGACY mesh_extraction: bool = True + mesh_extraction_extent_capping: MeshExtentCapping = MeshExtentCapping.NONE mesh_extraction_masking_options: MeshExtractionMaskingOptions = MeshExtractionMaskingOptions.INTERSECT mesh_extraction_fancy: Annotated[bool, deprecated("Old extraction method not in use anymore")] = True diff --git a/gempy_engine/core/data/raw_arrays_solution.py b/gempy_engine/core/data/raw_arrays_solution.py index b7f7706c..ef7f31ce 100644 --- a/gempy_engine/core/data/raw_arrays_solution.py +++ b/gempy_engine/core/data/raw_arrays_solution.py @@ -123,7 +123,7 @@ def meshes_to_subsurface(self, input_transform: Transform | None = None, element pd = require_pandas() vertex: list[np.ndarray] = self.vertices - simplex_list: list[np.ndarray] = self.edges + simplex_list: list[np.ndarray] = [triangles.copy() for triangles in self.edges] idx_max = 0 for i, simplex_array in enumerate(simplex_list): @@ -287,4 +287,3 @@ def _fill_finite_fault_scalar_fields_with_dense_grid( ] raw_arrays_solution.finite_fault_scalar_field_matrix = np.vstack(fields) raw_arrays_solution.finite_fault_stack_indices = stack_indices - diff --git a/gempy_engine/core/data/regular_grid.py b/gempy_engine/core/data/regular_grid.py index 80cabd71..eb510785 100644 --- a/gempy_engine/core/data/regular_grid.py +++ b/gempy_engine/core/data/regular_grid.py @@ -28,7 +28,7 @@ def __len__(self): def __post_init__(self): self.regular_grid_shape = BackendTensor.t.array(self.regular_grid_shape) - self.orthogonal_extent = BackendTensor.t.array(self.orthogonal_extent) + 1e-6 # * This to avoid some errors evaluating in 0 (e.g. bias in dual contouring) + self.orthogonal_extent = BackendTensor.t.array(self.orthogonal_extent) self._create_regular_grid_3d() @@ -48,15 +48,15 @@ def dz(self): @property def x_coord(self): - return BackendTensor.t.linspace(self.orthogonal_extent[0] + self.dx / 2, self.orthogonal_extent[1] - self.dx / 2, self.resolution[0]) + return BackendTensor.t.linspace(self.orthogonal_extent[0] + self.dx / 2, self.orthogonal_extent[1] - self.dx / 2, self.resolution[0], dtype=BackendTensor.dtype_obj) @property def y_coord(self): - return BackendTensor.t.linspace(self.orthogonal_extent[2] + self.dy / 2, self.orthogonal_extent[3] - self.dy / 2, self.resolution[1]) + return BackendTensor.t.linspace(self.orthogonal_extent[2] + self.dy / 2, self.orthogonal_extent[3] - self.dy / 2, self.resolution[1], dtype=BackendTensor.dtype_obj) @property def z_coord(self): - return BackendTensor.t.linspace(self.orthogonal_extent[4] + self.dz / 2, self.orthogonal_extent[5] - self.dz / 2, self.resolution[2]) + return BackendTensor.t.linspace(self.orthogonal_extent[4] + self.dz / 2, self.orthogonal_extent[5] - self.dz / 2, self.resolution[2], dtype=BackendTensor.dtype_obj) @classmethod def from_octree_level(cls, xyz_coords_octree: np.ndarray, previous_regular_grid: "RegularGrid", diff --git a/gempy_engine/modules/dual_contouring/_extent_capping.py b/gempy_engine/modules/dual_contouring/_extent_capping.py new file mode 100644 index 00000000..96e94211 --- /dev/null +++ b/gempy_engine/modules/dual_contouring/_extent_capping.py @@ -0,0 +1,368 @@ +"""Best-effort NumPy caps for a uniform dual-contouring lattice. + +This is not a mesh repair operation: only the six supplied extent faces are +sampled. Internal mask/fault holes remain untouched. Ambiguous/tied contours, +missing cells and unconstrained QEFs can leave holes; inspect ``capping_report``. +""" + +from collections import defaultdict +from itertools import product +import warnings + +import numpy as np + +from ._scalar_crossing import scalar_crossing_parameters + + +def boundary_lattice(shape, extent): + """Return integer corners, physical positions, and face metadata. + + ``shape`` is the three *cell* counts. ``extent`` is + ``(xmin, xmax, ymin, ymax, zmin, zmax)``. Corners are unique globally, + including box edges/corners. Metadata contains ``triangles`` (indices into + these arrays) and ``cells`` (one owning in-domain cell per triangle). + Faces are triangulated with a consistent low-to-high diagonal and outward + winding. Positions use the exact supplied bounds, not QEF-derived bounds. + + Storage scales with the full boundary area, not the active contour. The + Python lists/dictionaries used during construction have substantial peak + memory overhead on large grids. Reuse this tuple via ``boundary_geometry`` + rather than constructing it twice; this is not a streaming implementation. + """ + shape = np.asarray(shape) + bounds = np.asarray(extent, dtype=float).reshape(3, 2) + if shape.shape != (3,) or np.any(shape < 1) or np.any(shape != shape.astype(int)): + raise ValueError("shape must contain three positive integer cell counts") + if not np.isfinite(bounds).all() or np.any(bounds[:, 1] <= bounds[:, 0]): + raise ValueError("extent must contain finite increasing bounds") + shape = shape.astype(int) + coordinates, triangles, cells, lookup = [], [], [], {} + for axis in range(3): + u, v = (axis + 1) % 3, (axis + 2) % 3 + for side in (0, 1): + for i, j in product(range(shape[u]), range(shape[v])): + corners = [] + for du, dv in ((0, 0), (1, 0), (1, 1), (0, 1)): + c = [0, 0, 0] + c[axis], c[u], c[v] = side * shape[axis], i + du, j + dv + key = tuple(c) + if key not in lookup: + lookup[key] = len(coordinates) + coordinates.append(c) + corners.append(lookup[key]) + cell = [0, 0, 0] + cell[axis], cell[u], cell[v] = side * (shape[axis] - 1), i, j + for a, b, c in ((0, 1, 2), (0, 2, 3)): + triangles.append([corners[a], corners[b if side else c], corners[c if side else b]]) + cells.append(cell.copy()) + coordinates = np.asarray(coordinates, dtype=np.int64) + positions = bounds[:, 0] + coordinates / shape * (bounds[:, 1] - bounds[:, 0]) + for axis in range(3): + positions[coordinates[:, axis] == shape[axis], axis] = bounds[axis, 1] + return coordinates, positions, { + "triangles": np.asarray(triangles, dtype=np.int64), + "cells": np.asarray(cells, dtype=np.int64), + } + + +def cap_mesh(mesh, cell_coordinates, shape, extent, boundary_coordinates, + boundary_scalar, isovalue, *, boundary_geometry=None, skip_reason=None): + """Mutate NumPy ``mesh.vertices``/``mesh.edges`` and return ``mesh``. + + ``edges`` means triangle indices (M, 3). ``cell_coordinates`` is (N, 3) + in the exact row order of the original QEF vertices, not sorted order. + Boundary samples may be permuted, but must cover ``boundary_lattice``. + The retained solid is ``scalar <= isovalue``. Original vertices/triangles + keep their indices and winding. ``mesh.capping_report`` records counts and + whole-mesh topology audits; a RuntimeWarning summarizes limitations. + Invalid API array shapes raise ValueError; unsupported geometry does not. + ``boundary_geometry`` may be the unmodified tuple from ``boundary_lattice`` + for this same shape/extent, avoiding a second lattice construction. + + Report ``open_edges_before/after`` counts incidence-one edges. + ``non_extent_open_edges_before/after`` excludes edges whose endpoints share + an extent face: original QEF membership comes from owning cells, whereas + new cap points use physical positions. This is a cell-level classification, + not proof that an opening in a boundary cell was caused by the extent. + ``watertight`` audits topology only; ``closure_success`` additionally + requires valid finite geometry and no missing/escaped/unsupported cases. + ``added_cap_triangles`` and ``added_transition_triangles`` count retained + triangles after cleanup. ``cap_max_plane_error`` is the maximum, over cap + triangles, of the distance to their nearest common extent plane (physical + units, zero if no cap). Invalid original geometry is preserved and reported + unsuccessful without attempting additions. + """ + report = dict(added_vertices=0, added_triangles=0, contour_edges=0, + missing_qef=0, escaped_qef=0, unsupported_topology=0, + discarded_components=0, removed_triangles=0, warnings=[], + added_cap_triangles=0, added_transition_triangles=0, + cap_max_plane_error=0.0, closure_success=False) + mesh.capping_report = report + + def audit(): + result = {} + v = np.asarray(mesh.vertices) + triangles = np.asarray(mesh.edges) + valid = np.all(np.isfinite(triangles) & (triangles >= 0) & (triangles < len(v)) + & (triangles == np.floor(triangles)), axis=1) + result["invalid_triangles"] = int(np.count_nonzero(~valid)) + result["nonfinite_vertices"] = int(np.count_nonzero(~np.isfinite(v).all(axis=1))) + triangles = triangles[valid].astype(np.int64) + p = v[triangles] + with np.errstate(over="ignore", invalid="ignore"): + cross = np.cross(p[:, 1] - p[:, 0], p[:, 2] - p[:, 0]) + result["nonfinite_triangles"] = int(np.count_nonzero(~np.isfinite(cross).all(axis=1))) + result["degenerate_triangles"] = int(np.count_nonzero(np.all(cross == 0, axis=1))) + membership = np.zeros((len(v), 3, 2), dtype=bool) + membership[:len(cells), :, 0] = cells == 0 + membership[:len(cells), :, 1] = cells == shape - 1 + tolerance = 64 * np.finfo(float).eps * np.maximum(1, np.abs(bounds)) + membership[len(cells):] = np.abs(v[len(cells):, :, None] - bounds) <= tolerance + incidence = defaultdict(list) + links = defaultdict(list) + for a, b, c in triangles: + for x, y in ((a, b), (b, c), (c, a)): + incidence[tuple(sorted((int(x), int(y))))].append((int(x), int(y))) + for x, y, z in ((a, b, c), (b, c, a), (c, a, b)): + links[int(x)].append((int(y), int(z))) + openings = [edge for edge, entries in incidence.items() if len(entries) == 1] + result["open_edges"] = len(openings) + result["non_extent_open_edges"] = sum(not np.any(membership[a] & membership[b]) for a, b in openings) + result["nonmanifold_edges"] = sum(len(e) > 2 for e in incidence.values()) + result["orientation_conflicts"] = sum(len(e) == 2 and e[0] == e[1] for e in incidence.values()) + bad_links = 0 + for edges in links.values(): + graph = defaultdict(list) + for a, b in edges: + graph[a].append(b) + graph[b].append(a) + seen, todo = set(), [next(iter(graph))] + while todo: + node = todo.pop() + if node not in seen: + seen.add(node) + todo.extend(graph[node]) + degrees = [len(neighbors) for neighbors in graph.values()] + if len(seen) != len(graph) or max(degrees) > 2 or degrees.count(1) not in (0, 2): + bad_links += 1 + result["nonmanifold_vertices"] = bad_links + result["watertight"] = bool(incidence) and not any(result[k] for k in ( + "invalid_triangles", "open_edges", "nonmanifold_edges", "orientation_conflicts", "nonmanifold_vertices")) + return result + + def finish(): + after = audit() + report.update(after) + report["open_edges_after"] = after["open_edges"] + report["non_extent_open_edges_after"] = after["non_extent_open_edges"] + report["boundary_edges"] = after["open_edges"] + report["closure_success"] = report["watertight"] and not any(report[k] for k in ( + "missing_qef", "escaped_qef", "unsupported_topology", "invalid_triangles", + "nonfinite_vertices", "nonfinite_triangles", "degenerate_triangles")) + plane_tolerance = 64 * np.finfo(float).eps * max(1, np.max(np.abs(bounds))) + if report["cap_max_plane_error"] > plane_tolerance: + report["closure_success"] = False + report["warnings"].append("Cap triangles are off the supplied extent planes") + for key in ("missing_qef", "escaped_qef", "unsupported_topology", "boundary_edges", + "nonmanifold_edges", "orientation_conflicts", "nonmanifold_vertices", + "invalid_triangles", "nonfinite_vertices", "nonfinite_triangles", "degenerate_triangles"): + if report[key]: + report["warnings"].append(f"{key}: {report[key]}") + if report["warnings"]: + warnings.warn("Extent capping best effort: " + "; ".join(report["warnings"]), + RuntimeWarning, stacklevel=2) + return mesh + + if not isinstance(mesh.vertices, np.ndarray) or not isinstance(mesh.edges, np.ndarray): + report["unsupported_topology"] += 1 + report["warnings"].append("Only NumPy meshes are supported") + warnings.warn(report["warnings"][0], RuntimeWarning, stacklevel=2) + return mesh + vertices = np.asarray(mesh.vertices) + cells = np.asarray(cell_coordinates) + if vertices.ndim != 2 or vertices.shape[1] != 3 or cells.shape != vertices.shape: + raise ValueError("cell_coordinates must have one (x, y, z) row per QEF vertex") + if mesh.edges.ndim != 2 or mesh.edges.shape[1] != 3: + raise ValueError("mesh.edges must be triangle indices (M, 3)") + bounds = np.asarray(extent, dtype=float).reshape(3, 2) + shape = np.asarray(shape) + if shape.shape != (3,) or np.any(shape < 1) or not np.isfinite(shape).all() or np.any(shape != np.floor(shape)): + raise ValueError("shape must contain three positive integer cell counts") + if not np.isfinite(bounds).all() or np.any(bounds[:, 1] <= bounds[:, 0]): + raise ValueError("extent must contain finite increasing bounds") + shape = shape.astype(int) + before = audit() + report.update({key + "_before": value for key, value in before.items()}) + if skip_reason is not None: + report["unsupported_topology"] += 1 + report["warnings"].append(skip_reason) + report["skipped_reason"] = skip_reason + return finish() + if any(before[k] for k in ("invalid_triangles", "nonfinite_vertices", "nonfinite_triangles", "degenerate_triangles")): + return finish() + coords, positions, faces = (boundary_lattice(shape, extent) if boundary_geometry is None + else boundary_geometry) + supplied = np.asarray(boundary_coordinates) + scalar = np.asarray(boundary_scalar, dtype=float).reshape(-1) + if supplied.shape != (len(scalar), 3): + raise ValueError("boundary_coordinates and boundary_scalar must have matching rows") + samples = {tuple(c): s for c, s in zip(supplied, scalar)} + if len(samples) != len(supplied) or len(samples) != len(coords) or any(tuple(c) not in samples for c in coords): + report["unsupported_topology"] += 1 + report["warnings"].append("Incomplete or duplicate boundary lattice") + return finish() + values = np.asarray([samples[tuple(c)] for c in coords]) + # Canonical endpoint order gives every shared edge the identical parameter. + source = faces["triangles"] + ends = np.roll(source, -1, axis=1) + starts, ends = np.minimum(source, ends), np.maximum(source, ends) + try: + with np.errstate(over="ignore", invalid="ignore"): + crossing, parameters = scalar_crossing_parameters(values[starts], values[ends], isovalue) + except ValueError as error: + report["unsupported_topology"] += 1 + report["warnings"].append(f"Nonfinite boundary crossing: {error}") + return finish() + step = (bounds[:, 1] - bounds[:, 0]) / shape + qef = {} + seen_cells = set() + for i, cell in enumerate(cells): + key = tuple(cell) + if np.any(cell != np.floor(cell)) or np.any(cell < 0) or np.any(cell >= shape) or key in seen_cells: + report["unsupported_topology"] += 1 + continue + seen_cells.add(key) + low = bounds[:, 0] + cell * step + tolerance = 64 * np.finfo(float).eps * np.maximum(1, np.maximum(abs(low), abs(low + step))) + if not np.isfinite(vertices[i]).all() or np.any(vertices[i] < low - tolerance) or np.any(vertices[i] > low + step + tolerance): + report["escaped_qef"] += 1 + continue + qef[key] = i + + points = list(vertices.astype(float)) + point_ids = {} + + def point(a, b=None, t=None): + # Endpoint ties must share the lattice-point key, not an edge key. + if b is None or values[a] == isovalue: + key, p = ("p", a), positions[a] + elif values[b] == isovalue: + key, p = ("p", b), positions[b] + else: + a, b = sorted((a, b)) + if t == 0 or t == 1: + endpoint = a if t == 0 else b + key, p = ("p", endpoint), positions[endpoint] + else: + key = ("e", a, b) # Shared by primal edges and face diagonals. + p = positions[a] + t * (positions[b] - positions[a]) + if key not in point_ids: + point_ids[key] = len(points) + points.append(p) + return point_ids[key] + + cap, owners = [], [] + for i, (triangle, cell) in enumerate(zip(faces["triangles"], faces["cells"])): + polygon = [] + for j, (a, b) in enumerate(zip(triangle, np.roll(triangle, -1))): + if values[a] <= isovalue: + polygon.append(point(a)) + if crossing[i, j]: + polygon.append(point(a, b, parameters[i, j])) + polygon = list(dict.fromkeys(polygon)) + for j in range(1, len(polygon) - 1): + cap.append((polygon[0], polygon[j], polygon[j + 1])) + owners.append(tuple(cell)) + + # Connected cap patches are seeded only by contour edges with a usable QEF. + incidence = defaultdict(list) + for i, triangle in enumerate(cap): + for a, b in zip(triangle, np.roll(triangle, -1)): + incidence[tuple(sorted((a, b)))].append((i, a, b)) + adjacency = defaultdict(set) + contours, seeds = [], set() + for entries in incidence.values(): + if len(entries) == 2: + i, j = entries[0][0], entries[1][0] + adjacency[i].add(j) + adjacency[j].add(i) + elif len(entries) == 1: + i, a, b = entries[0] + report["contour_edges"] += 1 + q = qef.get(owners[i]) + if q is None: + report["missing_qef"] += 1 + else: + seeds.add(i) + contours.append((i, a, b, q)) + else: + report["unsupported_topology"] += 1 + retained, visited = set(), set() + for start in range(len(cap)): + if start in visited: + continue + component, todo = set(), [start] + while todo: + i = todo.pop() + if i not in component: + component.add(i) + todo.extend(adjacency[i] - component) + visited.update(component) + if component & seeds: + retained.update(component) + else: + report["discarded_components"] += 1 + added = [cap[i] for i in sorted(retained)] + cap_count = len(added) + fans = [(b, a, q) for i, a, b, q in contours if i in retained] + added.extend(fans) + + # The two in-domain cells around a boundary primal edge need one more + # triangle. Derive its winding from the existing directed face-fan edges. + spokes = defaultdict(list) + for b, a, q in fans: + spokes[(a, q)].append((a, q)) + spokes[(b, q)].append((q, b)) + residual = defaultdict(list) + for (p, q), edges in spokes.items(): + forward = sum(a == p for a, b in edges) + reverse = len(edges) - forward + if forward != reverse: + residual[p].extend([(p, q) if forward > reverse else (q, p)] * abs(forward - reverse)) + for p, edges in residual.items(): + if len(edges) == 2: + incoming = [a for a, b in edges if b == p] + outgoing = [b for a, b in edges if a == p] + if (len(incoming) == len(outgoing) == 1 + and np.sum(np.abs(cells[incoming[0]] - cells[outgoing[0]])) == 1): + added.append((p, incoming[0], outgoing[0])) + continue + report["unsupported_topology"] += 1 + + # Never remove/reindex original geometry, even if it already has defects. + seen = {tuple(sorted(t)) for t in mesh.edges} + clean = [] + for index, triangle in enumerate(added): + key = tuple(sorted(triangle)) + a, b, c = (points[i] for i in triangle) + cross = np.cross(b - a, c - a) + if len(set(triangle)) < 3 or key in seen or not np.any(cross) or not np.isfinite(cross).all(): + report["removed_triangles"] += 1 + continue + seen.add(key) + clean.append(triangle) + if index < cap_count: + report["added_cap_triangles"] += 1 + # A cap triangle must have all three vertices on the same plane. + distance = np.abs(np.asarray([a, b, c])[:, :, None] - bounds) + report["cap_max_plane_error"] = max(report["cap_max_plane_error"], float(distance.max(axis=0).min())) + else: + report["added_transition_triangles"] += 1 + used = sorted({i for t in clean for i in t if i >= len(vertices)}) + remap = {old: len(vertices) + i for i, old in enumerate(used)} + if clean: + mesh.vertices = np.concatenate((vertices, np.asarray([points[i] for i in used]).reshape(-1, 3))) + new_edges = np.asarray([[remap.get(i, i) for i in t] for t in clean], dtype=np.int64) + mesh.edges = np.concatenate((mesh.edges, new_edges)) + report["added_vertices"], report["added_triangles"] = len(used), len(clean) + return finish() diff --git a/gempy_engine/modules/dual_contouring/_gen_vertices.py b/gempy_engine/modules/dual_contouring/_gen_vertices.py index a622e25b..a0378975 100644 --- a/gempy_engine/modules/dual_contouring/_gen_vertices.py +++ b/gempy_engine/modules/dual_contouring/_gen_vertices.py @@ -29,21 +29,12 @@ def generate_dual_contouring_vertices(dc_data_per_stack: DualContouringData, sli edges_xyz[:, :12][valid_edges_bool] = xyz_on_edge edges_normals[:, :12][valid_edges_bool] = gradients - # Use nanmean directly without intermediate copy + # Zero coordinates are valid samples, not missing edge constraints. bias_xyz_slice = edges_xyz[:, :12] - - if BackendTensor.engine_backend == AvailableBackends.PYTORCH: - mask = bias_xyz_slice == 0 - bias_xyz_masked = BackendTensor.tfnp.where(mask, float('nan'), bias_xyz_slice) - mass_points = BackendTensor.tfnp.nanmean(bias_xyz_masked, axis=1) - else: - # NumPy: more efficient approach using sum and count - mask = bias_xyz_slice != 0 - sum_valid = (bias_xyz_slice * mask).sum(axis=1) - count_valid = mask.sum(axis=1) - # Avoid division by zero - count_valid = BackendTensor.tfnp.maximum(count_valid, 1) - mass_points = sum_valid / count_valid + mask = valid_edges_bool[:, :, None] + sum_valid = (bias_xyz_slice * mask).sum(axis=1) + count_valid = mask.sum(axis=1) + mass_points = sum_valid / count_valid # Assign mass points to bias positions edges_xyz[:, 12:15] = mass_points[:, None, :] @@ -96,7 +87,10 @@ def generate_dual_contouring_vertices(dc_data_per_stack: DualContouringData, sli # Solve ATA @ x = ATb (use solve instead of inv for numerical stability) import torch - reg = 1e-4 * torch.eye(3, device=ATA.device, dtype=ATA.dtype).unsqueeze(0) + # The mass-point constraints already make ATA positive definite. Extra + # origin-centered regularization biases the capped solid's volume. + strength = 0 if dc_data_per_stack.strict_crossings else 1e-4 + reg = strength * torch.eye(3, device=ATA.device, dtype=ATA.dtype).unsqueeze(0) vertices = torch.linalg.solve(ATA + reg, ATb) else: # NumPy: use efficient einsum diff --git a/gempy_engine/modules/dual_contouring/_scalar_crossing.py b/gempy_engine/modules/dual_contouring/_scalar_crossing.py new file mode 100644 index 00000000..e17c1a71 --- /dev/null +++ b/gempy_engine/modules/dual_contouring/_scalar_crossing.py @@ -0,0 +1,32 @@ +"""Shared scalar-only edge classification for contouring and NumPy cappers.""" + +import numpy as np + + +def scalar_crossing_parameters(start, end, iso, *, xp=np): + """Return broadcast ``(crossing, t)`` for ``(1-t)*start + t*end == iso``. + + Inputs are floating arrays in the supplied NumPy or PyTorch namespace. + ``scalar <= iso`` is inside: an entirely iso-valued edge never crosses, + and an equal endpoint crosses only if the other endpoint is outside. + Parameters are clamped to [0, 1], snapped within 1e-12 of an endpoint, + and zero on non-crossing edges. + Nonfinite inputs or interpolation arithmetic raise ValueError. + The default namespace is NumPy, independent of the engine backend. + """ + for name, value in (("start", start), ("end", end), ("iso", iso)): + if not bool(xp.isfinite(value).all()): + raise ValueError(f"Scalar crossing requires finite {name} values") + + crossing = (start <= iso) != (end <= iso) + # Neutralize unused edges before dividing, without perturbing crossings. + denominator = xp.where(crossing, end, 0) - xp.where(crossing, start, 0) + numerator = xp.where(crossing, iso, 0) - xp.where(crossing, start, 0) + if not bool((xp.isfinite(denominator) & xp.isfinite(numerator)).all()): + raise ValueError("Nonfinite scalar crossing interpolation arithmetic") + t = numerator / xp.where(crossing, denominator, 1) + if not bool(xp.isfinite(t).all()): + raise ValueError("Nonfinite scalar crossing interpolation parameter") + t = xp.clip(t, 0, 1) + t = xp.where(t <= 1e-12, 0, xp.where(t >= 1 - 1e-12, 1, t)) + return crossing, t diff --git a/gempy_engine/modules/dual_contouring/_support_report.py b/gempy_engine/modules/dual_contouring/_support_report.py index 08a32bac..20f8aac7 100644 --- a/gempy_engine/modules/dual_contouring/_support_report.py +++ b/gempy_engine/modules/dual_contouring/_support_report.py @@ -7,7 +7,7 @@ def mesh_support_report(coordinates, scalar_corners, isovalue, domain_shape, - mask=None, surface_index=0, ancestor_coordinates=()): + mask=None, surface_index=0, ancestor_coordinates=(), strict_crossings=False): """Classify missing incident cells for unique sampled sign-changing edges. This diagnoses sampled crossings, not unsampled components or later triangle @@ -30,7 +30,10 @@ def mesh_support_report(coordinates, scalar_corners, isovalue, domain_shape, ((0, 1), (2, 3), (4, 5), (6, 7)), )): for a, b in pairs: - crossing = (scalar[:, a] >= iso) != (scalar[:, b] >= iso) + if strict_crossings: + crossing = (scalar[:, a] <= iso) != (scalar[:, b] <= iso) + else: + crossing = (scalar[:, a] >= iso) != (scalar[:, b] >= iso) edges.update((direction, *p) for p in coords[crossing] + corners[a]) report = dict(surface_index=surface_index, crossing_edge_count=len(edges), missing_incident_cell_count=0, physical_boundary_edge_count=0, diff --git a/gempy_engine/modules/dual_contouring/dual_contouring_interface.py b/gempy_engine/modules/dual_contouring/dual_contouring_interface.py index 6bcb9ce1..d1f75aef 100644 --- a/gempy_engine/modules/dual_contouring/dual_contouring_interface.py +++ b/gempy_engine/modules/dual_contouring/dual_contouring_interface.py @@ -5,6 +5,7 @@ from ._find_vertex_overlap import find_repeated_voxels_across_stacks from ._apply_vertex_overlap_logic import apply_relations_to_overlaps +from ._scalar_crossing import scalar_crossing_parameters from .fancy_triangulation import get_left_right_array from ...config import AvailableBackends from ...core.backend_tensor import BackendTensor @@ -20,13 +21,44 @@ # region edges def find_intersection_on_edge(_xyz_corners, scalar_field_on_corners, - scalar_at_sp, masking=None) -> Tuple: + scalar_at_sp, masking=None, *, strict_crossings=False) -> Tuple: + """Find edge intersections, optionally using exact scalar-side classification.""" + if strict_crossings: + return _find_strict_intersections(_xyz_corners, scalar_field_on_corners, scalar_at_sp, masking) if BackendTensor.engine_backend == AvailableBackends.PYTORCH: return find_intersection_on_edge_torch(_xyz_corners, scalar_field_on_corners, scalar_at_sp, masking) else: return find_intersection_on_edge_numpy(_xyz_corners, scalar_field_on_corners, scalar_at_sp, masking) +def _find_strict_intersections(xyz_corners, scalars, isovalues, masking): + xp = BackendTensor.t + xyz = xyz_corners.reshape(-1, 8, 3) + scalars = scalars.reshape(1, -1, 8) + if masking is not None: + xyz = xyz[masking] + scalars = scalars[:, masking] + if not bool(xp.isfinite(xyz).all()): + raise ValueError("Strict edge crossings require finite corner coordinates") + + # Preserve the legacy x/y/z edge ordering and start-to-end direction. + start = [4, 5, 6, 7, 2, 3, 6, 7, 1, 3, 5, 7] + end = [0, 1, 2, 3, 0, 1, 4, 5, 0, 2, 4, 6] + valid, t = scalar_crossing_parameters( + scalars[:, :, start], scalars[:, :, end], isovalues.reshape(-1, 1, 1), xp=xp + ) + shape = (*valid.shape, 3) + a = xp.broadcast_to(xyz[None, :, start, :], shape)[valid] + b = xp.broadcast_to(xyz[None, :, end, :], shape)[valid] + weights = t[valid][:, None] + intersections = (1 - weights) * a + weights * b + if not bool(xp.isfinite(intersections).all()): + raise ValueError("Nonfinite strict edge intersection coordinates") + if BackendTensor.engine_backend == AvailableBackends.PYTORCH: + return intersections, valid.reshape(-1) + return intersections, valid.reshape(-1, 12) + + def find_intersection_on_edge_numpy(_xyz_corners, scalar_field_on_corners, scalar_at_sp, masking=None) -> Tuple: """This function finds all the intersections for multiple layers per series diff --git a/tests/test_common/test_modules/test_extent_capping.py b/tests/test_common/test_modules/test_extent_capping.py new file mode 100644 index 00000000..df24c748 --- /dev/null +++ b/tests/test_common/test_modules/test_extent_capping.py @@ -0,0 +1,342 @@ +"""Analytic interior DC meshes, independent of the capping triangulator.""" + +from itertools import product +from types import SimpleNamespace +import warnings + +import numpy as np +import pytest + +from gempy_engine.modules.dual_contouring._extent_capping import boundary_lattice, cap_mesh +from gempy_engine.modules.dual_contouring import _extent_capping + + +def plane_mesh(shape, extent, normal, level): + bounds = np.asarray(extent).reshape(3, 2) + shape = np.asarray(shape) + normal = np.asarray(normal) + step = (bounds[:, 1] - bounds[:, 0]) / shape + positions = lambda c: bounds[:, 0] + np.asarray(c) * step + scalar = lambda c: positions(c) @ normal - level + vertices, cells, lookup = [], [], {} + for cell in product(*(range(n) for n in shape)): + intersections = [] + for axis in range(3): + others = [i for i in range(3) if i != axis] + for offsets in product((0, 1), repeat=2): + a = np.array(cell) + a[others] += offsets + b = a.copy() + b[axis] += 1 + sa, sb = scalar(a), scalar(b) + if (sa <= 0) != (sb <= 0): + intersections.append(positions(a) + sa / (sa - sb) * (positions(b) - positions(a))) + if intersections: + lookup[cell] = len(vertices) + cells.append(cell) + vertices.append(np.mean(intersections, axis=0)) + triangles = [] + for axis in range(3): + u, v = (axis + 1) % 3, (axis + 2) % 3 + ranges = [range(n + 1) for n in shape] + ranges[axis] = range(shape[axis]) + for corner in product(*ranges): + a = np.array(corner) + b = a.copy() + b[axis] += 1 + if (scalar(a) <= 0) == (scalar(b) <= 0): + continue + ring = [] + for du, dv in ((-1, -1), (0, -1), (0, 0), (-1, 0)): + c = a.copy() + c[u] += du + c[v] += dv + if tuple(c) in lookup: + ring.append(lookup[tuple(c)]) + if len(ring) != 4: + continue + for triangle in ((ring[0], ring[1], ring[2]), (ring[0], ring[2], ring[3])): + p, q, r = [vertices[i] for i in triangle] + if np.dot(np.cross(q - p, r - p), normal) < 0: + triangle = triangle[::-1] + triangles.append(triangle) + return SimpleNamespace(vertices=np.asarray(vertices).reshape(-1, 3), + edges=np.asarray(triangles, dtype=int).reshape(-1, 3)), np.asarray(cells).reshape(-1, 3) + + +def assert_closed(mesh): + report = mesh.capping_report + assert report["watertight"], report + assert report["closure_success"], report + assert report["added_triangles"] == report["added_cap_triangles"] + report["added_transition_triangles"] + assert report["cap_max_plane_error"] == 0 + assert not report["warnings"], report + assert len({tuple(sorted(t)) for t in mesh.edges}) == len(mesh.edges) + triangles = mesh.vertices[mesh.edges] + assert np.all(np.linalg.norm(np.cross(triangles[:, 1] - triangles[:, 0], + triangles[:, 2] - triangles[:, 0]), axis=1) > 0) + + +def test_boundary_lattice_shared_points_and_outward_faces(): + shape, extent = (2, 3, 4), (-7, 2, 10, 16, -2, 6) + coordinates, positions, metadata = boundary_lattice(shape, extent) + assert coordinates.dtype.kind == "i" + assert len(coordinates) == 3 * 4 * 5 - 1 * 2 * 3 + assert len(np.unique(coordinates, axis=0)) == len(coordinates) + assert metadata["triangles"].shape == (4 * (2 * 3 + 3 * 4 + 4 * 2), 3) + center = np.asarray(extent).reshape(3, 2).mean(axis=1) + for triangle in positions[metadata["triangles"]]: + a, b, c = triangle + assert np.dot(np.cross(b - a, c - a), triangle.mean(axis=0) - center) > 0 + np.testing.assert_array_equal(positions.min(axis=0), np.asarray(extent)[::2]) + np.testing.assert_array_equal(positions.max(axis=0), np.asarray(extent)[1::2]) + + +@pytest.mark.parametrize("normal,level", [((1, .37, -.21), .731), ((1, 0, 0), .43), + ((1, 1, 1), .71), ((-.3, 1, .7), .52)]) +def test_plane_closed_and_original_indices_preserved(normal, level): + shape, extent = (4, 5, 3), (0, 1, 0, 1, 0, 1) + mesh, cells = plane_mesh(shape, extent, normal, level) + original_vertices, original_triangles = mesh.vertices.copy(), mesh.edges.copy() + # Neither the QEF rows nor boundary sample order is assumed sorted. + permutation = np.arange(len(cells))[::-1] + inverse = np.argsort(permutation) + mesh.vertices = mesh.vertices[permutation] + mesh.edges = inverse[mesh.edges] + coordinates, positions, _ = boundary_lattice(shape, extent) + assert cap_mesh(mesh, cells[permutation], shape, extent, coordinates[::-1], + (positions @ normal)[::-1], level) is mesh + assert_closed(mesh) + assert mesh.capping_report["open_edges_before"] > 0 + assert mesh.capping_report["non_extent_open_edges_before"] == 0 + assert mesh.capping_report["open_edges_after"] == 0 + assert mesh.capping_report["non_extent_open_edges_after"] == 0 + assert mesh.capping_report["added_cap_triangles"] > 0 + assert mesh.capping_report["added_transition_triangles"] > 0 + np.testing.assert_array_equal(mesh.vertices[:len(cells)], original_vertices[permutation]) + np.testing.assert_array_equal(mesh.edges[:len(original_triangles)], inverse[original_triangles]) + new_triangles = mesh.edges[len(original_triangles):] + assert np.count_nonzero(np.all(new_triangles >= len(cells), axis=1)) == mesh.capping_report["added_cap_triangles"] + assert np.count_nonzero(np.any(new_triangles < len(cells), axis=1)) == mesh.capping_report["added_transition_triangles"] + assert np.all(mesh.vertices >= 0) and np.all(mesh.vertices <= 1) + # Shared box edges have one vertex per position, not one per face. + new = mesh.vertices[len(cells):] + assert len(np.unique(new, axis=0)) == len(new) + t = mesh.vertices[mesh.edges] + assert np.sum(np.einsum("ij,ij->i", t[:, 0], np.cross(t[:, 1], t[:, 2]))) > 0 + + +def test_exact_extent_nonunit_plane(): + shape, extent = (3, 4, 5), (-2, 4, 10, 14, -8, -3) + normal, level = (1, .31, .17), 3.19 + mesh, cells = plane_mesh(shape, extent, normal, level) + coordinates, positions, _ = boundary_lattice(shape, extent) + cap_mesh(mesh, cells, shape, extent, coordinates, positions @ normal, level) + assert_closed(mesh) + np.testing.assert_array_equal(mesh.vertices.min(axis=0), [-2, 10, -8]) + np.testing.assert_array_equal(mesh.vertices.max(axis=0)[1:], [14, -3]) + + +@pytest.mark.parametrize("level", [1., 1.5]) +def test_corner_ties_are_shared_and_limitations_reported(level): + shape, extent = (4, 4, 4), (0, 1, 0, 1, 0, 1) + mesh, cells = plane_mesh(shape, extent, (1, 1, 1), level) + coordinates, positions, _ = boundary_lattice(shape, extent) + with warnings.catch_warnings(record=True) as emitted: + warnings.simplefilter("always") + cap_mesh(mesh, cells, shape, extent, coordinates, positions.sum(axis=1), level) + new = mesh.vertices[len(cells):] + assert len(np.unique(new, axis=0)) == len(new) + assert mesh.capping_report["added_triangles"] > 0 + assert mesh.capping_report["watertight"] or (emitted and mesh.capping_report["warnings"]) + + +@pytest.mark.parametrize("inside", [False, True]) +def test_no_boundary_contour_adds_no_shell(inside): + # An enclosed sphere represented by an outward-oriented octahedron. + vertices = np.vstack((np.eye(3), -np.eye(3))) * .2 + .5 + triangles = [] + for x, y, z in product((0, 3), (1, 4), (2, 5)): + a, b, c = vertices[[x, y, z]] + t = [x, y, z] + if np.dot(np.cross(b - a, c - a), a - .5) < 0: + t.reverse() + triangles.append(t) + mesh = SimpleNamespace(vertices=vertices.copy(), edges=np.asarray(triangles)) + coordinates, positions, _ = boundary_lattice((10, 10, 10), (0, 1, 0, 1, 0, 1)) + values = np.linalg.norm(positions - .5, axis=1) - .2 + cap_mesh(mesh, np.floor(vertices * 10).astype(int), (10, 10, 10), (0, 1, 0, 1, 0, 1), + coordinates, -values if inside else values, 0) + np.testing.assert_array_equal(mesh.vertices, vertices) + np.testing.assert_array_equal(mesh.edges, triangles) + assert mesh.capping_report["added_triangles"] == 0 + assert mesh.capping_report["contour_edges"] == 0 + assert_closed(mesh) + + +def test_missing_escaped_qef_and_internal_hole_are_not_repaired(): + shape, extent = (5, 5, 5), (0, 1, 0, 1, 0, 1) + mesh, cells = plane_mesh(shape, extent, (1, .37, -.21), .73) + mesh.vertices[0] = [-2, -2, -2] + # Remove an interior triangle, not a boundary feature. + mesh.edges = mesh.edges[1:].copy() + original = mesh.edges.copy() + coordinates, positions, _ = boundary_lattice(shape, extent) + with pytest.warns(RuntimeWarning, match="best effort"): + cap_mesh(mesh, cells, shape, extent, coordinates, positions @ (1, .37, -.21), .73) + report = mesh.capping_report + assert report["escaped_qef"] == 1 + assert report["missing_qef"] > 0 + assert report["boundary_edges"] > 0 + assert not report["watertight"] + np.testing.assert_array_equal(mesh.edges[:len(original)], original) + assert all(np.any(t >= len(cells)) for t in mesh.edges[len(original):]) + + +def test_no_available_qef_discards_all_cap_components(): + mesh = SimpleNamespace(vertices=np.empty((0, 3)), edges=np.empty((0, 3), dtype=int)) + coordinates, positions, _ = boundary_lattice((2, 2, 2), (0, 1, 0, 1, 0, 1)) + with pytest.warns(RuntimeWarning, match="missing_qef"): + cap_mesh(mesh, np.empty((0, 3)), (2, 2, 2), (0, 1, 0, 1, 0, 1), + coordinates, positions[:, 0], .3) + assert mesh.capping_report["discarded_components"] > 0 + assert mesh.vertices.shape == (0, 3) + assert mesh.edges.shape == (0, 3) + + +def test_nonfinite_boundary_is_diagnostic_not_geometry(): + mesh = SimpleNamespace(vertices=np.empty((0, 3)), edges=np.empty((0, 3), dtype=int)) + coordinates, positions, _ = boundary_lattice((1, 1, 1), (0, 1, 0, 1, 0, 1)) + values = positions[:, 0].copy() + values[0] = np.nan + with pytest.warns(RuntimeWarning, match="Nonfinite"): + cap_mesh(mesh, np.empty((0, 3)), (1, 1, 1), (0, 1, 0, 1, 0, 1), coordinates, values, .3) + assert mesh.capping_report["unsupported_topology"] == 1 + assert mesh.capping_report["added_triangles"] == 0 + + +def test_only_component_connected_to_available_qef_is_retained(): + shape, extent = (2, 2, 2), (0, 1, 0, 1, 0, 1) + coordinates, positions, _ = boundary_lattice(shape, extent) + values = np.ones(len(coordinates)) + values[np.all(coordinates == 0, axis=1) | np.all(coordinates == 2, axis=1)] = -1 + mesh = SimpleNamespace(vertices=np.asarray([[.1, .1, .1]]), + edges=np.empty((0, 3), dtype=int)) + with pytest.warns(RuntimeWarning, match="missing_qef"): + cap_mesh(mesh, np.asarray([[0, 0, 0]]), shape, extent, coordinates, values, 0) + assert mesh.capping_report["discarded_components"] == 1 + assert mesh.capping_report["watertight"] + assert not mesh.capping_report["closure_success"] + assert np.all(mesh.vertices <= .25) + + +def test_internal_hole_without_missing_qef_is_preserved(): + shape, extent = (6, 6, 6), (0, 1, 0, 1, 0, 1) + mesh, cells = plane_mesh(shape, extent, (1, .37, -.21), .731) + interior = np.flatnonzero(np.all((cells[mesh.edges] > 0) & (cells[mesh.edges] < 5), axis=(1, 2))) + assert len(interior) + removed = mesh.edges[interior[0]].copy() + mesh.edges = np.delete(mesh.edges, interior[0], axis=0) + coordinates, positions, _ = boundary_lattice(shape, extent) + with pytest.warns(RuntimeWarning, match="boundary_edges: 3"): + cap_mesh(mesh, cells, shape, extent, coordinates, positions @ (1, .37, -.21), .731) + assert mesh.capping_report["boundary_edges"] == 3 + assert mesh.capping_report["missing_qef"] == 0 + assert mesh.capping_report["escaped_qef"] == 0 + assert mesh.capping_report["non_extent_open_edges_before"] == 3 + assert mesh.capping_report["non_extent_open_edges_after"] == 3 + assert not mesh.capping_report["closure_success"] + assert tuple(sorted(removed)) not in {tuple(sorted(t)) for t in mesh.edges} + + +def test_prebuilt_geometry_and_single_batched_strict_crossing(monkeypatch): + shape, extent = (3, 4, 5), (0, 1, 0, 1, 0, 1) + geometry = boundary_lattice(shape, extent) + coordinates, positions, metadata = geometry + original_coordinates, original_positions = coordinates.copy(), positions.copy() + original_triangles, original_cells = metadata["triangles"].copy(), metadata["cells"].copy() + mesh, cells = plane_mesh(shape, extent, (1, 0, 0), .43) + # Every QEF on this plane is inset from the actual extent planes, including + # the rim. Its owning cell, not its position, identifies an extent opening. + assert np.all((mesh.vertices > 0) & (mesh.vertices < 1)) + original_crossing = _extent_capping.scalar_crossing_parameters + calls = [] + + def crossing(start, end, iso, *, xp=np): + calls.append(start.shape) + return original_crossing(start, end, iso, xp=xp) + + def unexpected_lattice(*args): + pytest.fail("prebuilt boundary geometry must not be regenerated") + + monkeypatch.setattr(_extent_capping, "boundary_lattice", unexpected_lattice) + monkeypatch.setattr(_extent_capping, "scalar_crossing_parameters", crossing) + # Tiny but nonzero scalar differences must not get tolerance-snapped. + values = (positions[:, 0] - .43) * 1e-200 + cap_mesh(mesh, cells, shape, extent, coordinates, values, 0, boundary_geometry=geometry) + assert calls == [metadata["triangles"].shape] + assert_closed(mesh) + assert mesh.capping_report["open_edges_before"] > 0 + assert mesh.capping_report["non_extent_open_edges_before"] == 0 + new = mesh.vertices[len(cells):] + assert np.any(np.isclose(new[:, 0], .43, atol=1e-15)) + np.testing.assert_array_equal(coordinates, original_coordinates) + np.testing.assert_array_equal(positions, original_positions) + np.testing.assert_array_equal(metadata["triangles"], original_triangles) + np.testing.assert_array_equal(metadata["cells"], original_cells) + + +def test_crossing_arithmetic_overflow_is_reported(): + shape, extent = (1, 1, 1), (0, 1, 0, 1, 0, 1) + coordinates, positions, _ = boundary_lattice(shape, extent) + values = np.where(positions[:, 0] == 0, -1e308, 1e308) + mesh = SimpleNamespace(vertices=np.empty((0, 3)), edges=np.empty((0, 3), dtype=int)) + with pytest.warns(RuntimeWarning, match="Nonfinite scalar crossing interpolation arithmetic"): + cap_mesh(mesh, np.empty((0, 3)), shape, extent, coordinates, values, 0) + assert mesh.capping_report["unsupported_topology"] == 1 + assert not mesh.capping_report["closure_success"] + assert mesh.capping_report["added_triangles"] == 0 + + +@pytest.mark.parametrize("defect,field", [ + ("nonfinite", "nonfinite_vertices"), + ("negative_index", "invalid_triangles"), + ("out_of_range", "invalid_triangles"), + ("fractional_index", "invalid_triangles"), + ("nonfinite_index", "invalid_triangles"), + ("degenerate", "degenerate_triangles"), + ("escaped", "escaped_qef"), + ("unsupported", "unsupported_topology"), +]) +def test_original_geometry_defects_never_report_success(defect, field): + shape, extent = (3, 3, 3), (0, 1, 0, 1, 0, 1) + vertices = np.asarray([[.2, .2, .2], [.7, .2, .2], [.2, .7, .2], [.2, .2, .7]]) + cells = np.floor(vertices * 3).astype(int) + triangles = np.asarray([[0, 2, 1], [0, 1, 3], [0, 3, 2], [1, 2, 3]]) + if defect == "nonfinite": + vertices[0, 0] = np.nan + elif defect == "negative_index": + triangles[0, 0] = -1 + elif defect == "out_of_range": + triangles[0, 0] = len(vertices) + elif defect in ("fractional_index", "nonfinite_index"): + triangles = triangles.astype(float) + triangles[0, 0] = .5 if defect == "fractional_index" else np.nan + elif defect == "degenerate": + vertices[0] = vertices[1] + elif defect == "escaped": + vertices[0, 0] = -.1 + elif defect == "unsupported": + cells[0] = cells[1] + mesh = SimpleNamespace(vertices=vertices.copy(), edges=triangles.copy()) + coordinates, positions, _ = boundary_lattice(shape, extent) + with pytest.warns(RuntimeWarning, match=field): + cap_mesh(mesh, cells, shape, extent, coordinates, np.ones(len(positions)), 0) + report = mesh.capping_report + assert report[field] > 0 + assert not report["closure_success"] + assert report["added_triangles"] == 0 + np.testing.assert_array_equal(mesh.vertices, vertices) + np.testing.assert_array_equal(mesh.edges, triangles) + if defect in ("nonfinite", "degenerate", "escaped", "unsupported"): + assert report["watertight"] # Topological closure alone is insufficient. diff --git a/tests/test_common/test_modules/test_extent_capping_integration.py b/tests/test_common/test_modules/test_extent_capping_integration.py new file mode 100644 index 00000000..9eaa5dca --- /dev/null +++ b/tests/test_common/test_modules/test_extent_capping_integration.py @@ -0,0 +1,413 @@ +"""Extent capping through the production model and boundary evaluation APIs.""" + +import importlib +import json +import warnings +from copy import deepcopy +from types import SimpleNamespace +from unittest.mock import Mock + +import numpy as np +import pytest + +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 TensorsStructure +from gempy_engine.core.data.engine_grid import EngineGrid +from gempy_engine.core.data.input_data_descriptor import InputDataDescriptor +from gempy_engine.core.data.interpolation_functions import CustomInterpolationFunctions +from gempy_engine.core.data.interpolation_input import InterpolationInput +from gempy_engine.core.data.kernel_classes.orientations import Orientations +from gempy_engine.core.data.kernel_classes.surface_points import SurfacePoints +from gempy_engine.core.data.options import InterpolationOptions +from gempy_engine.core.data.options.evaluation_options import MeshExtentCapping +from gempy_engine.core.data.regular_grid import RegularGrid +from gempy_engine.core.data.stack_relation_type import StackRelationType +from gempy_engine.core.data.stacks_structure import StacksStructure + + +dc = importlib.import_module("gempy_engine.API.dual_contouring.extent_capping") +EXTENT = np.array([-2., 4., 10., 14., -8., -3.]) + + +@pytest.fixture(params=[AvailableBackends.numpy, AvailableBackends.PYTORCH], ids=["numpy", "pytorch-cpu"]) +def backend(request, monkeypatch): + torch = pytest.importorskip("torch") if request.param is AvailableBackends.PYTORCH else None + grad_enabled = torch.is_grad_enabled() if torch is not None else None + saved = dict(engine_backend=BackendTensor.engine_backend, use_gpu=BackendTensor.use_gpu, + use_pykeops=BackendTensor.use_pykeops, dtype=BackendTensor.dtype, + grads=BackendTensor.COMPUTE_GRADS) + saved_pykeops = BackendTensor.pykeops_enabled + monkeypatch.setenv("GEMPY_SKIP_TRIANGULATION", "0") + try: + BackendTensor._change_backend(request.param, use_gpu=False, use_pykeops=False, dtype="float64") + BackendTensor.pykeops_enabled = False + yield request.param + finally: + BackendTensor._change_backend(**saved) + BackendTensor.COMPUTE_GRADS = saved["grads"] + BackendTensor.pykeops_enabled = saved_pykeops + if torch is not None: + torch.set_grad_enabled(grad_enabled) + + +@pytest.fixture +def plane_model(backend): + def make(levels=(-5.7,), normal=(0., 0., 1.), extent=EXTENT, resolution=(3, 3, 3)): + functions = CustomInterpolationFunctions( + scalar_field_at_surface_points=np.asarray(levels), + implicit_function=lambda xyz: xyz[:, 0] * normal[0] + xyz[:, 1] * normal[1] + xyz[:, 2] * normal[2], + gx_function=lambda xyz: xyz[:, 0] * 0 + normal[0], + gy_function=lambda xyz: xyz[:, 0] * 0 + normal[1], + gz_function=lambda xyz: xyz[:, 0] * 0 + normal[2], + ) + structure = InputDataDescriptor( + TensorsStructure(number_of_points_per_surface=np.array([], dtype=int)), + StacksStructure( + number_of_points_per_stack=np.array([0]), + number_of_orientations_per_stack=np.array([0]), + number_of_surfaces_per_stack=np.array([len(levels)]), + masking_descriptor=[StackRelationType.BASEMENT], + interp_functions_per_stack=[functions], + ), + ) + inputs = InterpolationInput( + SurfacePoints(np.empty((0, 3))), + Orientations(np.empty((0, 3)), np.empty((0, 3))), + EngineGrid(octree_grid=RegularGrid(np.array(extent, dtype=float), list(resolution))), + np.arange(len(levels) + 1), + ) + options = InterpolationOptions.from_args(10., 1.) + options.evaluation_options.number_octree_levels = 2 + options.evaluation_options.number_octree_levels_surface = 2 + return inputs, options, structure + + return make + + +@pytest.mark.parametrize("mode", list(MeshExtentCapping)) +def test_capping_option_json_roundtrip(mode): + options = InterpolationOptions.from_args(10., 1.) + assert options.evaluation_options.mesh_extraction_extent_capping is MeshExtentCapping.NONE + options.evaluation_options.mesh_extraction_extent_capping = mode + serialized = options.model_dump_json() + assert json.loads(serialized)["evaluation_options"]["mesh_extraction_extent_capping"] == mode.value + restored = InterpolationOptions.model_validate_json(serialized) + assert restored.evaluation_options.mesh_extraction_extent_capping is mode + + +def test_default_and_explicit_none_have_identical_arrays(plane_model, monkeypatch): + boundary = Mock(side_effect=AssertionError("Disabled capping must not evaluate the boundary")) + monkeypatch.setattr(dc, "_interp_on_boundary", boundary) + inputs, options, descriptor = plane_model() + default = compute_model(inputs, options, descriptor) + inputs, options, descriptor = plane_model() + options.evaluation_options.mesh_extraction_extent_capping = MeshExtentCapping.NONE + explicit = compute_model(inputs, options, descriptor) + assert len(default.dc_meshes) == len(explicit.dc_meshes) == 1 + np.testing.assert_array_equal(default.dc_meshes[0].vertices, explicit.dc_meshes[0].vertices) + np.testing.assert_array_equal(default.dc_meshes[0].edges, explicit.dc_meshes[0].edges) + assert default.dc_meshes[0].capping_report is None + for first, second in zip(default.octrees_output, explicit.octrees_output): + np.testing.assert_array_equal(BackendTensor.t.to_numpy(first.grid.octree_grid.values), + BackendTensor.t.to_numpy(second.grid.octree_grid.values)) + np.testing.assert_array_equal(BackendTensor.t.to_numpy(first.outputs[0].exported_fields.scalar_field), + BackendTensor.t.to_numpy(second.outputs[0].exported_fields.scalar_field)) + boundary.assert_not_called() + + +@pytest.mark.parametrize("mode", list(MeshExtentCapping)) +def test_all_modes_preserve_exact_extent_and_caller_grid(plane_model, backend, mode): + inputs, options, descriptor = plane_model() + grid = inputs.grid + root = grid.octree_grid + original_extent = BackendTensor.t.to_numpy(root.orthogonal_extent).copy() + original_values = BackendTensor.t.to_numpy(root.values).copy() + options.evaluation_options.mesh_extraction_extent_capping = mode + solution = compute_model(inputs, options, descriptor) + assert len(solution.octrees_output) == 2 + root_values = BackendTensor.t.to_numpy(solution.octrees_output[0].grid.octree_grid.values) + assert root_values.dtype == np.float64 + coordinate_tolerance = 2 * np.finfo(root_values.dtype).eps * np.max(np.abs(EXTENT)) + for level in solution.octrees_output: + octree = level.grid.octree_grid + np.testing.assert_array_equal(BackendTensor.t.to_numpy(octree.orthogonal_extent), EXTENT) + if backend is AvailableBackends.PYTORCH: + assert octree.orthogonal_extent.device.type == "cpu" + coordinates = BackendTensor.t.to_numpy(octree.integer_coordinates) + shape = BackendTensor.t.to_numpy(octree.regular_grid_shape) + expected = EXTENT[::2] + (coordinates + .5) * (EXTENT[1::2] - EXTENT[::2]) / shape + np.testing.assert_allclose(BackendTensor.t.to_numpy(octree.values), expected, + rtol=0, atol=coordinate_tolerance) + np.testing.assert_array_equal( + BackendTensor.t.to_numpy(solution.octrees_output[0].grid.octree_grid.orthogonal_extent), EXTENT) + assert inputs.grid is grid + assert grid.octree_grid is root + np.testing.assert_array_equal(BackendTensor.t.to_numpy(root.orthogonal_extent), EXTENT) + np.testing.assert_array_equal(BackendTensor.t.to_numpy(root.orthogonal_extent), original_extent) + np.testing.assert_array_equal(BackendTensor.t.to_numpy(root.values), original_values) + + +def test_enabled_planes_are_closed_with_analytic_volume_and_reuse_boundary(plane_model, monkeypatch): + levels = (-4.7, -6.3) + inputs, options, descriptor = plane_model(levels) + options.evaluation_options.mesh_extraction_extent_capping = MeshExtentCapping.SCALAR_LESS_EQUAL + boundary = Mock(wraps=dc._interp_on_boundary) + monkeypatch.setattr(dc, "_interp_on_boundary", boundary) + solution = compute_model(inputs, options, descriptor) + boundary.assert_called_once() + assert len(solution.dc_meshes) == len(levels) + volumes = [] + for index, (mesh, height) in enumerate(zip(solution.dc_meshes, levels)): + report = mesh.capping_report + assert report["watertight"], report + assert report["closure_success"], report + assert not report["warnings"], report + assert report["open_edges_before"] > 0 + assert report["open_edges_after"] == 0 + assert report["cap_max_plane_error"] == 0 + assert report["added_cap_triangles"] > 0 + assert report["added_transition_triangles"] > 0 + assert report["boundary_scalar_points"] == len(boundary.call_args.args[0]) + assert (mesh.stack_index, mesh.surface_index, mesh.exported_surface_index) == (0, index, index) + assert mesh.isovalue == height + assert mesh.inside_convention == "scalar <= isovalue" + np.testing.assert_array_equal(mesh.vertices.min(axis=0), EXTENT[::2]) + np.testing.assert_array_equal(mesh.vertices.max(axis=0)[:2], EXTENT[1:4:2]) + assert np.all(mesh.vertices <= EXTENT[1::2]) + triangles = mesh.vertices[mesh.edges] + edges = np.concatenate([mesh.edges[:, [0, 1]], mesh.edges[:, [1, 2]], mesh.edges[:, [2, 0]]]) + _, counts = np.unique(np.sort(edges, axis=1), axis=0, return_counts=True) + np.testing.assert_array_equal(counts, 2) + # Shift the origin to avoid cancellation for a translated model box. + triangles = triangles - EXTENT[::2] + volume = np.einsum("ij,ij->i", triangles[:, 0], np.cross(triangles[:, 1], triangles[:, 2])).sum() / 6 + volumes.append(volume) + np.testing.assert_allclose(volumes, 6 * 4 * (np.asarray(levels) + 8), rtol=1e-10, atol=0) + + +@pytest.mark.parametrize("fail", [False, True], ids=["success", "failure"]) +def test_boundary_chunks_restore_grid_stack_and_options(plane_model, monkeypatch, fail): + inputs, options, _ = plane_model() + descriptor = SimpleNamespace(stack_structure=SimpleNamespace(n_stacks=2, stack_number=1)) + saved_grid = inputs.grid + options.evaluation_options.evaluation_chunk_size = 3 + options.evaluation_options.compute_scalar = False + options.evaluation_options.compute_scalar_gradient = True + saved_options = options.model_dump_json() + points = np.arange(21, dtype=float).reshape(7, 3) + seen = [] + + def interpolate(actual_input, boundary_options, actual_descriptor): + assert actual_input is inputs + assert actual_descriptor is descriptor + assert boundary_options is not options + assert boundary_options.evaluation_options.compute_scalar is True + assert boundary_options.evaluation_options.compute_scalar_gradient is False + xyz = actual_input.grid.custom_grid.values + seen.append(BackendTensor.t.to_numpy(xyz).copy()) + descriptor.stack_structure.stack_number = 0 + if fail and len(seen) == 2: + raise RuntimeError("boundary evaluation failed") + # Include sentinel values outside the custom-grid slice. + return [SimpleNamespace( + grid=SimpleNamespace(custom_grid_slice=slice(1, len(xyz) + 1)), + exported_fields=SimpleNamespace(scalar_field=BackendTensor.t.concatenate([ + BackendTensor.t.array([-999.]), xyz[:, 2] + offset, BackendTensor.t.array([-999.]) + ])), + ) for offset in (0, 100)] + + interpolation = Mock(side_effect=interpolate) + monkeypatch.setattr(dc, "interpolate_all_fields_no_octree", interpolation) + if fail: + with pytest.raises(RuntimeError, match="boundary evaluation failed"): + dc._interp_on_boundary(points, inputs, options, descriptor) + else: + scalars = dc._interp_on_boundary(points, inputs, options, descriptor) + assert len(scalars) == 2 + for actual, offset in zip(scalars, (0, 100)): + np.testing.assert_array_equal(actual, points[:, 2] + offset) + assert interpolation.call_count == (2 if fail else 3) + np.testing.assert_array_equal(np.concatenate(seen), points[:6] if fail else points) + assert inputs.grid is saved_grid + assert descriptor.stack_structure.stack_number == 1 + assert options.model_dump_json() == saved_options + + +def assert_report_matches_mesh(mesh, emitted): + """Best-effort closure must not conceal topology or geometry defects.""" + report = mesh.capping_report + assert report is not None + assert np.isfinite(mesh.vertices).all() + assert mesh.edges.ndim == 2 and mesh.edges.shape[1] == 3 + assert np.all((mesh.edges >= 0) & (mesh.edges < len(mesh.vertices))) + edges = np.concatenate([mesh.edges[:, [0, 1]], mesh.edges[:, [1, 2]], mesh.edges[:, [2, 0]]]) + _, inverse, counts = np.unique(np.sort(edges, axis=1), axis=0, return_inverse=True, return_counts=True) + assert report["open_edges_after"] == np.count_nonzero(counts == 1), report + assert report["nonmanifold_edges"] == np.count_nonzero(counts > 2), report + directions = np.bincount(inverse, weights=np.where(edges[:, 0] < edges[:, 1], 1, -1)) + assert report["orientation_conflicts"] == np.count_nonzero((counts == 2) & (np.abs(directions) == 2)), report + triangles = mesh.vertices[mesh.edges] + cross = np.cross(triangles[:, 1] - triangles[:, 0], triangles[:, 2] - triangles[:, 0]) + assert report["degenerate_triangles"] == np.count_nonzero(np.all(cross == 0, axis=1)), report + defects = ("missing_qef", "escaped_qef", "unsupported_topology", "boundary_edges", + "nonmanifold_edges", "orientation_conflicts", "nonmanifold_vertices", + "invalid_triangles", "nonfinite_vertices", "nonfinite_triangles", "degenerate_triangles") + for key in defects: + if report[key]: + assert not report["closure_success"], report + assert f"{key}: {report[key]}" in report["warnings"], report + if report["warnings"]: + assert any(issubclass(w.category, RuntimeWarning) + and all(message in str(w.message) for message in report["warnings"]) + for w in emitted), report + if report["closure_success"]: + assert report["watertight"] and not report["warnings"], report + elif len(mesh.edges): + assert report["warnings"], report + if report["watertight"]: + assert len(counts) > 0 + np.testing.assert_array_equal(counts, 2) + assert report["added_triangles"] == report["added_cap_triangles"] + report["added_transition_triangles"] + assert report["cap_max_plane_error"] == 0 + + +@pytest.mark.parametrize("normal,level,volume", [ + ((-.2, -.1, 1.), .313, .463), + ((.2, .1, -1.), -.313, .537), +]) +def test_oblique_plane_production_volume(plane_model, normal, level, volume): + inputs, options, descriptor = plane_model((level,), normal, (0, 1, 0, 1, 0, 1)) + options.evaluation_options.mesh_extraction_extent_capping = MeshExtentCapping.SCALAR_LESS_EQUAL + with warnings.catch_warnings(record=True) as emitted: + warnings.simplefilter("always") + solution = compute_model(inputs, options, descriptor) + assert len(solution.dc_meshes) == 1 + mesh = solution.dc_meshes[0] + assert_report_matches_mesh(mesh, emitted) + assert mesh.capping_report["closure_success"], mesh.capping_report + original = mesh.vertices[:mesh.capping_report["original_vertex_count"]] + np.testing.assert_allclose(original @ normal, level, rtol=0, atol=2e-13) + assert np.all(mesh.vertices >= -2e-13) and np.all(mesh.vertices <= 1 + 2e-13) + triangles = mesh.vertices[mesh.edges] + actual_volume = np.einsum("ij,ij->i", triangles[:, 0], np.cross(triangles[:, 1], triangles[:, 2])).sum() / 6 + assert actual_volume == pytest.approx(volume, rel=1e-11) + + +@pytest.mark.parametrize("normal,level", [ + ((0., 0., 1.), .5), + ((1., 1., 1.), 1.), + ((1., 1., 1.), 1.5), +], ids=["lattice-plane", "box-corners", "lattice-points"]) +def test_exact_plane_ties_have_honest_diagnostics(plane_model, normal, level): + inputs, options, descriptor = plane_model((level,), normal, (0, 1, 0, 1, 0, 1), (2, 2, 2)) + options.evaluation_options.mesh_extraction_extent_capping = MeshExtentCapping.SCALAR_LESS_EQUAL + with warnings.catch_warnings(record=True) as emitted: + warnings.simplefilter("always") + solution = compute_model(inputs, options, descriptor) + assert len(solution.dc_meshes) == 1 + mesh = solution.dc_meshes[0] + assert len(mesh.vertices) > 0 and len(mesh.edges) > 0 + assert_report_matches_mesh(mesh, emitted) + report = mesh.capping_report + assert report["closure_success"] or report["warnings"], report + added = mesh.vertices[report["original_vertex_count"]:] + assert len(np.unique(added, axis=0)) == len(added) + + +@pytest.mark.parametrize("height", [-9., -2.], ids=["fully-outside", "fully-inside"]) +def test_nonintersecting_plane_produces_no_artificial_box(plane_model, height): + inputs, options, descriptor = plane_model((height,)) + options.evaluation_options.mesh_extraction_extent_capping = MeshExtentCapping.SCALAR_LESS_EQUAL + with warnings.catch_warnings(record=True) as emitted: + warnings.simplefilter("always") + solution = compute_model(inputs, options, descriptor) + assert len(solution.dc_meshes) == 1 + mesh = solution.dc_meshes[0] + assert mesh.vertices.shape == (0, 3) + assert mesh.edges.shape == (0, 3) + assert_report_matches_mesh(mesh, emitted) + assert mesh.capping_report["added_vertices"] == 0 + assert mesh.capping_report["added_triangles"] == 0 + assert not mesh.capping_report["closure_success"] + + +@pytest.mark.parametrize("skip", ["1", "true"]) +def test_skip_triangulation_skips_capping(plane_model, monkeypatch, skip): + inputs, options, descriptor = plane_model() + options.evaluation_options.mesh_extraction_extent_capping = MeshExtentCapping.SCALAR_LESS_EQUAL + monkeypatch.setenv("GEMPY_SKIP_TRIANGULATION", skip) + boundary = Mock(side_effect=AssertionError("Skipped triangulation must not evaluate caps")) + cap = Mock(side_effect=AssertionError("Skipped triangulation must not add triangles")) + monkeypatch.setattr(dc, "_interp_on_boundary", boundary) + monkeypatch.setattr(dc, "cap_mesh", cap) + solution = compute_model(inputs, options, descriptor) + assert len(solution.dc_meshes) == 1 + mesh = solution.dc_meshes[0] + assert len(mesh.vertices) > 0 + assert mesh.edges.size == 0 + assert mesh.capping_report is None + boundary.assert_not_called() + cap.assert_not_called() + + +@pytest.mark.parametrize("backend", [AvailableBackends.numpy], indirect=True) +@pytest.mark.parametrize("normal,level", [((0., 0., 1.), .43), ((-.2, -.1, 1.), .313)]) +def test_numpy_pytorch_mesh_equivalence(plane_model, normal, level): + torch = pytest.importorskip("torch") + saved_grad = torch.is_grad_enabled() + meshes = [] + try: + for engine in (AvailableBackends.numpy, AvailableBackends.PYTORCH): + BackendTensor._change_backend(engine, use_gpu=False, use_pykeops=False, + dtype="float64", grads=saved_grad) + inputs, options, descriptor = plane_model((level,), normal, (0, 1, 0, 1, 0, 1)) + options.evaluation_options.mesh_extraction_extent_capping = MeshExtentCapping.SCALAR_LESS_EQUAL + meshes.append(compute_model(inputs, options, descriptor).dc_meshes[0]) + first, second = meshes + assert first.capping_report["closure_success"], first.capping_report + assert second.capping_report["closure_success"], second.capping_report + np.testing.assert_array_equal(first.edges, second.edges) + np.testing.assert_allclose(first.vertices, second.vertices, rtol=0, atol=2e-13) + assert first.capping_report == second.capping_report + finally: + torch.set_grad_enabled(saved_grad) + + +def test_fault_model_capping_reports_all_limitations(graben_fault_model, monkeypatch): + inputs, descriptor, options = deepcopy(graben_fault_model) + options.evaluation_options.number_octree_levels = 3 + options.evaluation_options.number_octree_levels_surface = 3 + options.evaluation_options.mesh_extraction_extent_capping = MeshExtentCapping.SCALAR_LESS_EQUAL + monkeypatch.setenv("GEMPY_SKIP_TRIANGULATION", "0") + with warnings.catch_warnings(record=True) as emitted: + warnings.simplefilter("always") + solution = compute_model(inputs, options, descriptor) + assert len(solution.dc_meshes) == 6 + assert {mesh.stack_index for mesh in solution.dc_meshes} == {0, 1, 2} + for mesh in solution.dc_meshes: + assert_report_matches_mesh(mesh, emitted) + assert mesh.capping_report["boundary_scalar_points"] == 0 + assert mesh.capping_report["added_triangles"] == 0 + assert "boundary ownership" in mesh.capping_report["skipped_reason"] + assert mesh.inside_convention == "scalar <= isovalue" + + +def test_capped_export_offsets_do_not_mutate_meshes(plane_model, monkeypatch): + raw_module = importlib.import_module("gempy_engine.core.data.raw_arrays_solution") + monkeypatch.setattr(raw_module, "require_subsurface", lambda: SimpleNamespace( + UnstructuredData=SimpleNamespace(from_array=lambda **kwargs: kwargs))) + monkeypatch.setattr(raw_module, "require_pandas", lambda: SimpleNamespace(DataFrame=lambda data: data)) + inputs, options, descriptor = plane_model((-4.7, -6.3)) + options.evaluation_options.mesh_extraction_extent_capping = MeshExtentCapping.SCALAR_LESS_EQUAL + meshes = compute_model(inputs, options, descriptor).dc_meshes + raw = raw_module.RawArraysSolution(vertices=[mesh.vertices for mesh in meshes], + edges=[mesh.edges for mesh in meshes]) + original = [mesh.edges.copy() for mesh in meshes] + expected = np.concatenate((original[0], original[1] + len(meshes[0].vertices))) + for _ in range(2): + exported = raw.meshes_to_subsurface() + np.testing.assert_array_equal(exported["cells"], expected) + for mesh, triangles in zip(meshes, original): + np.testing.assert_array_equal(mesh.edges, triangles) diff --git a/tests/test_common/test_modules/test_grid_zero_handling.py b/tests/test_common/test_modules/test_grid_zero_handling.py new file mode 100644 index 00000000..3a852c9d --- /dev/null +++ b/tests/test_common/test_modules/test_grid_zero_handling.py @@ -0,0 +1,67 @@ +import numpy as np +import pytest + +from gempy_engine.config import AvailableBackends +from gempy_engine.core.backend_tensor import BackendTensor +from gempy_engine.core.data.dual_contouring_data import DualContouringData +from gempy_engine.core.data.engine_grid import EngineGrid +from gempy_engine.core.data.regular_grid import RegularGrid +from gempy_engine.modules.dual_contouring._gen_vertices import generate_dual_contouring_vertices +from gempy_engine.modules.octrees_topology._octree_common import _generate_next_level_centers + + +@pytest.fixture(params=['numpy', 'PYTORCH']) +def backend(request): + if request.param == 'PYTORCH': + pytest.importorskip('torch') + old = BackendTensor.engine_backend, BackendTensor.use_gpu, BackendTensor.dtype, BackendTensor.use_pykeops + BackendTensor._change_backend(engine_backend=AvailableBackends[request.param], use_gpu=False, dtype='float64') + yield BackendTensor.t + BackendTensor._change_backend(engine_backend=old[0], use_gpu=old[1], dtype=old[2], use_pykeops=old[3]) + + +def test_grid_extent_and_centers_are_not_translated(backend): + t = backend + extent = [-1., 1., -2., 2., -3., 3.] + dense = RegularGrid(extent, [3, 3, 3]) + np.testing.assert_array_equal(t.to_numpy(dense.values[13]), [0., 0., 0.]) + grid = EngineGrid.from_regular_grid(dense) + assert grid.dense_grid is dense + root = grid.octree_grid + np.testing.assert_array_equal(t.to_numpy(root.orthogonal_extent), extent) + np.testing.assert_array_equal(t.to_numpy(dense.orthogonal_extent), extent) + np.testing.assert_array_equal(t.to_numpy(root.values[0]), [-0.5, -1., -1.5]) + for _ in range(3): + active = root.integer_coordinates[:, 0] == 0 + xyz, bits = _generate_next_level_centers(root.values[active], root.dxdydz) + root = RegularGrid.from_octree_level(xyz, root, active, bits) + np.testing.assert_array_equal(t.to_numpy(root.orthogonal_extent), extent) + expected = root.orthogonal_extent[::2] + (root.integer_coordinates + 0.5) * ( + (root.orthogonal_extent[1::2] - root.orthogonal_extent[::2]) / root.regular_grid_shape + ) + np.testing.assert_array_equal(t.to_numpy(root.values), t.to_numpy(expected)) + + +@pytest.mark.parametrize('strict_crossings', [False, True]) +@pytest.mark.parametrize('points', [ + [[0., 0., 0.], [2., 0., 4.]], + [[0., 0., 0.], [0., 0., 0.]], +]) +def test_mass_points_include_zero_coordinates(backend, strict_crossings, points): + t = backend + valid = t.zeros((2, 12), dtype=bool) + valid[1, 0] = valid[1, 5] = True + data = DualContouringData( + xyz_on_edge=t.array(points, dtype=BackendTensor.dtype_obj), valid_edges=valid, + xyz_on_centers=t.zeros((2, 3)), dxdydz=(1., 1., 1.), + n_surfaces_to_export=1, left_right_codes=None, + gradients=t.array([[0., 1., 0.], [0., 1., 0.]], dtype=BackendTensor.dtype_obj), + strict_crossings=strict_crossings, + ) + vertices = generate_dual_contouring_vertices(data, debug=True) + expected = np.mean(points, axis=0) + np.testing.assert_array_equal(t.to_numpy(data.bias_center_mass), np.tile(expected, (3, 1))) + # Preserve the existing non-strict PyTorch origin regularization. + if BackendTensor.engine_backend == AvailableBackends.PYTORCH and not strict_crossings: + expected = expected / 1.0001 + np.testing.assert_allclose(t.to_numpy(vertices), expected[None, :], atol=1e-15) diff --git a/tests/test_common/test_modules/test_scalar_crossing.py b/tests/test_common/test_modules/test_scalar_crossing.py new file mode 100644 index 00000000..e892f70d --- /dev/null +++ b/tests/test_common/test_modules/test_scalar_crossing.py @@ -0,0 +1,111 @@ +import itertools + +import numpy as np +import pytest + +from gempy_engine.core.backend_tensor import BackendTensor, AvailableBackends +from gempy_engine.core.data.dual_contouring_data import DualContouringData +from gempy_engine.modules.dual_contouring._gen_vertices import generate_dual_contouring_vertices +from gempy_engine.modules.dual_contouring._scalar_crossing import scalar_crossing_parameters + + +def test_crossing_snaps_near_endpoints(): + crossing, parameters = scalar_crossing_parameters( + np.array([-1e-14, -1.]), np.array([1., 1e-14]), 0. + ) + np.testing.assert_array_equal(crossing, [True, True]) + np.testing.assert_array_equal(parameters, [0., 1.]) +from gempy_engine.modules.dual_contouring.dual_contouring_interface import find_intersection_on_edge + + +@pytest.fixture(params=[AvailableBackends.numpy, AvailableBackends.PYTORCH]) +def backend(request): + if request.param == AvailableBackends.PYTORCH: + pytest.importorskip("torch") + BackendTensor._change_backend(request.param, use_gpu=False) + return BackendTensor.t + + +def as_numpy(value): + return value.detach().cpu().numpy() if hasattr(value, "detach") else value + + +def test_parameters(backend): + start = backend.array([-1., 1., 0., 1., 0., -1., 0., 1., -1., -1e-14]) + end = backend.array([1., -1., 1., 0., -1., 0., 0., 2., -1., 3e-14]) + iso = backend.array(0.) + valid, t = scalar_crossing_parameters(start, end, iso, xp=backend) + np.testing.assert_array_equal(as_numpy(valid), [1, 1, 1, 1, 0, 0, 0, 0, 0, 1]) + np.testing.assert_allclose(as_numpy(t), [.5, .5, 0, 1, 0, 0, 0, 0, 0, .25]) + reverse_valid, reverse_t = scalar_crossing_parameters(end, start, iso, xp=backend) + np.testing.assert_array_equal(as_numpy(reverse_valid), as_numpy(valid)) + np.testing.assert_allclose(as_numpy(reverse_t)[as_numpy(valid)], 1 - as_numpy(t)[as_numpy(valid)]) + + +@pytest.mark.parametrize("bad", [np.nan, np.inf, -np.inf]) +@pytest.mark.parametrize("argument", [0, 1, 2]) +def test_nonfinite(backend, bad, argument): + args = [backend.array([-1.]), backend.array([1.]), backend.array([0.])] + args[argument] = backend.array([bad]) + with pytest.raises(ValueError, match="finite"): + scalar_crossing_parameters(*args, xp=backend) + + +def test_numpy_utility_independent_of_backend(backend): + valid, t = scalar_crossing_parameters(np.array([-1.]), np.array([3.]), np.array([0.])) + assert isinstance(t, np.ndarray) + np.testing.assert_array_equal(valid, [True]) + np.testing.assert_array_equal(t, [.25]) + + +def test_intersections_masking_and_multiple_surfaces(backend): + cube = np.array(list(itertools.product([0., 1.], repeat=3))) + xyz = backend.array(np.concatenate([cube, cube + 2])) + scalars = xyz[:, 0] + iso = backend.array([0., .25, 1.]) + points, valid = find_intersection_on_edge( + xyz, scalars, iso, masking=backend.array([True, False]), strict_crossings=True + ) + expected_mask = np.zeros((3, 12), dtype=bool) + expected_mask[:2, :4] = True + np.testing.assert_array_equal(as_numpy(valid).reshape(3, 12), expected_mask) + expected = np.concatenate([cube[:4], cube[:4] + [.25, 0, 0]]) + np.testing.assert_allclose(as_numpy(points), expected) + assert valid.ndim == (1 if BackendTensor.engine_backend == AvailableBackends.PYTORCH else 2) + + +def test_empty_and_iso_edges(backend): + xyz = backend.array(np.zeros((8, 3))) + for masking in (None, backend.array([False])): + points, valid = find_intersection_on_edge( + xyz, xyz[:, 0], backend.array([0.]), masking=masking, strict_crossings=True + ) + assert points.shape == (0, 3) + assert not as_numpy(valid).any() + + +def test_default_retains_tolerant_crossings(backend): + xyz = backend.array(np.array(list(itertools.product([0., 1.], repeat=3)))) + iso = backend.array([1.005]) + default = find_intersection_on_edge(xyz, xyz[:, 0], iso) + disabled = find_intersection_on_edge(xyz, xyz[:, 0], iso, strict_crossings=False) + for actual, expected in zip(default, disabled): + np.testing.assert_array_equal(as_numpy(actual), as_numpy(expected)) + assert as_numpy(default[1]).any() + strict = find_intersection_on_edge(xyz, xyz[:, 0], iso, strict_crossings=True) + assert not as_numpy(strict[1]).any() + + +def test_mass_point_includes_zero_coordinates(backend): + valid = np.zeros((2, 12), dtype=bool) + valid[0, :2] = True + data = DualContouringData( + xyz_on_edge=backend.array([[0., 0., 0.], [2., 4., 0.]], dtype=BackendTensor.dtype_obj), + valid_edges=backend.array(valid), xyz_on_centers=None, dxdydz=(1., 1., 1.), + n_surfaces_to_export=1, left_right_codes=None, + gradients=backend.array(np.zeros((2, 3)), dtype=BackendTensor.dtype_obj), strict_crossings=True + ) + vertices = generate_dual_contouring_vertices(data, debug=True) + np.testing.assert_allclose(as_numpy(data.bias_center_mass), [[1., 2., 0.]] * 3) + assert np.isfinite(as_numpy(vertices)).all() + assert DualContouringData.__dataclass_fields__["strict_crossings"].default is False