Merge branch 'master' of https://github.com/Loop3D/LoopStructural

Author: lachlangrose

.github/workflows/release-please.yml

@@ -9,7 +9,7 @@ jobs:
       fail-fast: false
       matrix:
         os: ["ubuntu-latest",  "windows-latest"] #"macos-latest",
-        python-version: ["3.7","3.8","3.9","3.10"]
+        python-version: ["3.8","3.9","3.10"]
     steps:
       - uses: actions/checkout@v2
       - uses: conda-incubator/setup-miniconda@v2
@@ -87,7 +87,7 @@ jobs:
       fail-fast: false
       matrix:
         os: ["ubuntu-latest",  "windows-latest"]
-        python-version: ["3.10","3.9","3.8","3.7"]
+        python-version: ["3.10","3.9","3.8"]
     steps:
       - uses: conda-incubator/setup-miniconda@v2
         with:

LoopStructural/api/__init__.py

@@ -0,0 +1,145 @@
+from ast import LShift
+import numpy as np
+import pandas as pd
+
+from typing import Optional
+from LoopStructural.interpolators import (
+    GeologicalInterpolator,
+    InterpolatorFactory,
+)
+from LoopStructural.utils import BoundingBox
+from LoopStructural.utils import getLogger
+
+logger = getLogger(__name__)
+
+
+class LoopInterpolator:
+    def __init__(
+        self,
+        bounding_box: BoundingBox,
+        dimensions: int = 3,
+        type: str = "FDI",
+        nelements: int = 1000,
+    ):
+        """Scikitlearn like interface for LoopStructural interpolators
+        useful for quickly building an interpolator to apply to a dataset
+        build a generic interpolation object speficying the bounding box
+        and then fit to constraints and evaluate the interpolator
+
+        Parameters
+        ----------
+        bounding_box : BoundingBox
+            Area where the interpolation will work
+        dimensions : int, optional
+            number of dimensions e.g. 3d or 2d, by default 3
+        type : str, optional
+            type of interpolation algorithm FDI- finite difference, PLI - linear finite elements,
+              by default "FDI"
+        nelements : int, optional
+            degrees of freedom for interpolator, by default 1000
+        """
+        self.dimensions = dimensions
+        self.type = "FDI"
+        self.interpolator: GeologicalInterpolator = (
+            InterpolatorFactory.create_interpolator(
+                type,
+                bounding_box,
+                nelements,
+            )
+        )
+
+    def fit(
+        self,
+        values: Optional[np.ndarray] = None,
+        tangent_vectors: Optional[np.ndarray] = None,
+        normal_vectors: Optional[np.ndarray] = None,
+        inequality_constraints: Optional[np.ndarray] = None,
+    ):
+        """_summary_
+
+        Parameters
+        ----------
+        values : Optional[np.ndarray], optional
+            Value constraints for implicit function, by default None
+        tangent_vectors : Optional[np.ndarray], optional
+            tangent constraints for implicit function, by default None
+        normal_vectors : Optional[np.ndarray], optional
+            gradient norm constraints for implicit function, by default None
+        inequality_constraints : Optional[np.ndarray], optional
+            _description_, by default None
+        """
+        if values is not None:
+            self.interpolator.set_value_constraints(values)
+        if tangent_vectors is not None:
+            self.interpolator.set_tangent_constraints(tangent_vectors)
+        if normal_vectors is not None:
+            self.interpolator.set_normal_constraints(normal_vectors)
+        if inequality_constraints:
+            pass
+
+        self.interpolator.setup()
+
+    def evaluate_scalar_value(self, locations: np.ndarray) -> np.ndarray:
+        """Evaluate the value of the interpolator at locations
+
+        Parameters
+        ----------
+        locations : np.ndarray
+            Nx3 array of locations to evaluate the interpolator at
+
+        Returns
+        -------
+        np.ndarray
+            value of implicit function at locations
+        """
+        self.interpolator.update()
+        return self.interpolator.evaluate_value(locations)
+
+    def evaluate_gradient(self, locations: np.ndarray) -> np.ndarray:
+        """Evaluate the gradient of the interpolator at locations
+
+        Parameters
+        ----------
+        locations : np.ndarray
+            Nx3 locations
+
+        Returns
+        -------
+        np.ndarray
+            Nx3 gradient of implicit function
+        """
+        self.interpolator.update()
+        return self.interpolator.evaluate_gradient(locations)
+
+    def fit_and_evaluate_value(
+        self,
+        values: Optional[np.ndarray] = None,
+        tangent_vectors: Optional[np.ndarray] = None,
+        normal_vectors: Optional[np.ndarray] = None,
+        inequality_constraints: Optional[np.ndarray] = None,
+    ):
+        # get locations
+        self.fit(
+            values=values,
+            tangent_vectors=tangent_vectors,
+            normal_vectors=normal_vectors,
+            inequality_constraints=inequality_constraints,
+        )
+        locations = self.interpolator.get_data_locations()
+        return self.evalute_scalar_value(locations)
+
+    def fit_and_evaluate_gradient(
+        self,
+        values: Optional[np.ndarray] = None,
+        tangent_vectors: Optional[np.ndarray] = None,
+        normal_vectors: Optional[np.ndarray] = None,
+        inequality_constraints: Optional[np.ndarray] = None,
+    ):
+        self.fit(
+            values=values,
+            tangent_vectors=tangent_vectors,
+            normal_vectors=normal_vectors,
+            inequality_constraints=inequality_constraints,
+        )
+        locations = self.interpolator.get_data_locations()
+        return self.evaluate_gradient(locations)

LoopStructural/api/_interpolate.py

@@ -0,0 +1,47 @@
+from typing import Optional, Union
+import numpy as np
+from LoopStructural.utils import getLogger
+
+logger = getLogger(__name__)
+try:
+    from skimage.measure import marching_cubes
+except ImportError:
+    logger.warning("Using deprecated version of scikit-image")
+    from skimage.measure import marching_cubes_lewiner as marching_cubes
+
+from LoopStructural.interpolators import GeologicalInterpolator
+from LoopStructural.utils import BoundingBox
+from LoopStructural.datatypes import Surface
+
+surface_list = dict[str, tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]]
+
+
+class LoopIsosurfacer:
+    def __init__(self, bounding_box: BoundingBox, interpolator: GeologicalInterpolator):
+        self.bounding_box = bounding_box
+        self.interpolator = interpolator
+
+    def fit(self, values: Union[list, int, float]) -> surface_list:
+        surfaces = {}
+        all_values = self.interpolator.evaluate_value(self.bounding_box.regular_grid())
+        if isinstance(values, list):
+            isovalues = values
+        elif isinstance(values, float):
+            isovalues = [values]
+        elif isinstance(values, int):
+            isovalues = np.linspace(np.min(all_values), np.max(all_values), values)
+        for isovalue in isovalues:
+            verts, faces, normals, values = marching_cubes(
+                all_values.reshape(self.bounding_box.nsteps, order="C"),
+                isovalue,
+                spacing=self.bounding_box.step_vector,
+            )
+            values = np.zeros(verts.shape[0]) + isovalue
+            surfaces[f"surface_{isovalue}"] = Surface(
+                vertices=verts,
+                triangles=faces,
+                normals=normals,
+                name=f"surface_{isovalue}",
+                values=values,
+            )
+        return surfaces

LoopStructural/api/_surface.py

@@ -0,0 +1,6 @@
+def CreateModelWithSingleScalarField():
+    pass
+
+
+def AddFaultToModel():
+    pass

LoopStructural/api/model.py

@@ -0,0 +1,2 @@
+from ._surface import Surface
+from ._bounding_box import BoundingBox

LoopStructural/datatypes/__init__.py

@@ -0,0 +1,147 @@
+from __future__ import annotations
+from typing import Optional
+from LoopStructural.utils.exceptions import LoopValueError
+import numpy as np
+
+
+class BoundingBox:
+    def __init__(
+        self,
+        dimensions: int = 3,
+        origin: Optional[np.ndarray] = None,
+        maximum: Optional[np.ndarray] = None,
+        nsteps: Optional[np.ndarray] = None,
+    ):
+        self._origin = origin
+        self._maximum = maximum
+        self.dimensions = dimensions
+        if nsteps is None:
+            self.nsteps = np.array([50, 50, 25])
+        self.name_map = {
+            "xmin": (0, 0),
+            "ymin": (0, 1),
+            "zmin": (0, 2),
+            "xmax": (1, 0),
+            "ymax": (1, 1),
+            "zmax": (1, 2),
+            "lower": (0, 2),
+            "upper": (1, 2),
+            "minx": (0, 0),
+            "miny": (0, 1),
+            "minz": (0, 2),
+            "maxx": (1, 0),
+            "maxy": (1, 1),
+            "maxz": (1, 2),
+        }
+
+    @property
+    def valid(self):
+        return self._origin is not None and self._maximum is not None
+
+    @property
+    def origin(self) -> np.ndarray:
+        if self._origin is None:
+            raise LoopValueError("Origin is not set")
+        return self._origin
+
+    @origin.setter
+    def origin(self, origin: np.ndarray):
+        self._origin = origin
+
+    @property
+    def maximum(self) -> np.ndarray:
+        if self._maximum is None:
+            raise LoopValueError("Maximum is not set")
+        return self._maximum
+
+    @maximum.setter
+    def maximum(self, maximum: np.ndarray):
+        self._maximum = maximum
+
+    @property
+    def nelements(self):
+        return self.nsteps.prod()
+
+    @property
+    def volume(self):
+        return np.product(self.maximum - self.origin)
+
+    @property
+    def bb(self):
+        return np.array([self.origin, self.maximum])
+
+    @nelements.setter
+    def nelements(self, nelements):
+        box_vol = self.volume
+        ele_vol = box_vol / nelements
+        # calculate the step vector of a regular cube
+        step_vector = np.zeros(3)
+        step_vector[:] = ele_vol ** (1.0 / 3.0)
+        # step_vector /= np.array([1,1,2])
+        # number of steps is the length of the box / step vector
+        nsteps = np.ceil((self.maximum - self.origin) / step_vector).astype(int)
+        self.nsteps = nsteps
+
+    @property
+    def step_vector(self):
+        return (self.maximum - self.origin) / self.nsteps
+
+    def fit(self, locations: np.ndarray):
+        if locations.shape[1] != self.dimensions:
+            raise LoopValueError(
+                f"locations array is {locations.shape[1]}D but bounding box is {self.dimensions}"
+            )
+        print("fitting")
+        self.origin = locations.min(axis=0)
+        self.maximum = locations.max(axis=0)
+        return self
+
+    def with_buffer(self, buffer: float = 0.2) -> BoundingBox:
+        if self.origin is None or self.maximum is None:
+            raise LoopValueError(
+                "Cannot create bounding box with buffer, no origin or maximum"
+            )
+        origin = self.origin - buffer * (self.maximum - self.origin)
+        maximum = self.maximum + buffer * (self.maximum - self.origin)
+        return BoundingBox(origin=origin, maximum=maximum)
+
+    def get_value(self, name):
+        ix, iy = self.name_map.get(name, (-1, -1))
+        if ix == -1 and iy == -1:
+            raise LoopValueError(f"{name} is not a valid bounding box name")
+        if iy == -1:
+            return self.origin[ix]
+
+        return self.bb[ix,]
+
+    def __getitem__(self, name):
+        if isinstance(name, str):
+            return self.get_value(name)
+        elif isinstance(name, tuple):
+            return self.origin
+        return self.get_value(name)
+
+    def is_inside(self, xyz):
+        inside = np.zeros(xyz.shape[0], dtype=bool)
+        inside = np.logical_and(inside, xyz[:, 0] > self.origin[0])
+        inside = np.logical_and(inside, xyz[:, 0] < self.maximum[0])
+        inside = np.logical_and(inside, xyz[:, 1] > self.origin[1])
+        inside = np.logical_and(inside, xyz[:, 1] < self.maximum[1])
+        inside = np.logical_and(inside, xyz[:, 2] > self.origin[2])
+        inside = np.logical_and(inside, xyz[:, 2] < self.maximum[2])
+        return inside
+
+    def regular_grid(self, nsteps=None, shuffle=False, order="C"):
+        if nsteps is None:
+            nsteps = self.nsteps
+        x = np.linspace(self.origin[0], self.maximum[0], nsteps[0])
+        y = np.linspace(self.origin[1], self.maximum[1], nsteps[1])
+        z = np.linspace(self.origin[2], self.maximum[2], nsteps[2])
+        xx, yy, zz = np.meshgrid(x, y, z, indexing="ij")
+        locs = np.array(
+            [xx.flatten(order=order), yy.flatten(order=order), zz.flatten(order=order)]
+        ).T
+        if shuffle:
+            # logger.info("Shuffling points")
+            np.random.shuffle(locs)
+        return locs

LoopStructural/datatypes/_bounding_box.py

@@ -0,0 +1,5 @@
+# from dataclasses import dataclass
+
+
+# @dataclass
+# class Value: