Transform slice() output to 2D numpy array

[bcjiang]

I was trying to obtain a cross-section image from a 3D volume using the slice() method with normal vector input, but the output of slice() method is an object of <class 'pyvista.core.pointset.PolyData'>. I'm wondering if there is way to efficiently transform the output to a 2D numpy array data?

[banesullivan]

Check out the first post in this thread: #70

as long as the normal is oriented with the cartesian axial directions then it is a breeze. Otherwise, you'd have to rotate the slice

[banesullivan]

Also, if you want a 2D array, the mesh that you are slicing must be rectilinear (ie. a pyvista.UniformGrid or pyvista.RectilinearGrid) and the slice must be oriented with the axes, otherwise the slice will run into irregular geometries that cannot be represented in a 2D array

[bcjiang]

Thanks for the comments! I do have the pyvista.UniformGrid as the volume, but I would need slicing planes with arbitrary normals so it is a bit tricky.

[banesullivan]

What exactly is your goal with the 2D array? Because as soon as you tilt the slice, it will have non-rectangular geometries that can’t be well translated to a 2d grid (array). Maybe some insight will help me find a better solution for you?

_{Sent with GitHawk}

[bcjiang]

We have a 3D volume image data and then we would like to get some 2D slicing plane images out of it for some image analysis tasks. Something like this.

[banesullivan]

@bcjiang, #101 (comment) has a routine for resampling planar data to a rectilinear grid which will give you a 2D array. Thing is, you have to project the arbitrarily oriented plane to an axial plane.

A workaround would be to create a structured grid around your slice like so.

no matter what though, you have to resample the slice because of the irregular geometries. and during that resampling, you have to make a choice about resolution.

import pyvista as pv
from pyvista import examples
import numpy as np
import matplotlib.pyplot as plt

normal = np.array([1, 1, 0.5])
slc = examples.load_uniform().slice(normal=(1,1,0.5))
slc.plot(show_grid=True, show_edges=True)

# Make structured grid
x = np.linspace(0, 20)
y = np.linspace(0, 20)
z = np.array([0])
plane = pv.StructuredGrid(*np.meshgrid(x,y,z))

# rotate and translate grid to be ontop of the slice
vx = np.array([0., 0., 1.])
direction = normal / np.linalg.norm(normal)
vx -= vx.dot(direction) * direction
vx /= np.linalg.norm(vx)
vy = np.cross(direction, vx)
rmtx = np.array([vx, vy, direction])
plane.points = plane.points.dot(rmtx)
plane.points -= plane.center
plane.points += slc.center

# resample the data
sampled = plane.sample(slc, )
# Fill bad data
sampled["Spatial Cell Data"][~sampled["vtkValidPointMask"].astype(bool)] = np.nan

# plot the 2D array
array = sampled["Spatial Cell Data"].reshape(sampled.dimensions[0:2])

plt.pcolormesh(array)

[bcjiang]

Thanks! I just saw the comment, I will try that.

[banesullivan]

@bcjiang, I had to do this today to make some figures of slices using MPL and made a better version of this into a helper function:

def slice_to_array(slc, normal, origin, name, ni=50, nj=50):
    """Converts a PolyData slice to a 2D NumPy array.
    
    It is crucial to have the true normal and origin of
    the slicing plane
    
    Parameters
    ----------
    slc : PolyData
        The slice to convert.
    normal : tuple(float)
        the normal of the original slice
    origin : tuple(float)
        the origin of the original slice
    name : str
        The scalar array to fetch from the slice
    ni : int
        The resolution of the array in the i-direction
    nj : int
        The resolution of the array in the j-direction
    
    """
    # Make structured grid
    x = np.linspace(slc.bounds[0], slc.bounds[1], ni)
    y = np.linspace(slc.bounds[2], slc.bounds[3], nj)
    z = np.array([0])
    plane = pv.StructuredGrid(*np.meshgrid(x,y,z))

    # rotate and translate grid to be ontop of the slice
    vx = np.array([0., 0., 1.])
    direction = normal / np.linalg.norm(normal)
    vx -= vx.dot(direction) * direction
    vx /= np.linalg.norm(vx)
    vy = np.cross(direction, vx)
    rmtx = np.array([vx, vy, direction])
    plane.points = plane.points.dot(rmtx)
    plane.points -= plane.center
    plane.points += origin

    # resample the data
    sampled = plane.sample(slc, tolerance=slc.length*0.5)
    # Fill bad data
    sampled[name][~sampled["vtkValidPointMask"].view(bool)] = np.nan

    # plot the 2D array
    array = sampled[name].reshape(sampled.dimensions[0:2])
    return array

[bcjiang]

This is pretty helpful! Thanks a lot for the update.

[BobHeidsieck]

Hello,
I am using finite element method to manipulate anatomical structure. I have found Pyvista quite convenient to manage this structure. I would like to simulate an x-ray image of these structures. This is could be accomplished by integrating all slices along one axis and perform the 2D np array transformation.
Is the operation of integrating of the 3D structure along one axis to get a 2D structure a Pyvista existing method ?

Using several slices and performing the 2D np transform would be an sub optimal option due to multiple sampling.

Thanks for any advice.

[banesullivan]

@BobHeidsieck, please open a new issue.