fix: updating bounding box project/reproject to just translate to origin

Author: lachlangrose

LoopStructural/datatypes/_bounding_box.py

@@ -17,6 +17,7 @@ def __init__(
         origin: Optional[np.ndarray] = None,
         maximum: Optional[np.ndarray] = None,
         global_origin: Optional[np.ndarray] = None,
+        global_maximum: Optional[np.ndarray] = None,
         nsteps: Optional[np.ndarray] = None,
         step_vector: Optional[np.ndarray] = None,
         dimensions: Optional[int] = 3,
@@ -39,9 +40,12 @@ def __init__(
         # we want the local coordinates to start at 0
         # otherwise uses provided origin. This is useful for having multiple bounding boxes rela
         if global_origin is not None and origin is None:
-            origin = np.zeros(global_origin.shape)
+            origin = np.zeros(np.array(global_origin).shape)
+        if global_maximum is not None and global_origin is not None:
+            maximum = np.array(global_maximum) - np.array(global_origin)
+        
         if maximum is None and nsteps is not None and step_vector is not None:
-            maximum = origin + nsteps * step_vector
+            maximum = np.array(origin) + np.array(nsteps) * np.array(step_vector)
         if origin is not None and global_origin is None:
             global_origin = np.zeros(3)
         self._origin = np.array(origin)
@@ -517,11 +521,8 @@ def project(self, xyz, inplace=False):
         """
         if inplace:
             xyz -= self.global_origin
-            xyz = np.clip(xyz, self.origin, self.maximum)
             return xyz
-        return (xyz - self.global_origin) / np.max(
-            (self.global_maximum - self.global_origin)
-        )  # np.clip(xyz, self.origin, self.maximum)
+        return (xyz - self.global_origin)  # np.clip(xyz, self.origin, self.maximum)
 
     def scale_by_projection_factor(self, value):
         return value / np.max((self.global_maximum - self.global_origin))
@@ -541,10 +542,9 @@ def reproject(self, xyz, inplace=False):
             reprojected point
         """
         if inplace:
-            xyz -= self.global_origin
-            xyz /= np.max((self.global_maximum - self.global_origin))
+            xyz += self.global_origin
             return xyz
-        return xyz * np.max((self.global_maximum - self.global_origin)) + self.global_origin
+        return xyz + self.global_origin
 
     def __repr__(self):
         return f"BoundingBox({self.origin}, {self.maximum}, {self.nsteps})"

LoopStructural/modelling/core/geological_model.py

@@ -111,8 +111,8 @@ def __init__(
                 raise ValueError("Must provide a bounding box")
             self.bounding_box = bounding_box
         if len(args) == 2:
-            origin = args[0]
-            maximum = args[1]
+            origin = np.array(args[0])
+            maximum = np.array(args[1])
             if not isinstance(origin, np.ndarray) or not isinstance(maximum, np.ndarray):
                 raise ValueError("Must provide origin and maximum as numpy arrays")
             self.bounding_box = BoundingBox(
@@ -510,7 +510,7 @@ def data(self, data: pd.DataFrame):
 
     def set_model_data(self, data):
         logger.warning("deprecated method. Model data can now be set using the data attribute")
-        self.data = data
+        self.data = data.copy()
 
     def set_stratigraphic_column(self, stratigraphic_column, cmap="tab20"):
         """

tests/integration/test_interpolator.py

@@ -14,8 +14,8 @@ def model_fit(model, data):
 def test_create_model():
     data, bb = load_claudius()
     model = GeologicalModel(bb[0, :], bb[1, :])
-    assert np.all(np.isclose(model.origin, bb[0, :]))
-    assert np.all(np.isclose(model.maximum, bb[1, :]))
+    assert np.all(np.isclose(model.bounding_box.global_origin, bb[0, :]))
+    assert np.all(np.isclose(model.bounding_box.global_maximum, bb[1, :]))
 
 
 def test_add_data():

tests/unit/datatypes/test_bounding_box.py

@@ -97,7 +97,11 @@ def test_regular_grid_2d():
     grid = bbox.regular_grid((10, 10))
     assert grid.shape == (10 * 10, 2)
 
-
+def test_project_to_local():
+    bbox = BoundingBox(global_origin=[10,10,10], global_maximum=[20,20,20])
+    point = np.array([15, 15, 15])
+    local_point = bbox.project(point)
+    assert np.all(local_point == np.array([5, 5, 5]))
 if __name__ == "__main__":
     test_create_bounding_box()
     test_create_bounding_box_from_points()

tests/unit/modelling/test_geological_model.py

@@ -1,17 +1,30 @@
 from LoopStructural import GeologicalModel
 from LoopStructural.datasets import load_claudius
 import numpy as np
+import pytest
 
+@pytest.mark.parametrize("origin, maximum", [([0,0,0],[5,5,5]), ([10,10,10],[15,15,15])])
+def test_create_geological_model(origin, maximum):
+    model = GeologicalModel(origin, maximum)
+    assert (model.bounding_box.global_origin - np.array(origin)).sum() == 0
+    assert (model.bounding_box.global_maximum - np.array(maximum)).sum() == 0
+    assert (model.bounding_box.origin - np.zeros(3)).sum() == 0
+    assert (model.bounding_box.maximum - np.ones(3)*5).sum() == 0
 
-def test_create_geological_model():
-    model = GeologicalModel([0, 0, 0], [5, 5, 5])
-    assert (model.origin - np.array([0, 0, 0])).sum() == 0
-    assert (model.maximum - np.array([5, 5, 5])).sum() == 0
-
-
+def test_rescale_model_data():
+    data, bb = load_claudius()
+    model = GeologicalModel(bb[0, :], bb[1, :])
+    model.set_model_data(data)
+    # Check that the model data is rescaled to local coordinates
+    expected = data[['X', 'Y', 'Z']].values - bb[None, 0, :]
+    actual = model.data[['X', 'Y', 'Z']].values
+    assert np.allclose(actual, expected, atol=1e-6)
 def test_access_feature_model():
     data, bb = load_claudius()
     model = GeologicalModel(bb[0, :], bb[1, :])
     model.set_model_data(data)
     s0 = model.create_and_add_foliation("strati")
     assert s0 == model["strati"]
+
+if __name__ == "__main__":
+    test_rescale_model_data()