Why does removing a transform fix a lot of our issues?

[alessandrofelder]

We've closed #255 by removing a transform - which fixes lots of other issues 🎉 but before we merge and release this we'd like to know and document why we do/don't need transformations.

This issue can be closed when we have at least full documentation (and ideally also know what change in brainrender, vedo or vtk made the transform redundant).

https://github.com/marcomusy/vedo/releases/tag/v2022.3.0 is a suspect, but @IgorTatarnikov and I couldn't understand why any of the changes between this and the previous release in vedo/volume.py would affect our code. It seems there are minor changes around how BaseVolume::getDataArray (now BaseVolume::tonumpy) handles the axis ordering, but these are not sufficient to explain why we don't need this z-axis flip, as we expect this function to behave the same way it did before despite the changes. 🤔

We ran git diff b340375 f0254ff vedo/volume.py (commits refer to vedo v2022.3.0 and the previous release).

[alessandrofelder]

To document this slightly better, part of the output of the git diff command above is

##########################################################################
 class BaseVolume:
     """Base class. Do not instantiate."""
-    def __init__(self, inputobj=None):
+
+    def __init__(self):
         self._data = None
         self._mapper = None

@@ -93,8 +97,19 @@ class BaseVolume:
         when all your modifications are completed.
         """
         narray_shape = tuple(reversed(self._data.GetDimensions()))
-        narray = utils.vtk2numpy(self._data.GetPointData().GetScalars()).reshape(narray_shape)
+
+        scals = self._data.GetPointData().GetScalars()
+        comps = scals.GetNumberOfComponents()
+        if comps == 1:
+            narray = utils.vtk2numpy(scals).reshape(narray_shape)
             narray = np.transpose(narray, axes=[2, 1, 0])
+        else:
+            narray = utils.vtk2numpy(scals).reshape(*narray_shape, comps)
+            narray = np.transpose(narray, axes=[2, 1, 0, 3])
+
+        # narray = utils.vtk2numpy(self._data.GetPointData().GetScalars()).reshape(narray_shape)
+        # narray = np.transpose(narray, axes=[2, 1, 0])
+
         return narray

which doesn't explain why there would be a change in the axis ordering?

May be clearer without the --ignore-all-whitespace option:

-    @deprecated(reason=colors.red+"Please use tonumpy()"+colors.reset)
+    @deprecated(reason=colors.red + "Please use tonumpy()" + colors.reset)
     def getDataArray(self):
         """Deprecated. Please use tonumpy()"""
         return self.tonumpy()
@@ -93,8 +97,19 @@ class BaseVolume:
         when all your modifications are completed.
         """
         narray_shape = tuple(reversed(self._data.GetDimensions()))
-        narray = utils.vtk2numpy(self._data.GetPointData().GetScalars()).reshape(narray_shape)
-        narray = np.transpose(narray, axes=[2, 1, 0])
+
+        scals = self._data.GetPointData().GetScalars()
+        comps = scals.GetNumberOfComponents()
+        if comps == 1:
+            narray = utils.vtk2numpy(scals).reshape(narray_shape)
+            narray = np.transpose(narray, axes=[2, 1, 0])
+        else:
+            narray = utils.vtk2numpy(scals).reshape(*narray_shape, comps)
+            narray = np.transpose(narray, axes=[2, 1, 0, 3])
+
+        # narray = utils.vtk2numpy(self._data.GetPointData().GetScalars()).reshape(narray_shape)
+        # narray = np.transpose(narray, axes=[2, 1, 0])
+
         return narray

[alessandrofelder]

Before and after have the same call to np.transpose for the first three axes...

[IgorTatarnikov]

I got some more insight into this today. It seems that the apply_transform function is acting differently depending on the mesh. When I run:

from brainrender import Scene
from brainrender.actors import Neuron, Cylinder
import numpy as np

angle = np.pi / 2
rot_matrix = [[np.cos(angle), 0, np.sin(angle)], [0, 1, 0], [-np.sin(angle), 0, np.cos(angle)]]

# Create a brainrender scene
scene = Scene(title="neurons")

# Rotate the root mesh
root_copy = scene.actors[0]._mesh.clone()
scene.add(root_copy.apply_transform(rot_matrix))

# Add a neuron from file
neuron = Neuron("examples/data/neuron1.swc")
scene.add(neuron)

# Rotate the neuron
neuron_copy = neuron._mesh.clone()
scene.add(neuron_copy.apply_transform(rot_matrix))

# Add a brain region
th = scene.add_brain_region(
    "TH",
    alpha=0.4,
)
scene.add(th)

# Rotate the brain region
th_copy = th._mesh.clone()
scene.add(th_copy.apply_transform(rot_matrix))

# create and add a cylinder actor
actor = Cylinder(
    th,  # center the cylinder at the center of mass of th
    scene.root,  # the cylinder actor needs information about the root mesh
)
scene.add(actor)
# Rotate the cylinder
actor_copy = actor._mesh.clone()
scene.add(actor_copy.apply_transform(rot_matrix))

# Render!
scene.render(axes=2)

I get:

Everything rotates correctly except for the Cylinder object, for some reason it snaps to the origin. Changing the rotation matrix to be [[1,0,0,0], [0,1,0,0], [0,0,1,0], [0,0,0,1]] has everything stay in place and overlap except the cylinder:

I think this was causing some elements to end up out of position when this line runs.

[marcomusy]

Hi with these fixes in brainrender mega-fix branch most mesh-related examples work with minor issues on vedo 2023.5.0+dev26a as in #246 :

diff --git a/brainrender/actor.py b/brainrender/actor.py
index b3f9cf8..c8bef43 100644
--- a/brainrender/actor.py
+++ b/brainrender/actor.py
@@ -45,7 +45,7 @@ def make_actor_label(
             color = [0.2, 0.2, 0.2]

         # Get mesh's highest point
-        points = actor.mesh.points().copy()
+        points = actor.mesh.vertices.copy()
         point = points[np.argmin(points[:, 1]), :]
         point += np.array(offset) + default_offset
         point[2] = -point[2]
@@ -64,12 +64,12 @@ def make_actor_label(

         # Mark a point on Mesh that corresponds to the label location
         if radius is not None:
-            pt = actor.closestPoint(point)
+            pt = actor.closest_point(point)
             pt[2] = -pt[2]
             sphere = Sphere(pt, r=radius, c=color, res=8)
             sphere.ancor = pt
             new_actors.append(sphere)
-            sphere.computeNormals()
+            sphere.compute_normals()

     return new_actors

@@ -159,7 +159,7 @@ class Actor(object):
         """
         Returns the coordinates of the mesh's center
         """
-        return self.mesh.points().mean(axis=0)
+        return self.mesh.center_of_mass()

     @classmethod
     def make_actor(cls, mesh, name, br_class):
@@ -218,7 +218,7 @@ class Actor(object):
             f"[{orange}]center of mass:[/{orange}][{amber}] {self.mesh.center_of_mass().astype(np.int32)}"
         )
         rep.add(
-            f"[{orange}]number of vertices:[/{orange}][{amber}] {len(self.mesh.points())}"
+            f"[{orange}]number of vertices:[/{orange}][{amber}] {self.mesh.npoints}"
         )
         rep.add(
             f"[{orange}]dimensions:[/{orange}][{amber}] {np.array(self.mesh.bounds()).astype(np.int32)}"
diff --git a/brainrender/actors/cylinder.py b/brainrender/actors/cylinder.py
index 715b268..0c1d222 100644
--- a/brainrender/actors/cylinder.py
+++ b/brainrender/actors/cylinder.py
@@ -20,7 +20,7 @@ class Cylinder(Actor):

         # Get pos
         if isinstance(pos, Mesh):
-            pos = pos.points().mean(axis=0)
+            pos = pos.center_of_mass()
         elif isinstance(pos, Actor):
             pos = pos.center
         logger.debug(f"Creating Cylinder actor at: {pos}")
diff --git a/brainrender/actors/points.py b/brainrender/actors/points.py
index 5568c09..22cbd08 100644
--- a/brainrender/actors/points.py
+++ b/brainrender/actors/points.py
@@ -140,7 +140,7 @@ class PointsDensity(Actor):
         volume = (
             vPoints(data)
             .density(dims=dims, radius=radius, **kwargs)
-            .c("Dark2")
+            .cmap("Dark2")
             .alpha([0, 0.9])
             .mode(1)
         )  # returns a vedo Volume
diff --git a/brainrender/actors/ruler.py b/brainrender/actors/ruler.py
index f4a98f9..735cec5 100644
--- a/brainrender/actors/ruler.py
+++ b/brainrender/actors/ruler.py
@@ -37,7 +37,7 @@ def ruler(p1, p2, unit_scale=1, units=None, s=50):
     dist = mag(p2 - p1) * unit_scale
     label = precision(dist, 3) + " " + units
     lbl = Text3D(label, pos=midpoint, s=s + 100, justify="center")
-    lbl.SetOrientation([0, 0, 180])
+    # lbl.SetOrientation([0, 0, 180])
     actors.append(lbl)

     # Add spheres add end
diff --git a/examples/cell_density.py b/examples/cell_density.py
index fcfa2bf..2a97fa7 100644
--- a/examples/cell_density.py
+++ b/examples/cell_density.py
@@ -27,7 +27,7 @@ def get_n_random_points_in_region(region, N):
     Z = np.random.randint(region_bounds[4], region_bounds[5], size=10000)
     pts = [[x, y, z] for x, y, z in zip(X, Y, Z)]

-    ipts = region.mesh.inside_points(pts).points()
+    ipts = region.mesh.inside_points(pts).vertices
     return np.vstack(random.choices(ipts, k=N))

diff --git a/examples/screenshot.py b/examples/screenshot.py
index 6fd7e46..efe6b50 100644
--- a/examples/screenshot.py
+++ b/examples/screenshot.py
@@ -32,7 +32,7 @@ camera = {
     "focalPoint": (7718, 4290, -3507),
     "distance": 40610,
 }
-zoom = 1.5
+zoom = 2.5

 # If you only want a screenshot and don't want to move the camera
 # around the scene, set interactive to False.
diff --git a/examples/user_volumetric_data.py b/examples/user_volumetric_data.py
index 3dd1c32..ddfd7f8 100644
--- a/examples/user_volumetric_data.py
+++ b/examples/user_volumetric_data.py
@@ -74,7 +74,7 @@ transformed_stack = source_space.map_stack_to(target_space, data)

 # 3. create a Volume vedo actor and smooth
 print("Creating volume")
-vol = Volume(transformed_stack, origin=scene.root.origin()).smooth_median()
+vol = Volume(transformed_stack).smooth_median()

I suggest to try out this vedo version as it will be released soon, and leaving aside the alignment issues in the volumetric examples just for the moment as they may need some extra work.

[IgorTatarnikov]

Sounds good. Thanks for your input! I've added these changes to the compatibility-with-vedo-2023.5.0+dev26a branch.

[marcomusy]

Nice! Adding volume.permute_axes() also fixes the volume examples in compatibility-with-vedo-2023.5.0+dev26a branch:

diff --git a/benchmark/bm_cells.py b/benchmark/bm_cells.py
index 6170d8e..4af2f5f 100644
--- a/benchmark/bm_cells.py
+++ b/benchmark/bm_cells.py
@@ -19,7 +19,7 @@ def get_n_random_points_in_region(region, N):
     Z = np.random.randint(region_bounds[4], region_bounds[5], size=10000)
     pts = [[x, y, z] for x, y, z in zip(X, Y, Z)]

-    ipts = region.mesh.insidePoints(pts).points()
+    ipts = region.mesh.inside_points(pts).coordinates
     return np.vstack(random.choices(ipts, k=N))

diff --git a/benchmark/bm_napari.py b/benchmark/bm_napari.py
index d802d45..7968ed8 100644
--- a/benchmark/bm_napari.py
+++ b/benchmark/bm_napari.py
@@ -21,7 +21,7 @@ for reg in regions[:400]:
 surfaces = []
 for act in scene.clean_actors:
     surfaces.append(
-        (act.points(), act.faces(), np.ones(len(act.points())) * 0.5)
+        (act.vertices, act.cells, np.ones(len(act.vertices)) * 0.5)
     )

 # render stuff in napar
diff --git a/benchmark/timer.py b/benchmark/timer.py
index e6eb560..2422a9d 100644
--- a/benchmark/timer.py
+++ b/benchmark/timer.py
@@ -37,7 +37,7 @@ class SceneContent:
         )
         for act in self.scene.clean_actors:
             try:
-                points = len(act.points())
+                points = act.npoints
             except AttributeError:
                 points = 0
             actors.add(
@@ -71,7 +71,7 @@ class Timer:
         points = []
         for act in self.scene.clean_actors:
             try:
-                points.append(len(act.points()))
+                points.append(act.npoints)
             except AttributeError:
                 pass
         tot_points = sum(points)
diff --git a/examples/add_cells.py b/examples/add_cells.py
index 42dcfb3..a52b0c5 100644
--- a/examples/add_cells.py
+++ b/examples/add_cells.py
@@ -22,7 +22,7 @@ def get_n_random_points_in_region(region, N):
     Z = np.random.randint(region_bounds[4], region_bounds[5], size=10000)
     pts = [[x, y, z] for x, y, z in zip(X, Y, Z)]

-    ipts = region.mesh.inside_points(pts).points()
+    ipts = region.mesh.inside_points(pts).coordinates
     return np.vstack(random.choices(ipts, k=N))

diff --git a/examples/user_volumetric_data.py b/examples/user_volumetric_data.py
index ddfd7f8..19b18e1 100644
--- a/examples/user_volumetric_data.py
+++ b/examples/user_volumetric_data.py
@@ -74,7 +74,8 @@ transformed_stack = source_space.map_stack_to(target_space, data)

 # 3. create a Volume vedo actor and smooth
 print("Creating volume")
-vol = Volume(transformed_stack).smooth_median()
+vol = Volume(transformed_stack).permute_axes(2, 1, 0)
+vol.smooth_median()

 # 4. Extract a surface mesh from the volume actor
diff --git a/tests/test_points.py b/tests/test_points.py
index 5191566..d6a6800 100644
--- a/tests/test_points.py
+++ b/tests/test_points.py
@@ -19,7 +19,7 @@ def get_n_random_points_in_region(region, N):
     Z = np.random.randint(region_bounds[4], region_bounds[5], size=10000)
     pts = [[x, y, z] for x, y, z in zip(X, Y, Z)]

-    ipts = region.mesh.inside_points(pts).points()
+    ipts = region.mesh.inside_points(pts).coordinates
     return np.vstack(random.choices(ipts, k=N))

[alessandrofelder]

Closing as we are confident transforms in vedo and brainrender are now working as we expect (even though I personally don't understand exactly why). Please open a new issue for any newly appearing transform problems!