Use infinite planes for slicing (#16)

bgheatmaps/plane.py

@@ -0,0 +1,111 @@
+from brainrender.actor import Actor
+from brainrender.scene import Scene
+from brainrender._utils import listify
+from typing import Dict, List, Self
+
+import numpy as np
+np.float = float # for compatibility with old vedo
+import vedo as vd
+
+try:
+    import vedo.vtkclasses as vtk
+except ImportError:
+    import vtkmodules.all as vtk
+
+vtk.vtkLogger.SetStderrVerbosity(vtk.vtkLogger.VERBOSITY_OFF)
+
+# from vedo 2023.4.6
+def intersect_with_plane(mesh: vd.Mesh, origin=(0, 0, 0), normal=(1, 0, 0)):
+    """
+    Intersect this Mesh with a plane to return a set of lines.
+
+    Example:
+        ```python
+        from vedo import *
+        sph = Sphere()
+        mi = sph.clone().intersect_with_plane().join()
+        print(mi.lines())
+        show(sph, mi, axes=1).close()
+        ```
+        ![](https://vedo.embl.es/images/feats/intersect_plane.png)
+    """
+    plane = vtk.vtkPlane()
+    plane.SetOrigin(origin)
+    plane.SetNormal(normal)
+
+    cutter = vtk.vtkPolyDataPlaneCutter()
+    cutter.SetInputData(mesh.polydata())
+    cutter.SetPlane(plane)
+    cutter.InterpolateAttributesOn()
+    cutter.ComputeNormalsOff()
+    cutter.Update()
+
+    msh = vd.Mesh(cutter.GetOutput(), "k", 1).lighting("off")
+    msh.GetProperty().SetLineWidth(3)
+    msh.name = "PlaneIntersection"
+    return msh
+
+class Plane:
+    def __init__(self, origin: np.ndarray, u: np.ndarray, v: np.ndarray) -> None:
+        self.center = origin
+        self.u = u / np.linalg.norm(u)
+        self.v = v / np.linalg.norm(v)
+        assert np.isclose(np.dot(self.u, self.v), 0), f"The plane vectors must be orthonormal to each other (u ⋅ v = {np.dot(self.u, self.v)})"
+        self.normal = np.cross(self.u, self.v)
+        self.M = np.vstack([u, v]).T
+
+    @staticmethod
+    def from_norm(origin: np.ndarray, norm: np.ndarray) -> Self:
+        u = np.zeros(3)
+        m = np.where(norm != 0)[0][0] # orientation can't be all-zeros
+        n = (m+1)%3
+        u[n] = norm[m]
+        u[m] = -norm[n]
+        norm = norm/np.linalg.norm(norm)
+        u = u/np.linalg.norm(u)
+        v = np.cross(norm, u)
+        return Plane(origin, u, v)
+
+    def to_mesh(self, actor: Actor):
+        bounds = actor.bounds()
+        length = max(
+            bounds[1] - bounds[0],
+            bounds[3] - bounds[2],
+            bounds[5] - bounds[4],
+        )
+        length += length / 3
+
+        plane_mesh = Actor(
+            vd.Plane(pos=self.center, normal=self.normal, sx=length, sy=length),
+            name=f"PlaneMesh at {self.center} norm: {self.normal}",
+            br_class="plane_mesh",
+        )
+        plane_mesh.width = length
+        return plane_mesh
+
+    def centerOfMass(self):
+        return self.center
+
+    def P3toP2(self, ps):
+        # ps is a list of 3D points
+        # returns a list of 2D point mapped on the plane (u -> x axis, v -> y axis)
+        return (ps - self.center) @ self.M
+
+    def intersectWith(self, mesh: vd.Mesh):
+        # mesh.intersect_with_plane(origin=self.center, normal=self.normal) in newer vedo
+        return intersect_with_plane(mesh, origin=self.center, normal=self.normal)
+
+    # for Slicer.get_structures_slice_coords()
+    def get_projections(self, actors: List[Actor]) -> Dict[str, np.ndarray]:
+        projected = {}
+        for actor in actors:
+            mesh: vd.Mesh = actor._mesh
+            intersection = self.intersectWith(mesh)
+            if not intersection.points().shape[0]:
+                continue
+            pieces = intersection.splitByConnectivity() # intersection.split() in newer vedo
+            for piece_n, piece in enumerate(pieces):
+                # sort coordinates
+                points = piece.join(reset=True).points()
+                projected[actor.name + f"_segment_{piece_n}"] = self.P3toP2(points)
+        return projected

bgheatmaps/planner.py

@@ -1,13 +1,14 @@
 from typing import Union
 import numpy as np
-from vedo import Arrow, Plane, Sphere
+from vedo import Arrow, Sphere
 from rich.table import Table
 from rich.panel import Panel
 from rich import print
 
-from myterial import pink_dark, blue_dark, amber_lighter, grey, amber, orange
+from myterial import pink_dark, blue_dark, green_dark, red_dark, amber_lighter, grey, amber, orange
 
 from bgheatmaps.heatmaps import heatmap
+from bgheatmaps.plane import Plane
 
 from brainrender import settings
 
@@ -23,21 +24,14 @@ def print_plane(name: str, plane: Plane, color: str):
     def fmt_array(x: np.ndarray) -> str:
         return str(tuple([round(v, 2) for v in x]))
 
-    # create a table to display the vertices posittion
-    vert_tb = Table(box=None)
-    vert_tb.add_column(style=f"{amber}", justify="right")
-    vert_tb.add_column(style="white")
-
-    for i in range(4):
-        vert_tb.add_row(f"({i})", fmt_array(plane.mesh.points()[i]))
-
     tb = Table(box=None)
     tb.add_column(style=f"bold {orange}", justify="right")
     tb.add_column(style="white")
 
-    tb.add_row("center point: ", fmt_array(plane.mesh.center))
-    tb.add_row("norm: ", str(tuple(plane.mesh.normal)))
-    tb.add_row("Vertices: ", vert_tb)
+    tb.add_row("center point: ", fmt_array(plane.center))
+    tb.add_row("norm: ", str(tuple(plane.normal)))
+    tb.add_row("u: ", str(tuple(plane.u)))
+    tb.add_row("v: ", str(tuple(plane.v)))
 
     print(
         Panel.fit(
@@ -91,21 +85,23 @@ def show(self):
             (blue_dark, pink_dark),
             (0.8, 0.3),
         ):
-            plane.alpha(alpha).color(color)
-
-            self.scene.add(plane, transform=False)
-            self.scene.add(
-                Arrow(
-                    plane.center,
-                    np.array(plane.center)
-                    + np.array(plane.mesh.normal) * self.arrow_scale,
-                    c=color,
-                ),
-                transform=False,
-            )
+            plane_mesh = plane.to_mesh(self.scene.root)
+            plane_mesh.alpha(alpha).color(color)
+
+            self.scene.add(plane_mesh, transform=False)
+            for vector, v_color in zip((plane.normal, plane.u, plane.v), (color, red_dark, green_dark)):
+                self.scene.add(
+                    Arrow(
+                        plane.center,
+                        np.array(plane.center)
+                        + np.array(vector) * self.arrow_scale,
+                        c=v_color,
+                    ),
+                    transform=False,
+                )
 
             self.scene.add(
-                Sphere(plane.center, r=plane.width / 125, c="k"),
+                Sphere(plane_mesh.center, r=plane_mesh.width / 125, c="k"),
                 transform=False,
             )
 

bgheatmaps/slicer.py

@@ -1,9 +1,9 @@
 from typing import List, Dict, Union, Optional
-from vedo import Plane
 import numpy as np
 
 from brainrender.actor import Actor
 from brainrender.scene import Scene
+from bgheatmaps.plane import Plane
 
 
 def get_ax_idx(orientation: str) -> int:
@@ -58,39 +58,39 @@ def __init__(
             p1 = position - shift
 
             # get the two planes
-            norm0, norm1 = np.zeros(3), np.zeros(3)
-            norm0[axidx] = 1
-            norm1[axidx] = -1
+            # assures that u0×v0 is all-positive -> it's for plane0
+            if orientation == "frontal":
+                u0, v0 = np.array([[0,0,-1], [0,1,0]])
+            elif orientation == "sagittal":
+                u0, v0 = np.array([[1,0,0], [0,1,0]])
+            else: # orientation == "horizontal"
+                u0, v0 = np.array([[0,0,-1], [-1,0,0]])
+            # u0, v0 = np.delete(np.array([[-1,0,0], [0,1,0], [0,0,-1]]), axidx, 0) # can't use this formula for backward compatibility
+            plane0 = Plane(position, u0, v0)
+            u1, v1 = u0.copy(), -v0.copy() # set u1:=u0 and v1:=-v0
+            # we could, alternatively set u1:=v0 and v1:=u0
+            # u1, v1 = v0.copy(), u0.copy()
+            plane1 = Plane(p1, u1, v1)
         else:
             orientation = np.array(orientation)
 
             p1 = position + orientation * thickness  # type: ignore
 
             norm0 = orientation  # type: ignore
+            plane0 = Plane.from_norm(position, norm0)
             norm1 = -orientation  # type: ignore
-
-        # get the length of the largest dimension of the atlas
-        bounds = root.bounds()
-        length = max(
-            bounds[1] - bounds[0],
-            bounds[3] - bounds[2],
-            bounds[5] - bounds[4],
-        )
-        length += length / 3
+            plane1 = Plane.from_norm(p1, norm1)
 
         self.plane0 = Actor(
-            Plane(pos=position, normal=norm0, sx=length, sy=length),
-            name=f"Plane at {position} norm: {norm0}",
+            plane0,
+            name=f"Plane at {plane0.center} norm: {plane0.normal}",
             br_class="plane",
         )
-        self.plane0.width = length
-
         self.plane1 = Actor(
-            Plane(pos=p1, normal=norm1, sx=length, sy=length),
-            name=f"Plane at {p1} norm: {norm1}",
+            plane1,
+            name=f"Plane at {plane1.center} norm: {plane1.normal}",
             br_class="plane",
         )
-        self.plane1.width = length
 
     def get_structures_slice_coords(self, regions: List[Actor], root: Actor):
         """
@@ -101,28 +101,7 @@ def get_structures_slice_coords(self, regions: List[Actor], root: Actor):
         """
         regions = regions + [root]
 
-        pts = self.plane0.points() - self.plane0.points()[0]
-        v = pts[1] / np.linalg.norm(pts[1])
-        w = pts[2] / np.linalg.norm(pts[2])
-
-        M = np.vstack([v, w]).T  # 3 x 2
-
-        projected: Dict[str, np.ndarray] = {}
-        for n, actor in enumerate(regions):
-            # get region/plane intersection
-            intersection = self.plane0.intersectWith(
-                actor._mesh.triangulate()
-            )  # points: (N x 3)
-
-            if not intersection.points().shape[0]:
-                continue  # no intersection
-
-            pieces = intersection.splitByConnectivity()
-            for piece_n, piece in enumerate(pieces):
-                # sort coordinates
-                points = piece.join(reset=True).points()
-
-                projected[actor.name + f"_segment_{piece_n}"] = points @ M
+        projected: Dict[str, np.ndarray] = self.plane0.get_projections(regions)
 
         # get output coordinates
         coordinates: Dict[str, List[np.ndarray]] = dict()
@@ -138,7 +117,7 @@ def show_plane_intersection(
         self, scene: Scene, regions: List[Actor], root: Actor
     ):
         """
-        Slices regions' meshjes with plane0 and adds the resulting intersection
+        Slices regions' meshes with plane0 and adds the resulting intersection
         to the brainrender scene.
         """
         for region in regions + [root]: