Revert "Use infinite planes for slicing"

bgheatmaps/planner.py

@@ -1,14 +1,13 @@
 from typing import Union
 import numpy as np
-from vedo import Arrow, Sphere
+from vedo import Arrow, Plane, Sphere
 from rich.table import Table
 from rich.panel import Panel
 from rich import print
 
-from myterial import pink_dark, blue_dark, green_dark, red_dark, amber_lighter, grey, amber, orange
+from myterial import pink_dark, blue_dark, amber_lighter, grey, amber, orange
 
 from bgheatmaps.heatmaps import heatmap
-from bgheatmaps.plane import Plane
 
 from brainrender import settings
 
@@ -24,14 +23,21 @@ 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.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)))
+    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)
 
     print(
         Panel.fit(
@@ -85,23 +91,21 @@ def show(self):
             (blue_dark, pink_dark),
             (0.8, 0.3),
         ):
-            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,
-                )
+            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,
+            )
 
             self.scene.add(
-                Sphere(plane_mesh.center, r=plane_mesh.width / 125, c="k"),
+                Sphere(plane.center, r=plane.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
-            # 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)
+            norm0, norm1 = np.zeros(3), np.zeros(3)
+            norm0[axidx] = 1
+            norm1[axidx] = -1
         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
-            plane1 = Plane.from_norm(p1, norm1)
+
+        # 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
 
         self.plane0 = Actor(
-            plane0,
-            name=f"Plane at {plane0.center} norm: {plane0.normal}",
+            Plane(pos=position, normal=norm0, sx=length, sy=length),
+            name=f"Plane at {position} norm: {norm0}",
             br_class="plane",
         )
+        self.plane0.width = length
+
         self.plane1 = Actor(
-            plane1,
-            name=f"Plane at {plane1.center} norm: {plane1.normal}",
+            Plane(pos=p1, normal=norm1, sx=length, sy=length),
+            name=f"Plane at {p1} norm: {norm1}",
             br_class="plane",
         )
+        self.plane1.width = length
 
     def get_structures_slice_coords(self, regions: List[Actor], root: Actor):
         """
@@ -101,7 +101,28 @@ def get_structures_slice_coords(self, regions: List[Actor], root: Actor):
         """
         regions = regions + [root]
 
-        projected: Dict[str, np.ndarray] = self.plane0.get_projections(regions)
+        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
 
         # get output coordinates
         coordinates: Dict[str, List[np.ndarray]] = dict()
@@ -117,7 +138,7 @@ def show_plane_intersection(
         self, scene: Scene, regions: List[Actor], root: Actor
     ):
         """
-        Slices regions' meshes with plane0 and adds the resulting intersection
+        Slices regions' meshjes with plane0 and adds the resulting intersection
         to the brainrender scene.
         """
         for region in regions + [root]: