Plotting parametric objects as a grid in PyVista

Question

I am stuck with probably a simple problem but after reading pyvista docs I am still looking for an answer. I am trying to plot a grid in which each cell will be a mesh defined as a parametric shape i.e. supertorus. In an early version of pyvista, I defined "my own" supertorus as below:

def supertorus(yScale, xScale, Height, InternalRadius, Vertical, Horizontal,
           deltaX=0, deltaY=0, deltaZ=0):

#  initial range for values used in parametric equation
n = 100
u = np.linspace(-np.pi, np.pi, n)
t = np.linspace(-np.pi, np.pi, n)
u, t = np.meshgrid(u, t)

# a1: Y Scale <0, 2>
a1 = yScale
# a2: X Scale <0, 2>
a2 = xScale
# a3: Height <0, 5>
a3 = Height
# a4: Internal radius <0, 5>
a4 = InternalRadius
# e1: Vertical squareness <0.25, 1>
e1 = Vertical
# e2: Horizontal squareness <0.25, 1>
e2 = Horizontal

# Definition of parametric equation for supertorus
x = a1 * (a4 + np.sign(np.cos(u)) * np.abs(np.cos(u)) ** e1) *\
    np.sign(np.cos(t)) * np.abs(np.cos(t)) ** e2
y = a2 * (a4 + np.sign(np.cos(u)) * np.abs(np.cos(u)) ** e1) *\
    np.sign(np.sin(t)) * np.abs(np.sin(t)) ** e2
z = a3 * np.sign(np.sin(u)) * np.abs(np.sin(u)) ** e1

grid = pyvista.StructuredGrid(x + deltaX + 5, y + deltaY + 5, z + deltaZ)
return grid 

I could manipulate with deltaX, deltaY and deltaZ to position supertori at the location of my choice. Unfortunately, this approach was not efficient and I am planning to use PyVista provided supertoroidal meshes (https://docs.pyvista.org/examples/00-load/create-parametric-geometric-objects.html?highlight=supertoroid). My question is: how I can place multiple meshes (like supertori) at the location defined by coordinates x, y, z?

Answer 1

I believe what you're looking for are glyphs. You can pass your own dataset as a glyph geometry that will in turn plot the dataset in each point of the supermesh. Without going into details of orienting your glyphs, colouring them according to scalars and whatnot, here's a simple "alien invasion" scenario as an example:

import numpy as np
import pyvista as pv

# get dataset for the glyphs: supertoroid in xy plane
saucer = pv.ParametricSuperToroid(ringradius=0.5, n2=1.5, zradius=0.5)
saucer.rotate_y(90)
# saucer.plot()  #  <-- check how a single saucer looks like

# get dataset where to put glyphs
x,y,z = np.mgrid[-1:2, -1:2, :2]
mesh = pv.StructuredGrid(x, y, z)

# construct the glyphs on top of the mesh
glyphs = mesh.glyph(geom=saucer, factor=0.3)
# glyphs.plot()  #  <-- simplest way to plot it

# create Plotter and add our glyphs with some nontrivial lighting
plotter = pv.Plotter(window_size=(1000, 800))
plotter.add_mesh(glyphs, color=[0.2, 0.2, 0.2], specular=1, specular_power=15)

plotter.show()

I've added some strong specular lighting to make the saucers look more menacing:

But the key point for your problem was creating the glyphs from your supermesh by passing it as the geom keyword of mesh.glyph. The other keywords such as orient and scale are useful for arrow-like glyphs where you can use the glyph to denote vectorial information of your dataset.

---

You've asked in comments whether it's possible to vary the glyphs along the dataset. I was certain that this was not possible, however the VTK docs clearly mention the possibility to define a collection of glyphs to use:

More than one glyph may be used by creating a table of source objects, each defining a different glyph. If a table of glyphs is defined, then the table can be indexed into by using either scalar value or vector magnitude.

It turns out that PyVista doesn't expose this functionality (yet), but the base vtk package lets us get our hands dirty. Here's a proof of concept based on DataSetFilters.glyph, which I'll float by the PyVista devs to see if there's interest in exposing this functionality.

import numpy as np
import pyvista as pv
from pyvista.core.filters import _get_output  # just for this standalone example
import vtk
pyvista = pv  # just for this standalone example

# below: adapted from core/filters.py
def multiglyph(dataset, orient=True, scale=True, factor=1.0,
          tolerance=0.0, absolute=False, clamping=False, rng=None,
          geom_datasets=None, geom_values=None):
    """Copy a geometric representation (called a glyph) to every point in the input dataset.
    The glyphs may be oriented along the input vectors, and they may be scaled according to scalar
    data or vector magnitude.
    Parameters
    ----------
    orient : bool
        Use the active vectors array to orient the glyphs
    scale : bool
        Use the active scalars to scale the glyphs
    factor : float
        Scale factor applied to sclaing array
    tolerance : float, optional
        Specify tolerance in terms of fraction of bounding box length.
        Float value is between 0 and 1. Default is 0.0. If ``absolute``
        is ``True`` then the tolerance can be an absolute distance.
    absolute : bool, optional
        Control if ``tolerance`` is an absolute distance or a fraction.
    clamping: bool
        Turn on/off clamping of "scalar" values to range.
    rng: tuple(float), optional
        Set the range of values to be considered by the filter when scalars
        values are provided.
    geom_datasets : tuple(vtk.vtkDataSet), optional
        The geometries to use for the glyphs in table mode
    geom_values : tuple(float), optional
        The value to assign to each geometry dataset, optional
    """
    # Clean the points before glyphing
    small = pyvista.PolyData(dataset.points)
    small.point_arrays.update(dataset.point_arrays)
    dataset = small.clean(point_merging=True, merge_tol=tolerance,
                          lines_to_points=False, polys_to_lines=False,
                          strips_to_polys=False, inplace=False,
                          absolute=absolute)
    # Make glyphing geometry
    if not geom_datasets:
        arrow = vtk.vtkArrowSource()
        arrow.Update()
        geom_datasets = arrow.GetOutput(),
        geom_values = 0,
    # check if the geometry datasets are consistent
    if not len(geom_datasets) == len(geom_values):
        raise ValueError('geom_datasets and geom_values must have the same length!')
        # TODO: other kinds of sanitization, e.g. check for sequences etc.
    # Run the algorithm
    alg = vtk.vtkGlyph3D()
    if len(geom_values) == 1:
        # use a single glyph
        alg.SetSourceData(geom_datasets[0])
    else:
        alg.SetIndexModeToScalar()
        # TODO: index by vectors?
        # TODO: SetInputArrayToProcess for arbitrary arrays, maybe?
        alg.SetRange(min(geom_values), max(geom_values))
        # TODO: different Range?
        for val, geom in zip(geom_values, geom_datasets):
            alg.SetSourceData(val, geom)
    if isinstance(scale, str):
        dataset.active_scalars_name = scale
        scale = True
    if scale:
        if dataset.active_scalars is not None:
            if dataset.active_scalars.ndim > 1:
                alg.SetScaleModeToScaleByVector()
            else:
                alg.SetScaleModeToScaleByScalar()
    else:
        alg.SetScaleModeToDataScalingOff()
    if isinstance(orient, str):
        dataset.active_vectors_name = orient
        orient = True
    if rng is not None:
        alg.SetRange(rng)
    alg.SetOrient(orient)
    alg.SetInputData(dataset)
    alg.SetVectorModeToUseVector()
    alg.SetScaleFactor(factor)
    alg.SetClamping(clamping)
    alg.Update()
    return _get_output(alg)

def example():
    """Small glyph example"""

    rng = np.random.default_rng()

    # get dataset for the glyphs: supertoroid in xy plane
    # use N random kinds of toroids over a mesh with 27 points
    N = 5
    values = np.arange(N)  # values for scalars to look up glyphs by
    geoms = [pv.ParametricSuperToroid(n1=n1, n2=n2) for n1,n2 in rng.uniform(0.5, 2, size=(N, 2))]
    for geom in geoms:
        # make the disks horizontal for aesthetics
        geom.rotate_y(90)

    # get dataset where to put glyphs
    x,y,z = np.mgrid[-1:2, -1:2, -1:2]
    mesh = pv.StructuredGrid(x, y, z)

    # add random scalars
    mesh.point_arrays['scalars'] = rng.integers(0, N, size=x.size)

    # construct the glyphs on top of the mesh; don't scale by scalars now
    glyphs = multiglyph(mesh, geom_datasets=geoms, geom_values=values, scale=False, factor=0.3)

    # create Plotter and add our glyphs with some nontrivial lighting
    plotter = pv.Plotter(window_size=(1000, 800))
    plotter.add_mesh(glyphs, specular=1, specular_power=15)

    plotter.show()

if __name__ == "__main__":
    example()

The multiglyph function in the above is mostly the same as mesh.glyph, but I've replaced the geom keyword with two keywords, geom_datasets and geom_values. These define an index -> geometry mapping that is then used to look up each glyph based on array scalars.

You asked whether you can colour the glyphs independently: you can. In the above proof of concept the choice of glyph is tied to the scalars (choosing vectors would be equally easy; I'm not so sure about arbitrary arrays). However you can easily choose what arrays to colour by when you call pv.Plotter.add_mesh, so my suggestion is to use something other than the proper scalars to colour your glyphs.

Here's a typical output:

I kept the scalars for colouring to make it easier to see the differences between the glyphs. You can see that there are five different kinds of glyphs being chosen randomly based on the random scalars. If you set non-integer scalars it will still work; I suspect vtk chooses the closest scalar or something similar for lookup.