From 5b5fcc76c78cdee1f874eb434c73ec2fc8052292 Mon Sep 17 00:00:00 2001 From: Lachlan Grose Date: Fri, 16 May 2025 13:18:29 +1000 Subject: [PATCH 1/3] add inplace arg to boundingbox --- LoopStructural/datatypes/_bounding_box.py | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/LoopStructural/datatypes/_bounding_box.py b/LoopStructural/datatypes/_bounding_box.py index b4a822d6..55da13db 100644 --- a/LoopStructural/datatypes/_bounding_box.py +++ b/LoopStructural/datatypes/_bounding_box.py @@ -500,20 +500,25 @@ def structured_grid( name=name, ) - def project(self, xyz): + def project(self, xyz, inplace=False): """Project a point into the bounding box Parameters ---------- xyz : np.ndarray point to project + inplace : bool, optional + Whether to modify the input array in place, by default False Returns ------- np.ndarray projected point """ - + 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) @@ -521,20 +526,24 @@ def project(self, xyz): def scale_by_projection_factor(self, value): return value / np.max((self.global_maximum - self.global_origin)) - def reproject(self, xyz): + def reproject(self, xyz, inplace=False): """Reproject a point from the bounding box to the global space Parameters ---------- xyz : np.ndarray point to reproject - + inplace : bool, optional + Whether to modify the input array in place, by default False Returns ------- np.ndarray reprojected point """ - + if inplace: + xyz -= self.global_origin + xyz /= np.max((self.global_maximum - self.global_origin)) + return xyz return xyz * np.max((self.global_maximum - self.global_origin)) + self.global_origin def __repr__(self): From 829b869a366f240ff9f18b84dfebb516033db999 Mon Sep 17 00:00:00 2001 From: Lachlan Grose Date: Fri, 16 May 2025 13:18:49 +1000 Subject: [PATCH 2/3] move to rely on boundingbox class for boundingbox --- .../modelling/core/geological_model.py | 200 ++++-------------- 1 file changed, 47 insertions(+), 153 deletions(-) diff --git a/LoopStructural/modelling/core/geological_model.py b/LoopStructural/modelling/core/geological_model.py index 72233b09..018de253 100644 --- a/LoopStructural/modelling/core/geological_model.py +++ b/LoopStructural/modelling/core/geological_model.py @@ -2,7 +2,7 @@ Main entry point for creating a geological model """ -from ...utils import getLogger, log_to_file +from ...utils import getLogger import numpy as np import pandas as pd @@ -61,33 +61,24 @@ class GeologicalModel: the origin of the model box parameters : dict a dictionary tracking the parameters used to build the model - scale_factor : double - the scale factor used to rescale the model - + """ def __init__( self, - origin: np.ndarray, - maximum: np.ndarray, - data=None, - nsteps=(50, 50, 25), - reuse_supports=False, - logfile=None, - loglevel="info", + *args ): """ Parameters ---------- - origin : numpy array - specifying the origin of the model - maximum : numpy array - specifying the maximum extent of the model - rescale : bool - whether to rescale the model to between 0/1 - epsion : float - a fudge factor for isosurfacing, used to make sure surfaces appear + bounding_box : BoundingBox + the bounding box of the model + origin : np.array(3,dtype=doubles) + the origin of the model + maximum : np.array(3,dtype=doubles) + the maximum of the model + Examples -------- Demo data @@ -111,38 +102,32 @@ def __init__( """ - if logfile: - self.logfile = logfile - log_to_file(logfile, level=loglevel) - + args = list(args) + if len(args) == 0: + raise ValueError("Must provide either bounding_box or origin and maximum") + if len(args) == 1: + bounding_box = args[0] + if not isinstance(bounding_box, BoundingBox): + raise ValueError("Must provide a bounding box") + self.bounding_box = bounding_box + if len(args) == 2: + origin = args[0] + maximum = 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( + dimensions=3, + origin=np.zeros(3), + maximum=maximum - origin, + global_origin=origin, + ) logger.info("Initialising geological model") self.features = [] self.feature_name_index = {} self._data = pd.DataFrame() # None - if data is not None: - self.data = data - self.nsteps = nsteps - - # we want to rescale the model area so that the maximum length is - # 1 - self.origin = np.array(origin).astype(float) - originstr = f"Model origin: {self.origin[0]} {self.origin[1]} {self.origin[2]}" - logger.info(originstr) - self.maximum = np.array(maximum).astype(float) - maximumstr = "Model maximum: {} {} {}".format( - self.maximum[0], self.maximum[1], self.maximum[2] - ) - logger.info(maximumstr) - - self.scale_factor = 1.0 - - self.bounding_box = BoundingBox( - dimensions=3, - origin=np.zeros(3), - maximum=self.maximum - self.origin, - global_origin=self.origin, - ) + + self.stratigraphic_column = None self.tol = 1e-10 * np.max(self.bounding_box.maximum - self.bounding_box.origin) @@ -160,10 +145,7 @@ def to_dict(self): json = {} json["model"] = {} json["model"]["features"] = [f.name for f in self.features] - # json["model"]["data"] = self.data.to_json() - # json["model"]["origin"] = self.origin.tolist() - # json["model"]["maximum"] = self.maximum.tolist() - # json["model"]["nsteps"] = self.nsteps + json['model']['bounding_box'] = self.bounding_box.to_dict() json["model"]["stratigraphic_column"] = self.stratigraphic_column # json["features"] = [f.to_json() for f in self.features] return json @@ -192,86 +174,12 @@ def to_dict(self): # model.features.append(GeologicalFeature.from_json(feature,model)) # return model def __str__(self): - lengths = self.maximum - self.origin - _str = "GeologicalModel - {} x {} x {}\n".format(*lengths) - _str += "------------------------------------------ \n" - _str += "The model contains {} GeologicalFeatures \n".format(len(self.features)) - _str += "" - _str += "------------------------------------------ \n" - _str += "" - _str += "Model origin: {} {} {}\n".format(self.origin[0], self.origin[1], self.origin[2]) - _str += "Model maximum: {} {} {}\n".format( - self.maximum[0], self.maximum[1], self.maximum[2] - ) - _str += "Model rescale factor: {} \n".format(self.scale_factor) - _str += "------------------------------------------ \n" - _str += "Feature list: \n" - for feature in self.features: - _str += " {} \n".format(feature.name) - return _str + return f"GeologicalModel with {len(self.features)} features" def _ipython_key_completions_(self): return self.feature_name_index.keys() - @classmethod - def from_map2loop_directory( - cls, - m2l_directory, - foliation_params={}, - fault_params={}, - use_thickness=True, - vector_scale=1, - gradient=False, - **kwargs, - ): - """Alternate constructor for a geological model using m2l output - - Uses the information saved in the map2loop files to build a geological model. - You can specify kwargs for building foliation using foliation_params and for - faults using fault_params. faults is a flag that allows for the faults to be - skipped. - - Parameters - ---------- - m2l_directory : string - path to map2loop directory - - Returns - ------- - (GeologicalModel, dict) - the created geological model and a dictionary of the map2loop data - - Notes - ------ - For additional information see :class:`LoopStructural.modelling.input.Map2LoopProcessor` - and :meth:`LoopStructural.GeologicalModel.from_processor` - """ - from LoopStructural.modelling.input.map2loop_processor import Map2LoopProcessor - - log_to_file(f"{m2l_directory}/loopstructural_log.txt") - logger.info("Creating model from m2l directory") - processor = Map2LoopProcessor(m2l_directory, use_thickness) - processor._gradient = gradient - processor.vector_scale = vector_scale - for foliation_name in processor.stratigraphic_column.keys(): - if foliation_name != "faults": - if foliation_name in foliation_params.keys(): - processor.foliation_properties[foliation_name] = foliation_params[ - foliation_name - ] - else: - processor.foliation_properties[foliation_name] = foliation_params - - for fault_name in processor.fault_names: - if fault_name in fault_params.keys(): - for param_name, value in fault_params[fault_name].items(): - processor.fault_properties.loc[fault_name, param_name] = value - else: - for param_name, value in fault_params.items(): - processor.fault_properties.loc[fault_name, param_name] = value - - model = GeologicalModel.from_processor(processor) - return model, processor + @classmethod def from_processor(cls, processor): @@ -565,13 +473,9 @@ def data(self, data: pd.DataFrame): raise BaseException("Cannot load data") logger.info(f"Adding data to GeologicalModel with {len(data)} data points") self._data = data.copy() - - self._data["X"] -= self.origin[0] - self._data["Y"] -= self.origin[1] - self._data["Z"] -= self.origin[2] - self._data["X"] /= self.scale_factor - self._data["Y"] /= self.scale_factor - self._data["Z"] /= self.scale_factor + self._data[['X','Y','Z']] = self.bounding_box.project(self._data[['X','Y','Z']].to_numpy()) + + if "type" in self._data: logger.warning("'type' is deprecated replace with 'feature_name' \n") self._data.rename(columns={"type": "feature_name"}, inplace=True) @@ -1358,7 +1262,7 @@ def create_and_add_fault( if "data_region" in kwargs: kwargs.pop("data_region") logger.error("kwarg data_region currently not supported, disabling") - displacement_scaled = displacement / self.scale_factor + displacement_scaled = displacement fault_frame_builder = FaultBuilder( interpolatortype, bounding_box=self.bounding_box, @@ -1375,11 +1279,11 @@ def create_and_add_fault( if fault_center is not None and ~np.isnan(fault_center).any(): fault_center = self.scale(fault_center, inplace=False) if minor_axis: - minor_axis = minor_axis / self.scale_factor + minor_axis = minor_axis if major_axis: - major_axis = major_axis / self.scale_factor + major_axis = major_axis if intermediate_axis: - intermediate_axis = intermediate_axis / self.scale_factor + intermediate_axis = intermediate_axis fault_frame_builder.create_data_from_geometry( fault_frame_data=fault_frame_data, fault_center=fault_center, @@ -1434,11 +1338,9 @@ def rescale(self, points: np.ndarray, inplace: bool = False) -> np.ndarray: points : np.array((N,3),dtype=double) """ - if not inplace: - points = points.copy() - points *= self.scale_factor - points += self.origin - return points + + return self.bounding_box.reproject(points,inplace=inplace) + # TODO move scale to bounding box/transformer def scale(self, points: np.ndarray, inplace: bool = False) -> np.ndarray: @@ -1456,16 +1358,8 @@ def scale(self, points: np.ndarray, inplace: bool = False) -> np.ndarray: points : np.a::rray((N,3),dtype=double) """ - points = np.array(points).astype(float) - if not inplace: - points = points.copy() - # if len(points.shape) == 1: - # points = points[None,:] - # if len(points.shape) != 2: - # logger.error("cannot scale array of dimensions".format(len(points.shape))) - points -= self.origin - points /= self.scale_factor - return points + return self.bounding_box.project(np.array(points).astype(float),inplace=inplace) + def regular_grid(self, nsteps=None, shuffle=True, rescale=False, order="C"): """ @@ -1614,7 +1508,7 @@ def evaluate_fault_displacements(self, points, scale=True): if f.type == FeatureType.FAULT: disp = f.displacementfeature.evaluate_value(points) vals[~np.isnan(disp)] += disp[~np.isnan(disp)] - return vals * -self.scale_factor # convert from restoration magnutude to displacement + return vals # convert from restoration magnutude to displacement def get_feature_by_name(self, feature_name) -> GeologicalFeature: """Returns a feature from the mode given a name From 13f3118e8293c7a288dfb025882ea0d4759bd404 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Mon, 19 May 2025 14:14:08 +1000 Subject: [PATCH 3/3] fix: updating bounding box project/reproject to just translate to origin --- LoopStructural/datatypes/_bounding_box.py | 18 ++++++------- .../modelling/core/geological_model.py | 6 ++--- tests/integration/test_interpolator.py | 4 +-- tests/unit/datatypes/test_bounding_box.py | 6 ++++- tests/unit/modelling/test_geological_model.py | 25 ++++++++++++++----- 5 files changed, 38 insertions(+), 21 deletions(-) diff --git a/LoopStructural/datatypes/_bounding_box.py b/LoopStructural/datatypes/_bounding_box.py index 55da13db..a2e1eb25 100644 --- a/LoopStructural/datatypes/_bounding_box.py +++ b/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})" diff --git a/LoopStructural/modelling/core/geological_model.py b/LoopStructural/modelling/core/geological_model.py index 018de253..28012269 100644 --- a/LoopStructural/modelling/core/geological_model.py +++ b/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"): """ diff --git a/tests/integration/test_interpolator.py b/tests/integration/test_interpolator.py index bcb326fc..c3c1a046 100644 --- a/tests/integration/test_interpolator.py +++ b/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(): diff --git a/tests/unit/datatypes/test_bounding_box.py b/tests/unit/datatypes/test_bounding_box.py index 1b3bb2cf..ecf9f594 100644 --- a/tests/unit/datatypes/test_bounding_box.py +++ b/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() diff --git a/tests/unit/modelling/test_geological_model.py b/tests/unit/modelling/test_geological_model.py index 0318181a..c09405e9 100644 --- a/tests/unit/modelling/test_geological_model.py +++ b/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() \ No newline at end of file