plane direction

docs/changes.md

@@ -31,6 +31,7 @@
 - fixing axes type 10 by @Poisoned
 - improvements to input/output functionality for Assembly @ttsesm
 - added `mesh.remove_all_lines()` method
+- added keyword `Plane(edge_direction=...)` by @smoothumut
 
 
 
@@ -54,6 +55,8 @@ examples/pyplot/fit_curve2.py
 
 tests/issues/issue_1146.py
 tests/issues/discussion_1190.py
+tests/issues/test_sph_harm2.py
+
 tests/snippets/test_interactive_plotxy.py
 tests/snippets/test_elastic_pendulum.py
 ```

tests/issues/test_sph_harm2.py

@@ -0,0 +1,94 @@
+"""
+Morph one shape into another using spherical harmonics package shtools.
+
+In this example we morph a sphere into a octahedron and vice-versa.
+"""
+import numpy as np
+from vedo import settings, Plotter, Points, Sphere, cos, dataurl, mag, sin, Mesh
+
+try:
+    import pyshtools
+    print(__doc__)
+except ModuleNotFoundError:
+    print("Please install pyshtools to run this example")
+    print("Follow instructions at https://shtools.oca.eu/shtools")
+    exit(0)
+
+
+##########################################################
+N = 100  # number of sample points on the unit sphere
+lmax = 15  # maximum degree of the sph. harm. expansion
+rbias = 0.5  # subtract a constant average value
+x0 = [0, 0, 0]  # set object at this position
+##########################################################
+
+
+def makeGrid(shape, N):
+    rmax = 2.0  # line length
+    agrid, pts = [], []
+    for th in np.linspace(0, np.pi, N, endpoint=True):
+        lats = []
+        for ph in np.linspace(0, 2 * np.pi, N, endpoint=True):
+            p = np.array([sin(th) * cos(ph), sin(th) * sin(ph), cos(th)]) * rmax
+            intersections = shape.intersect_with_line([0, 0, 0], p)
+            if len(intersections):
+                value = mag(intersections[0])
+                lats.append(value - rbias)
+                pts.append(intersections[0])
+            else:
+                lats.append(rmax - rbias)
+                pts.append(p)
+        agrid.append(lats)
+    agrid = np.array(agrid)
+    actor = Points(pts, c="k", alpha=0.4, r=1)
+    return agrid, actor
+
+
+def morph(clm1, clm2, t, lmax):
+    """Interpolate linearly the two sets of sph harm. coeeficients."""
+    clm = (1 - t) * clm1 + t * clm2
+    grid_reco = clm.expand(lmax=lmax)  # cut "high frequency" components
+    agrid_reco = grid_reco.to_array()
+    pts = []
+    for i, longs in enumerate(agrid_reco):
+        ilat = grid_reco.lats()[i]
+        for j, value in enumerate(longs):
+            ilong = grid_reco.lons()[j]
+            th = np.deg2rad(90 - ilat)
+            ph = np.deg2rad(ilong)
+            r = value + rbias
+            p = np.array([sin(th) * cos(ph), sin(th) * sin(ph), cos(th)]) * r
+            pts.append(p)
+    return pts
+
+settings.use_depth_peeling = True
+
+plt = Plotter(shape=[2, 2], axes=3, interactive=0)
+
+shape1 = Sphere(alpha=0.2)
+shape2 = Mesh(dataurl + "icosahedron.vtk").normalize().linewidth(1)
+plt += shape2
+
+agrid1, actorpts1 = makeGrid(shape1, N)
+
+plt.at(0).show(shape1, actorpts1)
+
+agrid2, actorpts2 = makeGrid(shape2, N)
+plt.at(1).show(shape2, actorpts2)
+
+plt.camera.Zoom(1.2)
+
+clm1 = pyshtools.SHGrid.from_array(agrid1).expand()
+clm2 = pyshtools.SHGrid.from_array(agrid2).expand()
+# clm1.plot_spectrum2d()  # plot the value of the sph harm. coefficients
+# clm2.plot_spectrum2d()
+
+for t in np.arange(0, 1, 0.005):
+    act21 = Points(morph(clm2, clm1, t, lmax), c="r", r=4)
+    act12 = Points(morph(clm1, clm2, t, lmax), c="g", r=4)
+
+    plt.at(2).show(act21, resetcam=False)
+    plt.at(3).show(act12)
+    plt.azimuth(2)
+
+plt.interactive().close()

tests/snippets/test_interactive_plotxy1.py

@@ -0,0 +1,23 @@
+"""Interactive plot with a slider
+to change the k value of a sigmoid function."""
+import numpy as np
+from vedo import Plotter, settings
+from vedo.pyplot import plot
+
+kinit = 0.75
+n = 3
+x = np.linspace(-1, 1, 100)
+
+def update_plot(widget=None, event=""):
+    k = widget.value if widget else kinit
+    # y = 1/(1 + (k/x)**n) # hill function
+    # y = np.abs(k*x)**n *np.sign(k*x) # power law
+    y = (2 / (1 + np.exp(-np.abs(n*k*x))) - 1) *np.sign(k*x) # sigmoid
+    p = plot(x, y, c='red5', lw=4, xlim=(-1,1), ylim=(-1,1), aspect=1)
+    plt.remove("PlotXY").add(p)
+
+settings.default_font = "Roboto"
+plt = Plotter(size=(900, 1050))
+plt.add_slider(update_plot, -1, 1, value=kinit, title="k value", c="red3")
+update_plot()
+plt.show(__doc__, mode="2d", zoom='tight')

vedo/pointcloud.py

@@ -2106,7 +2106,13 @@ def warp(self, source, target, sigma=1.0, mode="3d") -> Self:
         self.pipeline = utils.OperationNode("warp", parents=parents)
         return self
 
-    def cut_with_plane(self, origin=(0, 0, 0), normal=(1, 0, 0), invert=False) -> Self:
+    def cut_with_plane(
+            self,
+            origin=(0, 0, 0),
+            normal=(1, 0, 0),
+            invert=False,
+            # generate_ids=False,
+    ) -> Self:
         """
         Cut the mesh with the plane defined by a point and a normal.
 
@@ -2115,6 +2121,8 @@ def cut_with_plane(self, origin=(0, 0, 0), normal=(1, 0, 0), invert=False) -> Se
                 the cutting plane goes through this point
             normal : (array)
                 normal of the cutting plane
+            invert : (bool)
+                select which side of the plane to keep
 
         Example:
             ```python
@@ -2153,13 +2161,29 @@ def cut_with_plane(self, origin=(0, 0, 0), normal=(1, 0, 0), invert=False) -> Se
         clipper.SetInputData(self.dataset)
         clipper.SetClipFunction(plane)
         clipper.GenerateClippedOutputOff()
-        clipper.GenerateClipScalarsOff()
+        clipper.SetGenerateClipScalars(0)
         clipper.SetInsideOut(invert)
         clipper.SetValue(0)
         clipper.Update()
 
-        self._update(clipper.GetOutput())
+        # if generate_ids:
+        #     saved_scalars = None # otherwise the scalars are lost
+        #     if self.dataset.GetPointData().GetScalars():
+        #         saved_scalars = self.dataset.GetPointData().GetScalars()
+        #     varr = clipper.GetOutput().GetPointData().GetScalars()
+        #     if varr.GetName() is None:
+        #         varr.SetName("DistanceToCut")
+        #     arr = utils.vtk2numpy(varr)
+        #     # array of original ids
+        #     ids = np.arange(arr.shape[0]).astype(int)
+        #     ids[arr == 0] = -1
+        #     ids_arr = utils.numpy2vtk(ids, dtype=int)
+        #     ids_arr.SetName("OriginalIds")
+        #     clipper.GetOutput().GetPointData().AddArray(ids_arr)
+        #     if saved_scalars:
+        #         clipper.GetOutput().GetPointData().AddArray(saved_scalars)
 
+        self._update(clipper.GetOutput())
         self.pipeline = utils.OperationNode("cut_with_plane", parents=[self])
         return self
 

vedo/shapes.py

@@ -2980,18 +2980,20 @@ def __init__(self, pos=(0, 0, 0), s=(1, 1), res=(10, 10), lw=1, c="k3", alpha=1.
 
 class Plane(Mesh):
     """Create a plane in space."""
-
     def __init__(
-            self,
-            pos=(0, 0, 0),
-            normal=(0, 0, 1),
-            s=(1, 1),
-            res=(1, 1),
-            c="gray5", alpha=1.0,
-        ) -> None:
+        self,
+        pos=(0, 0, 0),
+        normal=(0, 0, 1),
+        s=(1, 1),
+        res=(1, 1),
+        edge_direction=(),
+        c="gray5",
+        alpha=1.0,
+    ) -> None:
         """
         Create a plane of size `s=(xsize, ysize)` oriented perpendicular
-        to vector `normal` so that it passes through point `pos`.
+        to vector `normal` so that it passes through point `pos`, optionally
+        aligning an edge with `direction`.
 
         Arguments:
             pos : (list)
@@ -3002,36 +3004,58 @@ def __init__(
                 size of the plane along x and y
             res : (list)
                 resolution of the plane along x and y
+            edge_direction : (list)
+                direction vector to align one edge of the plane
         """
         if isinstance(pos, vtki.vtkPolyData):
             super().__init__(pos, c, alpha)
-            # self.transform = LinearTransform().translate(pos)
 
         else:
             ps = vtki.new("PlaneSource")
             ps.SetResolution(res[0], res[1])
             tri = vtki.new("TriangleFilter")
             tri.SetInputConnection(ps.GetOutputPort())
             tri.Update()
-
             super().__init__(tri.GetOutput(), c, alpha)
 
             pos = utils.make3d(pos)
             normal = np.asarray(normal, dtype=float)
             axis = normal / np.linalg.norm(normal)
+
+            # Calculate orientation using normal
             theta = np.arccos(axis[2])
             phi = np.arctan2(axis[1], axis[0])
 
             t = LinearTransform()
             t.scale([s[0], s[1], 1])
+
+            # Rotate to align normal
             t.rotate_y(np.rad2deg(theta))
             t.rotate_z(np.rad2deg(phi))
+
+            # Additional direction alignment
+            if len(edge_direction) >= 2:
+                direction = utils.make3d(edge_direction).astype(float)
+                direction /= np.linalg.norm(direction)
+
+                if s[0] <= s[1]:
+                    current_direction = np.asarray([0,1,0])
+                else:
+                    current_direction = np.asarray([1,0,0])
+
+                transformed_current_direction = t.transform_point(current_direction)
+                n = transformed_current_direction / np.linalg.norm(transformed_current_direction)
+
+                if np.linalg.norm(transformed_current_direction) >= 1e-6:
+                    angle = np.arccos(np.dot(n, direction))
+                    t.rotate(axis=axis, angle=np.rad2deg(angle))
+
             t.translate(pos)
             self.apply_transform(t)
 
         self.lighting("off")
         self.name = "Plane"
-        self.variance = 0
+        self.variance = 0 # used by pointcloud.fit_plane()
 
     def clone(self, deep=True) -> "Plane":
         newplane = Plane()

vedo/version.py

@@ -1 +1 @@
-_version = '2024.5.2+dev18'
+_version = '2024.5.2+dev19'