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 new file mode 100644 index 00000000..50daab9b --- /dev/null +++ b/tests/unit/modelling/test_scalar_field_grid.py @@ -0,0 +1,96 @@ +"""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, 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. +""" +import numpy as np +import pytest + +from LoopStructural.geometry import BoundingBox +from LoopStructural.modelling.features._analytical_feature import ( + AnalyticalGeologicalFeature, +) + + +@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]) + + +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) + + +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) + + +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) + + +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)