From 2dd06aaa4092cfc9a49abd0b072a692a3996e5bd Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 09:44:47 +1030 Subject: [PATCH 1/2] test: pin scalar_field grid registration (strict xfail) feature.scalar_field stores values computed at bounding_box.regular_grid() / cell_centres() on a grid built from step_vector and nsteps. Those disagree, so the values are misregistered and the grid stops short of the box. The tests fail today and are marked strict xfail so they flip to XPASS when the geometry is fixed. Co-Authored-By: Claude Sonnet 5.5 --- .../unit/modelling/test_scalar_field_grid.py | 71 +++++++++++++++++++ 1 file changed, 71 insertions(+) create mode 100644 tests/unit/modelling/test_scalar_field_grid.py diff --git a/tests/unit/modelling/test_scalar_field_grid.py b/tests/unit/modelling/test_scalar_field_grid.py new file mode 100644 index 00000000..fc7ecb76 --- /dev/null +++ b/tests/unit/modelling/test_scalar_field_grid.py @@ -0,0 +1,71 @@ +"""The grid returned by ``feature.scalar_field`` must be registered to the values stored on it. + +``scalar_field`` evaluates the feature at ``bounding_box.regular_grid()`` / +``bounding_box.cell_centres()`` and stores the results on a grid built from +``bounding_box.origin``, ``step_vector`` and ``nsteps``. Those two descriptions of the box +disagree (``step_vector`` is ``extent / nsteps`` but ``regular_grid`` places ``nsteps`` points +across the full extent), so the stored values are not at the node positions the grid reports. + +The feature is the analytical plane ``f(x, y, z) = x``, so the right value at any location is +its x coordinate and no interpolation is involved. + +These are strict xfails: they fail today. When the bounding box geometry is made consistent +they will XPASS, which fails the run - remove the markers then. +""" +import numpy as np +import pytest + +from LoopStructural.geometry import BoundingBox +from LoopStructural.modelling.features._analytical_feature import ( + AnalyticalGeologicalFeature, +) + +GEOMETRY_BUG = pytest.mark.xfail( + strict=True, + reason="BoundingBox.step_vector (extent / nsteps) is inconsistent with regular_grid " + "(nsteps points across the full extent), so scalar_field grids are misregistered", +) + + +@pytest.fixture +def feature(): + # value == x coordinate + return AnalyticalGeologicalFeature( + name="x", + vector=np.array([1.0, 0.0, 0.0]), + origin=np.array([0.0, 0.0, 0.0]), + ) + + +@pytest.fixture +def bounding_box(): + return BoundingBox(origin=[0.0, 0.0, 0.0], maximum=[10.0, 10.0, 10.0], nsteps=[10, 10, 10]) + + +@GEOMETRY_BUG +def test_scalar_field_values_are_at_the_node_positions(feature, bounding_box): + grid = feature.scalar_field(bounding_box) + nodes = np.asarray(grid.nodes) + np.testing.assert_allclose(grid.properties["x"], nodes[:, 0], atol=1e-8) + + +@GEOMETRY_BUG +def test_scalar_field_grid_spans_the_bounding_box(feature, bounding_box): + grid = feature.scalar_field(bounding_box) + nodes = np.asarray(grid.nodes) + np.testing.assert_allclose(nodes.min(axis=0), bounding_box.origin, atol=1e-8) + np.testing.assert_allclose(nodes.max(axis=0), bounding_box.maximum, atol=1e-8) + + +@GEOMETRY_BUG +def test_scalar_field_cell_values_are_at_the_cell_centres(feature, bounding_box): + grid = feature.scalar_field(bounding_box) + centres = np.asarray(grid.cell_centres) + np.testing.assert_allclose(grid.cell_properties["x"], centres[:, 0], atol=1e-8) + + +@GEOMETRY_BUG +def test_cell_centres_are_inside_the_bounding_box(bounding_box): + centres = bounding_box.cell_centres() + assert np.all(centres >= bounding_box.origin - 1e-8) + assert np.all(centres <= bounding_box.maximum + 1e-8) From 7f2e94cb88373e4ed190d8c878df5b468aa52e49 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 10:34:27 +1030 Subject: [PATCH 2/2] fix: make scalar_field/block-model grids line up with their values BoundingBox.nsteps is a cell count for interpolation supports (step_vector = extent / nsteps) but regular_grid() puts nsteps points across the full extent. scalar_field, gradient_norm_scalar_field and get_block_model built their result grid from step_vector, so values were stored on nodes up to one step away from where they were computed and the grid stopped short of the box. - add BoundingBox.node_step_vector = extent / (nsteps - 1) and use it for those grids - cell_centres now returns the midpoints of the nsteps-1 cells between the regular_grid points (it overshot the box by half a cell before) - fix the x upper bound of the result grid's cell_centres (+ instead of -) - document the cell-based and node-based conventions on BoundingBox - turn the strict-xfail registration tests into passing tests and add offset / anisotropic box cases Co-Authored-By: Claude Sonnet 5.5 --- COMPAT.md | 1 + .../modelling/core/_model_exporter.py | 2 +- .../features/_base_geological_feature.py | 4 +- .../src/loop_common/geometry/_bounding_box.py | 46 +++++++++++++-- .../loop_common/geometry/_structured_grid.py | 2 +- .../unit/modelling/test_scalar_field_grid.py | 59 +++++++++++++------ 6 files changed, 87 insertions(+), 27 deletions(-) diff --git a/COMPAT.md b/COMPAT.md index 0de642d4..7a5441b0 100644 --- a/COMPAT.md +++ b/COMPAT.md @@ -38,6 +38,7 @@ unchanged, but the new path is preferred going forward. | `LoopStructural.modelling.features._feature_converters.add_fold_to_feature` (private module, imported directly by the QGIS plugin) | `GeologicalModel.add_fold_to_feature` | 2026-07-24 | | `LoopStructural.modelling.features._feature_converters.convert_feature_to_structural_frame` | `GeologicalModel.convert_feature_to_structural_frame` | 2026-07-24 | | `LoopStructural.geometry.StructuredGrid` / `loop_common.geometry.StructuredGrid` (results container with `properties`/`cell_properties`) | `PropertyGrid` (same module; `StructuredGrid` stays as an alias). `LoopStructural.geometry` now re-exports the `loop_common` class instead of carrying a duplicate copy. Requires `loop-common>=0.1.1,<0.2.0`; publish the updated shared package before releasing LoopStructural | 2026-10-29 | +| `BoundingBox.cell_centres()` positions (`regular_grid(nsteps-1) + 0.5*step_vector`, which overshot the box by up to half a cell) | Midpoints of the `nsteps-1` cells between the `regular_grid()` points (always inside the box). Count unchanged. New `BoundingBox.node_step_vector` = `extent/(nsteps-1)`; `scalar_field`, `gradient_norm_scalar_field` and `get_block_model` grids now use it, so their values line up with their node and cell positions | 2026-10-29 | | `GeologicalModel.create_and_add_fault(..., faults=[])` | `GeologicalModel.create_and_add_fault(..., faults=None)` (list built internally, same effective default) | 2026-07-27 | | `GeologicalModel.create_and_add_intrusion(..., intrusion_frame_parameters={}, geometric_scaling_parameters={})` | `GeologicalModel.create_and_add_intrusion(..., intrusion_frame_parameters=None, geometric_scaling_parameters=None)` (dicts built internally, same effective default) | 2026-07-27 | | `GeologicalModel.get_fault_surfaces(faults=[])` | `GeologicalModel.get_fault_surfaces(faults=None)` (list built internally, same effective default) | 2026-07-27 | diff --git a/LoopStructural/modelling/core/_model_exporter.py b/LoopStructural/modelling/core/_model_exporter.py index 8276930e..7098e048 100644 --- a/LoopStructural/modelling/core/_model_exporter.py +++ b/LoopStructural/modelling/core/_model_exporter.py @@ -76,7 +76,7 @@ def get_block_model( # the PropertyGrid result container for storing values. grid = PropertyGrid( origin=model.bounding_box.origin, - step_vector=model.bounding_box.step_vector, + step_vector=model.bounding_box.node_step_vector, nsteps=model.bounding_box.nsteps, name=name, ) diff --git a/LoopStructural/modelling/features/_base_geological_feature.py b/LoopStructural/modelling/features/_base_geological_feature.py index 6e8bfdd1..73956c21 100644 --- a/LoopStructural/modelling/features/_base_geological_feature.py +++ b/LoopStructural/modelling/features/_base_geological_feature.py @@ -389,7 +389,7 @@ def scalar_field(self, bounding_box=None): # the PropertyGrid result container for storing values. grid = PropertyGrid( origin=bounding_box.origin, - step_vector=bounding_box.step_vector, + step_vector=bounding_box.node_step_vector, nsteps=bounding_box.nsteps, name=self.name, ) @@ -419,7 +419,7 @@ def gradient_norm_scalar_field(self, bounding_box=None): bounding_box = self.model.bounding_box grid = PropertyGrid( origin=bounding_box.origin, - step_vector=bounding_box.step_vector, + step_vector=bounding_box.node_step_vector, nsteps=bounding_box.nsteps, name=self.name, ) diff --git a/packages/loop_common/src/loop_common/geometry/_bounding_box.py b/packages/loop_common/src/loop_common/geometry/_bounding_box.py index d4978129..65d584fb 100644 --- a/packages/loop_common/src/loop_common/geometry/_bounding_box.py +++ b/packages/loop_common/src/loop_common/geometry/_bounding_box.py @@ -30,6 +30,20 @@ def __init__( """A bounding box for a model, defined by the origin, maximum and number of steps in each direction + Notes + ----- + ``nsteps`` is used in two ways, so two spacings are provided: + + * Cell-based (interpolation): ``nsteps`` is the number of cells in each + direction and ``step_vector`` is ``(maximum - origin) / nsteps``. Interpolation + supports are built from these. + * Node-based (values on a grid): ``regular_grid`` returns ``nsteps`` points spanning + the full extent, so the spacing between them is ``node_step_vector``, + ``(maximum - origin) / (nsteps - 1)``. ``cell_centres`` returns the + ``nsteps - 1`` centres of the cells between those points. Use this + convention (as pyvista/VTK ``ImageData`` and ``numpy.linspace`` do) + when storing values computed at ``regular_grid`` or ``cell_centres``. + Parameters ---------- dimensions : int, optional @@ -393,6 +407,18 @@ def step_vector(self): raise LoopValueError("Cannot compute step_vector: nsteps contains zero values") return (self.maximum - self.origin) / self.nsteps + @property + def node_step_vector(self): + """Spacing between the ``nsteps`` points returned by :meth:`regular_grid`. + + ``(maximum - origin) / (nsteps - 1)``; compare ``step_vector``, which divides + by ``nsteps`` and is the cell size used to build interpolation supports. + """ + nodes = np.asarray(self.nsteps, dtype=float) + if np.any(nodes < 2): + raise LoopValueError("Cannot compute node_step_vector: nsteps must be at least 2") + return (self.maximum - self.origin) / (nodes - 1) + @property def length(self): return self.maximum - self.origin @@ -573,21 +599,29 @@ def regular_grid( return locs def cell_centres(self, order: str = "F") -> np.ndarray: - """Get the cell centres of a regular grid + """Get the centres of the cells between the points of :meth:`regular_grid` + + There are ``nsteps - 1`` cells in each direction, each ``node_step_vector`` wide, + so every centre is inside the bounding box. Parameters ---------- order : str, optional - order of the grid, by default "C" + order of the grid, by default "F" Returns ------- np.ndarray - array of cell centres + array of cell centres, in the local coordinate system """ - locs = self.regular_grid(order=order, nsteps=self.nsteps - 1) - - return locs + 0.5 * self.step_vector + half = 0.5 * self.node_step_vector + coordinates = [ + np.linspace(self.origin[i] + half[i], self.maximum[i] - half[i], self.nsteps[i] - 1) + for i in range(self.dimensions) + ] + coordinate_grid = np.meshgrid(*coordinates, indexing="ij") + locs = np.array([coord.flatten(order=order) for coord in coordinate_grid]).T + return self.project(locs) def to_dict(self) -> dict: """Export the defining characteristics of the bounding diff --git a/packages/loop_common/src/loop_common/geometry/_structured_grid.py b/packages/loop_common/src/loop_common/geometry/_structured_grid.py index 06a7e8d4..a6622ce5 100644 --- a/packages/loop_common/src/loop_common/geometry/_structured_grid.py +++ b/packages/loop_common/src/loop_common/geometry/_structured_grid.py @@ -131,7 +131,7 @@ def cell_centres(self): """ x = np.linspace( self.origin[0] + self.step_vector[0] * 0.5, - self.maximum[0] + self.step_vector[0] * 0.5, + self.maximum[0] - self.step_vector[0] * 0.5, self.nsteps[0] - 1, ) y = np.linspace( diff --git a/tests/unit/modelling/test_scalar_field_grid.py b/tests/unit/modelling/test_scalar_field_grid.py index fc7ecb76..50daab9b 100644 --- a/tests/unit/modelling/test_scalar_field_grid.py +++ b/tests/unit/modelling/test_scalar_field_grid.py @@ -1,16 +1,13 @@ """The grid returned by ``feature.scalar_field`` must be registered to the values stored on it. ``scalar_field`` evaluates the feature at ``bounding_box.regular_grid()`` / -``bounding_box.cell_centres()`` and stores the results on a grid built from -``bounding_box.origin``, ``step_vector`` and ``nsteps``. Those two descriptions of the box -disagree (``step_vector`` is ``extent / nsteps`` but ``regular_grid`` places ``nsteps`` points -across the full extent), so the stored values are not at the node positions the grid reports. +``bounding_box.cell_centres()`` and stores the results on a grid, so the grid's node +and cell positions must be the points the values were computed at. ``regular_grid`` puts +``nsteps`` points across the full extent (spacing ``node_step_vector``); ``step_vector`` is the +cell size used for interpolation supports and must not be used to describe that grid. The feature is the analytical plane ``f(x, y, z) = x``, so the right value at any location is its x coordinate and no interpolation is involved. - -These are strict xfails: they fail today. When the bounding box geometry is made consistent -they will XPASS, which fails the run - remove the markers then. """ import numpy as np import pytest @@ -20,12 +17,6 @@ AnalyticalGeologicalFeature, ) -GEOMETRY_BUG = pytest.mark.xfail( - strict=True, - reason="BoundingBox.step_vector (extent / nsteps) is inconsistent with regular_grid " - "(nsteps points across the full extent), so scalar_field grids are misregistered", -) - @pytest.fixture def feature(): @@ -42,14 +33,12 @@ def bounding_box(): return BoundingBox(origin=[0.0, 0.0, 0.0], maximum=[10.0, 10.0, 10.0], nsteps=[10, 10, 10]) -@GEOMETRY_BUG def test_scalar_field_values_are_at_the_node_positions(feature, bounding_box): grid = feature.scalar_field(bounding_box) nodes = np.asarray(grid.nodes) np.testing.assert_allclose(grid.properties["x"], nodes[:, 0], atol=1e-8) -@GEOMETRY_BUG def test_scalar_field_grid_spans_the_bounding_box(feature, bounding_box): grid = feature.scalar_field(bounding_box) nodes = np.asarray(grid.nodes) @@ -57,15 +46,51 @@ def test_scalar_field_grid_spans_the_bounding_box(feature, bounding_box): np.testing.assert_allclose(nodes.max(axis=0), bounding_box.maximum, atol=1e-8) -@GEOMETRY_BUG def test_scalar_field_cell_values_are_at_the_cell_centres(feature, bounding_box): grid = feature.scalar_field(bounding_box) centres = np.asarray(grid.cell_centres) np.testing.assert_allclose(grid.cell_properties["x"], centres[:, 0], atol=1e-8) -@GEOMETRY_BUG def test_cell_centres_are_inside_the_bounding_box(bounding_box): centres = bounding_box.cell_centres() assert np.all(centres >= bounding_box.origin - 1e-8) assert np.all(centres <= bounding_box.maximum + 1e-8) + + +@pytest.fixture +def offset_box(): + # non-zero origin and different node counts per axis: catches mixed-up axes and offsets + return BoundingBox(origin=[100.0, -20.0, 5.0], maximum=[112.0, 4.0, 35.0], nsteps=[7, 5, 4]) + + +def test_node_step_vector_is_the_spacing_of_regular_grid(offset_box): + points = offset_box.regular_grid(local=False) + for axis in range(3): + coords = np.unique(np.round(points[:, axis], 10)) + assert len(coords) == offset_box.nsteps[axis] + np.testing.assert_allclose(np.diff(coords), offset_box.node_step_vector[axis]) + + +def test_cell_centres_count_and_midpoints(offset_box): + centres = offset_box.cell_centres() + assert centres.shape == (np.prod(offset_box.nsteps - 1), 3) + for axis in range(3): + coords = np.unique(np.round(centres[:, axis], 10)) + nodes = np.linspace( + offset_box.origin[axis], offset_box.maximum[axis], offset_box.nsteps[axis] + ) + np.testing.assert_allclose(coords, 0.5 * (nodes[:-1] + nodes[1:])) + + +def test_scalar_field_on_offset_anisotropic_box(feature, offset_box): + feature = AnalyticalGeologicalFeature( + name="f", + vector=np.array([1.0, 2.0, 3.0]), + origin=np.array([0.0, 0.0, 0.0]), + ) + grid = feature.scalar_field(offset_box) + expected_nodes = feature.evaluate_value(np.asarray(grid.nodes)) + np.testing.assert_allclose(grid.properties["f"], expected_nodes, atol=1e-8) + expected_cells = feature.evaluate_value(np.asarray(grid.cell_centres)) + np.testing.assert_allclose(grid.cell_properties["f"], expected_cells, atol=1e-8)