Image offset when using `Volume.slice_plane`

[olivierdelree]

I am working with a volume that I am cutting using Volume.slice_plane to display only the slice to the user. I can retrieve the slice just fine and can easily display it to the user. However, it seems the image is offset by some value compared to what I am expecting.

The coordinate framework I am working with has the X axis come towards the camera, the Y axis point downward, and the Z axis point to the right. Using controls for pitch and yaw, the user can rotate clockwise around the Z axis by pitch-degrees and around the new Y' axis by yaw-degrees. When only working with pitch, the centre of the image corresponds to the centre of the cutting plane. However, when yaw is introduced, the slice returned is offset from the centre. By this I mean that a plane cutting the volume in two halves with the same volume (plane parallel to the XY plane) does not cut the slice in two halves with the same area.

This is not much of a problem when simply displaying the slice but I am looking to show the user the volume coordinates of their cursor when hovering over the slice. Unfortunately, this offset makes it so I can't work backward from the pixel coordinates of the cursor since I don't know where on the image the plane centre with which I cut the volume is.

My question is: is it possible to figure out what this offset is? If so, how?
Additionally, this might be a naive approach and there is a better way of doing this I've not thought about. If so, I'd love to hear about it. (keep in mind this would need to be computed many times per second, preferably 30, so it would need to be fast)

Here is as minimal an example I could make:

import matplotlib.pyplot as plt
import numpy as np
from scipy.spatial.transform import Rotation
import vedo

def compute_normal(pitch: int, yaw: int) -> np.ndarray:
    normal = [-1, 0, 0]
    rotation = Rotation.from_euler("ZY", [pitch, yaw], degrees=True)

    return rotation.apply(normal)

def slice_image(volume: vedo.Volume, offset: int, pitch: int, yaw: int) -> np.ndarray:
    slice_mesh = volume.slice_plane(
        origin=volume.center() + [-offset, 0, 0],
        normal=compute_normal(pitch, yaw),
    )

    # I would also normally rotate the image by pitch so that the intersection line
    # points up on the image.
    return slice_mesh.pointdata["ImageScalars"].reshape(slice_mesh.metadata["shape"])

offset = 50
pitch = 10
yaw = 40

array = np.zeros(shape=(200, 150, 180), dtype=np.uint8)
array[..., (array.shape[2] - 1) // 2] = 1

volume = vedo.Volume(array)

slice_plane = vedo.Plane(
    pos=volume.center() + [-offset, 0, 0],
    normal=compute_normal(pitch, yaw),
    s=(array.shape[2], array.shape[1]),
)

camera = dict(
    position=volume.center() + [500, 0, 0],
    focal_point=volume.center(),
    viewup=(0, -1, 0),
)

plotter = vedo.Plotter(axes=3)
plotter += [volume, slice_plane]
plotter.show(camera=camera)

image = slice_image(volume, offset, pitch, yaw)

# Show the center of the image
image[(image.shape[0] - 1) // 2] = 2
image[:, (image.shape[1] - 1) // 2] = 2

plt.imshow(image)
plt.show()

[marcomusy]

Hi you also have slice_mesh.metadata["original_bounds"] which gives you what it says and then compute the center of it.. or just slice_mesh.bounds().

Or maybe you like this alternative solutuion?

from vedo import *

normal = [0, 0, 1]
cmap = "gist_stern_r"

def func(w, _):
    c, n = pcutter.origin, pcutter.normal
    vslice = vol.slice_plane(c, n, autocrop=True).cmap(cmap)
    vslice.name = "Slice"
    plt.at(1).remove("Slice").add(vslice)

def keypress(event):
    if event.keypress == "q":
        plt.close()
        return
    if event.keypress == "s":
        print("Plane origin:", pcutter.origin, "normal:", pcutter.normal)

vol = Volume(dataurl + "embryo.slc").cmap(cmap)
vslice = vol.slice_plane(vol.center(), normal).cmap(cmap)
vslice.name = "Slice"

settings.enable_default_keyboard_callbacks = False

plt = Plotter(axes=0, N=2, bg="k", bg2="bb", interactive=False)
plt.add_callback("keypress", keypress)

pcutter = PlaneCutter(
    vslice,
    normal=normal,
    alpha=0,
    c="white",
    padding=0,
    can_translate=False,
    can_scale=False,
)
pcutter.add_observer("interaction", func)
plt.at(0).show(vol, zoom=1.5)
plt.at(1).add(pcutter)
plt.interactive()
plt.close()