diff --git a/bgheatmaps/plane.py b/bgheatmaps/plane.py deleted file mode 100644 index 815a987..0000000 --- a/bgheatmaps/plane.py +++ /dev/null @@ -1,111 +0,0 @@ -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 \ No newline at end of file diff --git a/bgheatmaps/planner.py b/bgheatmaps/planner.py index d4033e0..92f8398 100644 --- a/bgheatmaps/planner.py +++ b/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, ) diff --git a/bgheatmaps/slicer.py b/bgheatmaps/slicer.py index 348cd2f..17b4ed0 100644 --- a/bgheatmaps/slicer.py +++ b/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]: