Merge branch 'gempy1.2-alpha' of https://github.com/CGR-Aachen/gempy …

…into gempy1.2-NewTheano

gempy/core/data.py

@@ -940,7 +940,7 @@ def read_surface_points(self, file_path, debug=False, inplace=False, append=Fals
         coord_x_name = kwargs.get('coord_x_name', "X")
         coord_y_name = kwargs.get('coord_y_name', "Y")
         coord_z_name = kwargs.get('coord_z_name', "Z")
-        surface_name = kwargs.get('surface_name', "formation")
+        surface_name = kwargs.get('surface_name', "surface")
         if 'sep' not in kwargs_pandas:
             kwargs_pandas['sep'] = ','
 
@@ -1280,7 +1280,7 @@ def read_orientations(self, filepath, debug=False, inplace=True, append=False, k
         azimuth_name = kwargs.get('azimuth_name', 'azimuth')
         dip_name = kwargs.get('dip_name', 'dip')
         polarity_name = kwargs.get('polarity_name', 'polarity')
-        surface_name = kwargs.get('surface_name', "formation")
+        surface_name = kwargs.get('surface_name', "surface")
 
         table = pn.read_csv(filepath, **kwargs_pandas)
 

gempy/core/theano/theano_graph.py

@@ -1360,23 +1360,51 @@ def compare(self, a, b, slice_init, Z_x, l, n_surface, drift):
             return T.le(Zx, a) * T.ge(Zx, b) * n_surface_0
 
     def select_finite_faults(self):
+        # get data points of fault
         fault_points = T.vertical_stack(T.stack(self.ref_layer_points[0]), self.rest_layer_points).T
-        ctr = T.mean(fault_points, axis=1)
-        x = fault_points - ctr.reshape((-1, 1))
-        M = T.dot(x, x.T)
-        U = T.nlinalg.svd(M)[2]
-        rotated_x = T.dot(self.x_to_interpolate(), U)
-        rotated_fault_points = T.dot(fault_points.T, U)
-        rotated_ctr = T.mean(rotated_fault_points, axis=0)
-        a_radio = (rotated_fault_points[:, 0].max() - rotated_fault_points[:, 0].min())/2 + self.inf_factor[self.n_surface_op[0]-1]
-        b_radio = (rotated_fault_points[:, 1].max() - rotated_fault_points[:, 1].min())/2 + self.inf_factor[self.n_surface_op[0]-1]
-        sel = T.lt((rotated_x[:, 0] - rotated_ctr[0])**2/a_radio**2 + (rotated_x[:, 1] - rotated_ctr[1])**2/b_radio**2,
-                   1)
+        # compute centroid of fault points
+        centroid = T.mean(fault_points, axis=1)
+        # compute difference of fault points from centroid
+        x = fault_points - centroid.reshape((-1, 1))
+        M = T.dot(x, x.T)  # same as np.cov(x) * 2
+        U = T.nlinalg.svd(M)  # is this the normal of the plane?
+        # overall this looks like some sort of plane fit to me
+        rotated_x = T.dot(self.x_to_interpolate(), U[0])  # this rotates ALL grid points that need to be interpolated
+        # rotated_x = T.dot(rotated_x, U[-1])  # rotate them with both rotation matrices
+        rotated_fault_points = T.dot(fault_points.T, U[0])  # same with fault points
+        # rotated_fault_points = T.dot(rotated_fault_points, U[-1])  # same
+        rotated_ctr = T.mean(rotated_fault_points, axis=0)  # and compute centroid of rotated points
+        # a factor: horizontal vector of ellipse of normal fault
+        a_radio = (rotated_fault_points[:, 0].max() - rotated_fault_points[:, 0].min()) / 2 \
+                  + self.inf_factor[self.n_surface_op[0] - 1]
+        # b_factor: vertical vector of ellipse
+        b_radio = (rotated_fault_points[:, 1].max() - rotated_fault_points[:, 1].min()) / 2 \
+                  + self.inf_factor[self.n_surface_op[0] - 1]
+
+        # sel = T.lt((rotated_x[:, 0] - rotated_ctr[0])**2 / a_radio**2 +
+        #            (rotated_x[:, 1] - rotated_ctr[1])**2 / b_radio**2,
+        #            1)
+
+        # ellipse equation: (x, c_x)^2 / a^2 +  (y - c_y)^2 / b^2 <= 1 if in ellipse
+        ellipse_factor = (rotated_x[:,0] - rotated_ctr[0])**2 / a_radio**2 + \
+            (rotated_x[:, 1] - rotated_ctr[1])**2 / b_radio**2
 
         if "select_finite_faults" in self.verbose:
-            sel = theano.printing.Print("scalar_field_iter")(sel)
+            ellipse_factor = theano.printing.Print("h")(ellipse_factor)
 
-        return sel
+        # h_factor = 1 - h
+        # if "select_finite_faults" in self.verbose:
+        #     h_factor = theano.printing.Print("h_factor")(h_factor)
+
+        # because we select all grid points as rotated_x, the selection here is
+        # a boolean for all grid points: True if in ellipse, False if outside ellipse
+
+        # if "select_finite_faults" in self.verbose:
+        #     sel = theano.printing.Print("scalar_field_iter")(sel)
+            # sum of boolean array sel is in my example: 38301
+            # so I guess this selects all grid points affected by this finite fault
+
+        return ellipse_factor  # sel
 
     def block_series(self, slope=5000, weights=None):
         """
@@ -1459,7 +1487,7 @@ def block_series(self, slope=5000, weights=None):
 
         return Z_x, partial_block
 
-    def block_fault(self, slope=50):
+    def block_fault(self, slope=50):  #
         """
         Compute the part of the block model of a given series (dictated by the bool array yet to be computed)
 
@@ -1483,32 +1511,47 @@ def block_fault(self, slope=50):
         # max_pot = theano.printing.Print("max_pot")(max_pot)
 
         min_pot = T.min(Z_x)
-        #     min_pot = theano.printing.Print("min_pot")(min_pot)
+        # min_pot = theano.printing.Print("min_pot")(min_pot)
 
         # max_pot_sigm = 2 * max_pot - self.scalar_field_at_surface_points_values[0]
         # min_pot_sigm = 2 * min_pot - self.scalar_field_at_surface_points_values[-1]
 
-        boundaty_pad = (max_pot - min_pot) * 0.01
+        boundary_pad = (max_pot - min_pot) * 0.01
         #l = slope / (max_pot - min_pot)  # (max_pot - min_pot)
 
-        # This is the different line with respect layers
-        l = T.switch(self.select_finite_faults(), 5000 / (max_pot - min_pot), 50 / (max_pot - min_pot))
-        #  l = theano.printing.Print("l")(l)
+        ellipse_factor = self.select_finite_faults()
+        ellipse_factor_rectified = T.switch(ellipse_factor < 1., ellipse_factor, 1.)
+
+        if "select_finite_faults" in self.verbose:
+            ellipse_factor_rectified = theano.printing.Print("h_factor_rectified")(ellipse_factor_rectified)
+
+        if "select_finite_faults" in self.verbose:
+            min_pot = theano.printing.Print("min_pot")(min_pot)
+            max_pot = theano.printing.Print("max_pot")(max_pot)
+
+        self.not_l = theano.shared(50.)
+        self.ellipse_factor_exponent = theano.shared(2)
+        # sigmoid_slope = (self.not_l * (1 / ellipse_factor_rectified)**3) / (max_pot - min_pot)
+        sigmoid_slope = 950 - 950 * ellipse_factor_rectified ** self.ellipse_factor_exponent + self.not_l
+        # l = T.switch(self.select_finite_faults(), 5000 / (max_pot - min_pot), 50 / (max_pot - min_pot))
+
+        if "select_finite_faults" in self.verbose:
+            sigmoid_slope = theano.printing.Print("l")(sigmoid_slope)
 
         # A tensor with the values to segment
         scalar_field_iter = T.concatenate((
-            T.stack([max_pot + boundaty_pad]),
+            T.stack([max_pot + boundary_pad]),
             self.scalar_field_at_surface_points_values,
-            T.stack([min_pot - boundaty_pad])
+            T.stack([min_pot - boundary_pad])
         ))
 
         if "scalar_field_iter" in self.verbose:
             scalar_field_iter = theano.printing.Print("scalar_field_iter")(scalar_field_iter)
 
         n_surface_op_float_sigmoid = T.repeat(self.n_surface_op_float[[0], :], 2, axis=1)
-
         # TODO: instead -1 at the border look for the average distance of the input!
-        n_surface_op_float_sigmoid = T.set_subtensor(n_surface_op_float_sigmoid[:, 0], -1)
+
+        n_surface_op_float_sigmoid = T.set_subtensor(n_surface_op_float_sigmoid[:, 1], -1)
         # - T.sqrt(T.square(n_surface_op_float_sigmoid[0] - n_surface_op_float_sigmoid[2])))
 
         n_surface_op_float_sigmoid = T.set_subtensor(n_surface_op_float_sigmoid[:, -1], -1)
@@ -1525,7 +1568,7 @@ def block_fault(self, slope=50):
             outputs_info=None,
             sequences=[dict(input=scalar_field_iter, taps=[0, 1]),
                        T.arange(0, n_surface_op_float_sigmoid.shape[1], 2, dtype='int64')],
-            non_sequences=[Z_x, l, n_surface_op_float_sigmoid, drift],
+            non_sequences=[Z_x, sigmoid_slope, n_surface_op_float_sigmoid, drift],
             name='Looping compare',
             profile=False,
             return_list=False)

gempy/plot/plot.py

@@ -32,7 +32,7 @@
 import numpy as _np
 import pandas as _pn
 import copy
-from .visualization import PlotData2D, steno3D, vtkVisualization
+from .visualization import PlotData2D, steno3D, vtkVisualization, ipyvolumeVisualization, cut_finite_fault_surfaces
 from .colors import cmap, norm, color_lot
 import gempy as _gempy
 
@@ -587,6 +587,11 @@ def plot_surfaces_3D(geo_data, vertices_l=None, simplices_l=None,
     return vv
 
 
+def plot_surfaces_3d_ipv(geo_model, ver, sim):
+    vv = ipyvolumeVisualization(geo_model, ver, sim)
+    vv.plot_ipyvolume()
+
+
 def export_to_vtk(geo_data, path=None, name=None, voxels=True, surfaces=True):
     """
       Export data to a vtk file for posterior visualizations

gempy/plot/visualization.py

@@ -38,6 +38,12 @@
 except ImportError:
     STENO_IMPORT = False
 
+try:
+    import ipyvolume as ipv
+    IPV_IMPORT = True
+except ImportError:
+    IPV_IMPORT = False
+
 import matplotlib.pyplot as plt
 import seaborn as sns
 import numpy as np
@@ -529,7 +535,7 @@ def annotate_plot(frame, label_col, x, y, **kwargs):
 
 class steno3D():
     def __init__(self, geo_data, project, **kwargs ):
-        if VTK_IMPORT is False:
+        if STENO_IMPORT is False:
             raise ImportError( 'Steno 3D package is not installed. No 3D online visualization available.')
         description = kwargs.get('description', 'Nothing')
 
@@ -1635,3 +1641,124 @@ def export_vtk_surfaces(vertices:dict, simplices, path=None, name='_surfaces', a
                 writer.SetInputData(polydata)
             writer.Write()
 
+
+class ipyvolumeVisualization:
+    def __init__(self, geo_model, ver, sim):
+        if VTK_IMPORT is False:
+            raise ImportError('ipyvolume package is not installed.')
+
+        self.geo_model = geo_model
+        self.ver = ver
+        self.sim = sim
+
+    def get_color_id(self, surface):
+        """Get id of given surface (str)."""
+        filter_ = self.geo_model.surfaces.df.surface == surface
+        color_id = self.geo_model.surfaces.df.id[filter_].values[0]
+        return color_id
+
+    def get_color(self, surface):
+        """Get color code of given gempy surface."""
+        return gp.plot.color_lot[self.get_color_id(surface)]
+
+    def plot_ipyvolume(self):
+        """Plot gempy surface model."""
+        ipv.figure()
+        meshes = []
+        for surf in self.ver.keys():
+            points = self.ver[surf]
+            triangles = self.sim[surf]
+            # color
+
+            mesh = ipv.plot_trisurf(points[:, 0] + self.geo_model.grid.extent[0],
+                                    points[:, 1] + self.geo_model.grid.extent[2],
+                                    points[:, 2] + self.geo_model.grid.extent[4],
+                                    triangles=triangles,
+                                    color=self.get_color(surf))
+            meshes.append(mesh)
+
+        ipv.xlim(self.geo_model.grid.extent[0], self.geo_model.grid.extent[1])
+        ipv.ylim(self.geo_model.grid.extent[2], self.geo_model.grid.extent[3])
+        ipv.zlim(self.geo_model.grid.extent[4], self.geo_model.grid.extent[5])
+        ipv.show()
+        return None
+
+
+def get_fault_ellipse_params(fault_points:np.ndarray):
+    """Get the fault ellipse parameters a and b from griven fault points (should
+    be the rotated ones.
+
+    Args:
+        fault_points (np.ndarray): Fault points
+
+    Returns:
+        (tuple) main axis scalars of fault ellipse a,b
+    """
+    a = (fault_points[:, 0].max() - fault_points[:, 0].min()) / 2
+    b = (fault_points[:, 1].max() - fault_points[:, 1].min()) / 2
+    return a, b
+
+
+def get_fault_rotation_objects(geo_model, fault:str):
+    """Gets fault rotation objects: rotation matrix U, the rotated fault points,
+    rotated centroid, and the ellipse parameters a and b.
+
+    Args:
+        geo_model (gempy.core.model.Model): gempy geo_model object
+        fault (str): Name of the fault surface.
+
+    Returns:
+        U (np.ndarray): Rotation matrix.
+        rfpts (np.ndarray): Rotated fault points.
+        rctr (np.array): Centroid of the rotated fault points.
+        a (float): Horizontal ellipse parameter.
+        b (float): Vertical ellipse parameter.
+    """
+    filter_ = geo_model.surface_points.df.surface == fault
+    fpts = geo_model.surface_points.df[filter_][["X", "Y", "Z"]].values.T
+    ctr = np.mean(fpts, axis=1)
+    x = fpts - ctr.reshape((-1, 1))
+    M = np.dot(x, x.T)
+    U = np.linalg.svd(M)
+    rfpts = np.dot(fpts.T, U[0])
+    # rfpts = np.dot(rfpts, U[-1])
+    rctr = np.mean(rfpts, axis=0)
+
+    a, b = get_fault_ellipse_params(rfpts)
+    return U, rfpts, rctr, a, b
+
+
+def cut_finite_fault_surfaces(geo_model, ver:dict, sim:dict):
+    """Cut vertices and simplices for finite fault surfaces to finite fault ellipse
+
+    Args:
+        geo_model (gempy.core.model.Model): gempy geo_model object
+        ver (dict): Dictionary with surfaces as keys and vertices ndarray as values.
+        sim (dict): Dictionary with surfaces as keys and simplices ndarray as values.
+
+    Returns:
+        ver, sim (dict, dict): Updated vertices and simplices with finite fault
+            surfaces cut to ellipses.
+    """
+    from scipy.spatial import Delaunay
+    from copy import copy
+
+    finite_ver = copy(ver)
+    finite_sim = copy(sim)
+
+    finite_fault_series = list(geo_model.faults.df[geo_model.faults.df["isFinite"] == True].index)
+    finite_fault_surfaces = list(
+        geo_model.surfaces.df[geo_model.surfaces.df.series == finite_fault_series].surface.unique())
+
+    for fault in finite_fault_surfaces:
+        U, fpoints_rot, fctr_rot, a, b = get_fault_rotation_objects(geo_model, "Fault 1")
+        rpoints = np.dot(ver[fault], U[0])
+        # rpoints = np.dot(rpoints, U[-1])
+        r = (rpoints[:, 0] - fctr_rot[0]) ** 2 / a ** 2 + (rpoints[:, 1] - fctr_rot[1]) ** 2 / b ** 2
+
+        finite_ver[fault] = finite_ver[fault][r < 1]
+        delaunay = Delaunay(finite_ver[fault])
+        finite_sim[fault] = delaunay.simplices
+        # finite_sim[fault] = finite_sim[fault][np.isin(sim[fault], np.argwhere(r<0.33))]
+
+    return finite_ver, finite_sim