From 03406c08807672c2cf8be1a3810d7d9f275eb9b1 Mon Sep 17 00:00:00 2001 From: benjamink Date: Tue, 10 Feb 2026 14:20:49 -0800 Subject: [PATCH 1/6] mask topo surface before active cells compute --- simpeg_drivers/components/topography.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/simpeg_drivers/components/topography.py b/simpeg_drivers/components/topography.py index fadcf2a6..394e51f9 100644 --- a/simpeg_drivers/components/topography.py +++ b/simpeg_drivers/components/topography.py @@ -101,13 +101,13 @@ def active_cells(self, mesh: InversionMesh, data: InversionData) -> np.ndarray: self.params.active_cells.active_model, mesh.entity ) else: + topo_surface = self.params.active_cells.topography_object.copy() + topo_surface.mask_by_extent(mesh.entity.extent) active_cells = active_from_xyz( mesh.entity, self.locations, grid_reference="bottom" if forced_to_surface else "center", - triangulation=getattr( - self.params.active_cells.topography_object, "cells", None - ), + triangulation=getattr(topo_surface, "cells", None), ) active_cells = (mesh.permutation @ active_cells).astype(bool) From 23219b8fe8c675ccd892a18aeb9558b22b23a756 Mon Sep 17 00:00:00 2001 From: benjamink Date: Wed, 11 Feb 2026 08:49:00 -0800 Subject: [PATCH 2/6] restrict topography cells to the mesh extents when calculating the active cells --- simpeg_drivers/components/topography.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/simpeg_drivers/components/topography.py b/simpeg_drivers/components/topography.py index 394e51f9..294522db 100644 --- a/simpeg_drivers/components/topography.py +++ b/simpeg_drivers/components/topography.py @@ -101,13 +101,22 @@ def active_cells(self, mesh: InversionMesh, data: InversionData) -> np.ndarray: self.params.active_cells.active_model, mesh.entity ) else: - topo_surface = self.params.active_cells.topography_object.copy() - topo_surface.mask_by_extent(mesh.entity.extent) + centers = self.params.active_cells.topography_object.centroids + xmin, xmax, ymin, ymax, zmin, zmax = mesh.entity.extent.ravel(order="F") + mask = ( + (centers[:, 0] > xmin) + & (centers[:, 0] < xmax) + & (centers[:, 1] > ymin) + & (centers[:, 1] < ymax) + & (centers[:, 2] > zmin) + & (centers[:, 2] < zmax) + ) + cells = self.params.active_cells.topography_object.cells[mask, :] active_cells = active_from_xyz( mesh.entity, self.locations, grid_reference="bottom" if forced_to_surface else "center", - triangulation=getattr(topo_surface, "cells", None), + triangulation=cells, ) active_cells = (mesh.permutation @ active_cells).astype(bool) From a68968418ec6cea6c425028e3316d80a37b10647 Mon Sep 17 00:00:00 2001 From: benjamink Date: Wed, 11 Feb 2026 09:58:16 -0800 Subject: [PATCH 3/6] switching to dev --- simpeg_drivers/components/topography.py | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/simpeg_drivers/components/topography.py b/simpeg_drivers/components/topography.py index 294522db..acf83a5a 100644 --- a/simpeg_drivers/components/topography.py +++ b/simpeg_drivers/components/topography.py @@ -101,22 +101,24 @@ def active_cells(self, mesh: InversionMesh, data: InversionData) -> np.ndarray: self.params.active_cells.active_model, mesh.entity ) else: - centers = self.params.active_cells.topography_object.centroids + topography = self.params.active_cells.topography_object + locations = getattr(topography, "centroids", None) or self.locations xmin, xmax, ymin, ymax, zmin, zmax = mesh.entity.extent.ravel(order="F") mask = ( - (centers[:, 0] > xmin) - & (centers[:, 0] < xmax) - & (centers[:, 1] > ymin) - & (centers[:, 1] < ymax) - & (centers[:, 2] > zmin) - & (centers[:, 2] < zmax) + (locations[:, 0] > xmin) + & (locations[:, 0] < xmax) + & (locations[:, 1] > ymin) + & (locations[:, 1] < ymax) + & (locations[:, 2] > zmin) + & (locations[:, 2] < zmax) ) - cells = self.params.active_cells.topography_object.cells[mask, :] + + cells = getattr(topography, "cells", None) active_cells = active_from_xyz( mesh.entity, - self.locations, + locations, grid_reference="bottom" if forced_to_surface else "center", - triangulation=cells, + triangulation=cells[mask, :] or None, ) active_cells = (mesh.permutation @ active_cells).astype(bool) From 61f3311571715d98d2a6dde7b0cfb48f04318013 Mon Sep 17 00:00:00 2001 From: benjamink Date: Thu, 12 Feb 2026 09:36:15 -0800 Subject: [PATCH 4/6] tests passing --- simpeg_drivers/components/topography.py | 31 +++++++------ simpeg_drivers/utils/utils.py | 59 ++++++++++++++++++++++++- 2 files changed, 75 insertions(+), 15 deletions(-) diff --git a/simpeg_drivers/components/topography.py b/simpeg_drivers/components/topography.py index acf83a5a..832a0165 100644 --- a/simpeg_drivers/components/topography.py +++ b/simpeg_drivers/components/topography.py @@ -38,6 +38,8 @@ floating_active, get_containing_cells, get_neighbouring_cells, + mask_vertices_and_cells, + octree_extents, ) @@ -100,25 +102,26 @@ def active_cells(self, mesh: InversionMesh, data: InversionData) -> np.ndarray: active_cells = InversionModel.obj_2_mesh( self.params.active_cells.active_model, mesh.entity ) + else: - topography = self.params.active_cells.topography_object - locations = getattr(topography, "centroids", None) or self.locations - xmin, xmax, ymin, ymax, zmin, zmax = mesh.entity.extent.ravel(order="F") - mask = ( - (locations[:, 0] > xmin) - & (locations[:, 0] < xmax) - & (locations[:, 1] > ymin) - & (locations[:, 1] < ymax) - & (locations[:, 2] > zmin) - & (locations[:, 2] < zmax) - ) + if any(k in self.params.inversion_type for k in ["2d", "p3d"]): + vertices = self.locations + cells = getattr( + self.params.active_cells.topography_object, "cells", None + ) + else: + extent = octree_extents(mesh.entity) + vertices, cells = mask_vertices_and_cells( + extent.ravel(order="F"), + self.locations, + getattr(self.params.active_cells.topography_object, "cells", None), + ) - cells = getattr(topography, "cells", None) active_cells = active_from_xyz( mesh.entity, - locations, + vertices, grid_reference="bottom" if forced_to_surface else "center", - triangulation=cells[mask, :] or None, + triangulation=cells, ) active_cells = (mesh.permutation @ active_cells).astype(bool) diff --git a/simpeg_drivers/utils/utils.py b/simpeg_drivers/utils/utils.py index 64bb1cd9..fdc2bdc2 100644 --- a/simpeg_drivers/utils/utils.py +++ b/simpeg_drivers/utils/utils.py @@ -11,6 +11,7 @@ from __future__ import annotations +from collections.abc import Sequence from copy import deepcopy from typing import TYPE_CHECKING @@ -43,6 +44,62 @@ from simpeg_drivers.driver import InversionDriver +def octree_extents(octree: Octree) -> np.ndarray: + """ + Get the true extents of an octree (min/max of the perimeter). + + The octree.extents property returns min/max of the centroids + + :param octree: Octree mesh object. + + :returns: Array of [xmin, xmax, ymin, ymax]. + """ + + def half_cell(axis): + return ( + getattr(octree, f"{axis}_cell_size") * octree.octree_cells["NCells"] + ) / 2 + + xmin = (octree.centroids[:, 0] - half_cell("u")).min() + xmax = (octree.centroids[:, 0] + half_cell("u")).max() + ymin = (octree.centroids[:, 1] - half_cell("v")).min() + ymax = (octree.centroids[:, 1] + half_cell("v")).max() + + return np.array([xmin, xmax, ymin, ymax]) + + +def mask_vertices_and_cells( + extent: Sequence, vertices: np.ndarray, cells: np.ndarray | None +) -> tuple[np.ndarray, np.ndarray]: + """ + Mask vertices and remove cells whose vertices are all outside the extent. + + :param extent: Array-like object of [xmin, xmax, ymin, ymax]. + :param vertices: Array of shape (n_vertices, 3) containing the x, y, z coordinates. + :param cells: Array of shape (n_cells, 3) containing the indices of the vertices + that make up each cell. + """ + + vertex_mask = ( + (vertices[:, 0] >= extent[0]) + & (vertices[:, 0] <= extent[1]) + & (vertices[:, 1] >= extent[2]) + & (vertices[:, 1] <= extent[3]) + ) + if cells is None: + return vertices[vertex_mask], None + + cell_mask = np.any(vertex_mask[cells], axis=1) + vertex_mask = np.zeros_like(vertex_mask, dtype=bool) + vertex_mask[cells[cell_mask].flatten()] = True + + new_cells = cells.copy()[cell_mask] + cell_map = np.arange(len(vertices))[vertex_mask] + new_cells = np.searchsorted(cell_map, new_cells) + + return vertices[vertex_mask], new_cells + + def calculate_2D_trend( points: np.ndarray, values: np.ndarray, order: int = 0, method: str = "all" ): @@ -495,7 +552,7 @@ def active_from_xyz( raise ValueError("'grid_reference' must be one of 'center', 'top', or 'bottom'") # Return the active cell array - return mask_under_horizon(locations, topo, triangulation=triangulation) + return mask_under_horizon(locations, horizon=topo, triangulation=triangulation) def truncate_locs_depths(locs: np.ndarray, depth_core: float) -> np.ndarray: From fdf7fd00b3c8d8c3639184714d3e29686d3302d8 Mon Sep 17 00:00:00 2001 From: benjamink Date: Fri, 13 Feb 2026 14:08:40 -0800 Subject: [PATCH 5/6] Cleaner octree_extents implementation --- simpeg_drivers/components/topography.py | 2 +- simpeg_drivers/utils/utils.py | 18 ++++--- tests/utils_test.py | 64 +++++++++++++++++++++++++ 3 files changed, 73 insertions(+), 11 deletions(-) create mode 100644 tests/utils_test.py diff --git a/simpeg_drivers/components/topography.py b/simpeg_drivers/components/topography.py index 832a0165..9e3f3663 100644 --- a/simpeg_drivers/components/topography.py +++ b/simpeg_drivers/components/topography.py @@ -110,7 +110,7 @@ def active_cells(self, mesh: InversionMesh, data: InversionData) -> np.ndarray: self.params.active_cells.topography_object, "cells", None ) else: - extent = octree_extents(mesh.entity) + extent = octree_extents(mesh.entity)[:4] vertices, cells = mask_vertices_and_cells( extent.ravel(order="F"), self.locations, diff --git a/simpeg_drivers/utils/utils.py b/simpeg_drivers/utils/utils.py index fdc2bdc2..59996371 100644 --- a/simpeg_drivers/utils/utils.py +++ b/simpeg_drivers/utils/utils.py @@ -55,17 +55,15 @@ def octree_extents(octree: Octree) -> np.ndarray: :returns: Array of [xmin, xmax, ymin, ymax]. """ - def half_cell(axis): - return ( - getattr(octree, f"{axis}_cell_size") * octree.octree_cells["NCells"] - ) / 2 - - xmin = (octree.centroids[:, 0] - half_cell("u")).min() - xmax = (octree.centroids[:, 0] + half_cell("u")).max() - ymin = (octree.centroids[:, 1] - half_cell("v")).min() - ymax = (octree.centroids[:, 1] + half_cell("v")).max() + origin = np.array(list(octree.origin.tolist())) + span = np.array( + [ + getattr(octree, f"{axis}_cell_size") * getattr(octree, f"{axis}_count") + for axis in "uvw" + ] + ) - return np.array([xmin, xmax, ymin, ymax]) + return np.stack([origin, origin + span]).flatten(order="F") def mask_vertices_and_cells( diff --git a/tests/utils_test.py b/tests/utils_test.py new file mode 100644 index 00000000..fc11f220 --- /dev/null +++ b/tests/utils_test.py @@ -0,0 +1,64 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2026 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import numpy as np +from geoh5py import Workspace +from geoh5py.objects import Octree, Points +from grid_apps.octree_creation.driver import OctreeDriver +from grid_apps.octree_creation.options import OctreeOptions, RefinementOptions + +from simpeg_drivers.utils.utils import mask_vertices_and_cells, octree_extents + + +def test_octree_extents(tmp_path): + with Workspace(tmp_path / "test.geoh5") as ws: + X, Y = np.meshgrid(np.linspace(0, 1000, 51), np.linspace(0, 1000, 51)) + Z = np.zeros_like(X) + vertices = np.column_stack([X.ravel(), Y.ravel(), Z.ravel()]) + pts = Points.create(ws, name="points", vertices=vertices) + options = OctreeOptions( + geoh5=ws, + objects=pts, + refinements=[ + RefinementOptions( + refinement_object=pts, levels=[4, 2], horizon=False, distance=100 + ), + ], + ) + octree = OctreeDriver.octree_from_params(options) + + extents = octree_extents(octree) + assert np.allclose(extents, [-1112.5, 2087.5, -1112.5, 2087.5]) + + +def test_mask_vertices_and_cells(): + X, Y = np.meshgrid(np.arange(3), np.arange(3)) + Z = np.zeros_like(X) + vertices = np.column_stack([X.ravel(), Y.ravel(), Z.ravel()]) + cells = np.array( + [ + [0, 1, 3], + [3, 1, 4], + [1, 2, 4], + [4, 2, 5], + [3, 4, 6], + [6, 4, 7], + [4, 5, 7], + [7, 5, 8], + ] + ) + extent = [0.5, 2, 0, 2, 0, 1] + masked_vertices, masked_cells = mask_vertices_and_cells(extent, vertices, cells) + assert len(masked_vertices) == len(vertices) + assert len(masked_cells) == len(cells) + extent = [1.5, 2, 0, 2, 0, 1] + masked_vertices, masked_cells = mask_vertices_and_cells(extent, vertices, cells) + assert len(masked_vertices) == 6 + assert len(masked_cells) == 4 From 785d24626c97b2100c61aca9ddf34f0980d0b636 Mon Sep 17 00:00:00 2001 From: benjamink Date: Fri, 13 Feb 2026 14:17:19 -0800 Subject: [PATCH 6/6] fix tests --- tests/utils_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/utils_test.py b/tests/utils_test.py index fc11f220..fdfb8014 100644 --- a/tests/utils_test.py +++ b/tests/utils_test.py @@ -35,7 +35,7 @@ def test_octree_extents(tmp_path): octree = OctreeDriver.octree_from_params(options) extents = octree_extents(octree) - assert np.allclose(extents, [-1112.5, 2087.5, -1112.5, 2087.5]) + assert np.allclose(extents, [-1112.5, 2087.5, -1112.5, 2087.5, -1062.5, 537.5]) def test_mask_vertices_and_cells():