3D deformation grids visualisation.

[antonfrancois]

Hello !

I am trying to implement a tool or API to visualize deformation grids in 3D. In 2D the visualization is straight forward. You can see some examples in one of my repository jupyters. I apologize in advance for the blurriness of my question but I am a bit lost here.

After battling with matplotlib, I got the advice from my supervisors to try other libraries as the 3D of matplotlib isn't really 3D. I spotted vedo, was amazed by the beautiful plots and tried to wrap my head around it.
I really liked the one in the thinplate grid example. But I can't manage to plot the deformation grids I generate on my own. Actually I'm a bit lost on how to make basic things.

I would like to plot different transparent surfaces as in the thinplate example, but I cannot find how to do it from a given numpy array. Adding some steamlines may also add some understanding of the deformation.

I am working with deformations that are functions in the PyTorch convention, and can be generated like so :

# Dimension definition
Tau,D,H,W = (1,100,150,200) # we will use a temporal vector field at a given time
x,y,z = np.mgrid[:D,:H,:W]
grid = np.stack([x,y,z],axis=-1)[None]
v = np.zeros((Tau,D,H,W,3)) # vector field initialization 
coef_mat = np.exp(-5*((grid[0,...,0])**2 +
                         (grid[0,...,1])**2 +
                         (grid[0,...,2])**2))[...,None]
v[0] = coef_mat*np.stack(
    (grid[0,...,1],
     -grid[0,...,2],
     np.ones(grid[0,...,0].shape)),
axis=3)
deformation = (grid + v[0][None])

Thank you for reading !
Anton

[marcomusy]

Hi Anton,
Thanks for your nice words! You can create a Mesh by defining vertices and faces with Mesh([pts,faces]). There are classes that simplify this for simple objects ike "Grid". E.g.:

from vedo import *
import numpy as np

xcoords = np.arange(0, 2, 0.2)
ycoords = np.arange(0, 1, 0.2)
sqrtx = sqrt(xcoords) # to make them non-uniform

grid1 = Grid(sx=sqrtx, sy=ycoords).lineWidth(3)
print(grid1.points())
grid2 = Grid(sx=sqrtx, sy=ycoords).lineWidth(3)
grid2.rotateY(15).shift(0,0.5,0.8)

arrows = Arrows(grid1, grid2, s=0.3).lighting('off')

grid = merge(grid1, grid2) # optionally fuse them into a single object

plt = show(grid, arrows, axes=1)  # returns a Plotter object

For deformations you need to specify where specific points in space displace to interpolate all space displacement. It does not need to be a point of the Mesh.
If you already know the displacements for all points then it's trivial - you just need to set the new coordinates with something like mymesh.points(my_new_points).

Hope this helps.

[antonfrancois]

Hello and thank you for the fast answer.

So yes I know the displacements for all points (I deal with custom non linear deformation).
I'm not sure if the Grid object can help me in my case. I tied to initialize a regular grid and then changing it's points with only getting messy results.
I then tried to use the Mesh object which gets me closer but I have only a bunch of unconnected points. I didn't understood how to set the connectivity matrix, and is there a simple way to do it ?.

I'm getting the updated positions of a grid by applying a given vector field v that I compute, and it look like :

from vedo import *
import numpy as np
#
# Generating a given 3D deformation
Tau,D,H,W = (1,10,15,20) # we will use a temporal vector field at a given time
x,y,z = np.meshgrid(np.linspace(-1,1,H),
                    np.linspace(-1,1,D),
                    np.linspace(-1,1,W))
grid = np.stack([x,y,z],axis=-1)[None]
v = np.zeros((Tau,D,H,W,3)) # vector field initialization
coef_mat = np.exp(-5*((grid[0,...,0])**2 +
                         (grid[0,...,1])**2 +
                         (grid[0,...,2])**2))[...,None]
v[0] = coef_mat*np.stack(
    (grid[0,...,1],
     -grid[0,...,2],
     np.ones(grid[0,...,0].shape)),
axis=3)
deformation = (grid + v[0][None]) # deformation shape is [1,D,H,W,3]

# Plots
deformed_grid = Mesh(deformation[0,5].reshape((H*W,3))).lw(3)
pts = Points(deformation[0,5].reshape((H*W,3)),c='red')

plt = show(deformed_grid,pts, __doc__, bg='k', bg2='lg', axes=11)  # returns a Plotter object

[marcomusy]

You need to find the mapping between your meshgrid numpy array and the Mesh object.

from vedo import *
import numpy as np

# Generating a given 3D deformation
Tau,D,H,W = (1,10,15,20)  # we will use a temporal vector field at a given time
x,y,z = np.meshgrid(np.linspace(-1,1,H),
                    np.linspace(-1,1,D),
                    np.linspace(-1,1,W))
grid = np.stack([x,y,z],axis=-1)[None]
v = np.zeros((Tau,D,H,W,3)) # vector field initialization
coef_mat = np.exp(-5*((grid[0,...,0])**2 +
                         (grid[0,...,1])**2 +
                         (grid[0,...,2])**2))[...,None]
v[0] = coef_mat*np.stack(
     (grid[0,...,1],
     -grid[0,...,2],
     np.ones(grid[0,...,0].shape)), axis=3)

deformation = (grid + v[0][None]) # deformation shape is [1,D,H,W,3]

# Plots
points = deformation[0,5].reshape((H*W,3))
faces = [[0,1,21,20], [1,2,22,21], [21,22,42,41], [83,84,104,103]]
deformed_grid = Mesh([points, faces]).lineWidth(1)

show(deformed_grid, 
     Points(points, c='red6', r=5),
     deformed_grid.labels('id', rotX=-45, scale=.02).c('yellow'),
     deformed_grid.labels('id', rotX=-90, scale=.03, justify='center', cells=1).shift(dy=.005),
     bg='k', bg2='lg', axes=9,
) 

this loop might be useful (or might try to do something more clever in numpy):

vedo/vedo/shapes.py

Line 2216 in 263fbf7

faces = []

PS: if you need triangular faces just add deformed_grid.triangulate():

[antonfrancois]

Thanks a lot for the help ! I understand better how the library works now. I added here one loop with numpy to mark the faces, here is a working example :

from vedo import *
import numpy as np

# Generating a given 3D deformation
Tau,D,H,W = (1,10,15,20)  # we will use a temporal vector field at a given time
x,y,z = np.meshgrid(np.linspace(-1,1,H),
                    np.linspace(-1,1,D),
                    np.linspace(-1,1,W))
grid = np.stack([x,y,z],axis=-1)[None]
v = np.zeros((Tau,D,H,W,3)) # vector field initialization
coef_mat = np.exp(-5*((grid[0,...,0])**2 +
                         (grid[0,...,1])**2 +
                         (grid[0,...,2])**2))[...,None]
v[0] = coef_mat*np.stack(
     (grid[0,...,1],
     -grid[0,...,2],
     np.ones(grid[0,...,0].shape)), axis=3)

deformation = (grid + v[0][None]) # deformation shape is [1,D,H,W,3]

# Plots
points = deformation[0,3].reshape((H*W,3))

# construct list containing face indexes
coord = np.arange(H*W).reshape((H,W))
faces = np.zeros((H*W,4))
for h in range(H-1):
    faces[h*W:h*W + (W-1),:] = np.stack([coord[h,:-1],coord[h,1:],coord[h+1,1:],coord[h+1,:-1]],axis=1)
print(faces.shape)

deformed_grid = Mesh([points, faces]).lineWidth(1).alpha(.5)

show(deformed_grid,
     Points(points, c='red6', r=5),
     deformed_grid.labels('id', rotX=-45, scale=.02).c('yellow'),
     bg='k', bg2='lg', axes=9,
)

I'll complete this answer soon with a full deformation visualization tool.
Thanks again for the help.

[marcomusy]

Cool ! that was fast :) for cosmetics you may add deformed_grid.computeNormals() it will make it look smooth.

[antonfrancois]

Yes, it looks better indeed, but in real cases I have much more points and so we may not need to smooth the surface. I work with diffeomorphic deformation, so smooth by nature =)

[antonfrancois]

Hello again,

So here is an implementation of the viewer as I described it in the question above. I am sure we can still improve it. What I have noted is :

• There is strange behavior with the box cutter. It appears only when I close the window, (which does not close the window but make the box cutter appears) then when I re-close the window the cutter exit normally, and I can finally close the window.
  Maybe putting a slider that we can activate/desactivate
  to show only one grid at a time would be a solution
• It would be cool to add a button to save the current view in a selected format (ex :png,pdf)
• I would be super cool to add some streamlines. I guess it would not be to hard to implement, but I didn't had the time to go trough your methods.
• EDIT : The cutter tool do not works in jupyters notebooks. with this error : AttributeError: 'ViewInteractiveWidget' object has no attribute 'addCutterTool'
  Even with the modification to make it works in a jupyter notebook, it seems to do not work inside a ipython evironement.

Don't hesitate to give me some feedback on how to improve the class :)

import vedo
import numpy as np


vedo.embedWindow('ipyvtk')
class deformation_grid3D_vedo:

    def __init__(self,max_resolution=30,color='yellow',alpha=0.5):
        """ When you initialize the class you have the opportunity
        to select some constants that will define the plot.

        :param max_resolution: int or 3sized tuple of int. Maximum number of
        mesh summits in each dimension. If type(max_resolution) == int then
        its value will be the maximum in each dimension.
        :param color: Color of the surfaces
        :param alpha: alpha value of the surfaces
        """
        self.max_resolution = max_resolution
        self.color = color
        self.alpha = alpha

    def _construct_surface_coord(self, deformation_slice):
        """ takes a deformation slice of the selected dimensions and
        prepare it to be plotted by vedo

        :param defomation_slice: [dim_1,dim_2,3]
        :return: an array with point coordinates and face links.
        """
        if deformation_slice.shape[-1] != 3:
            raise TypeError("deformation slice must have shape (dim_1,dim_2,3) got "+
                            str(deformation_slice.shape))
        dim_1,dim_2,_ = deformation_slice.shape

        points = deformation_slice.reshape((dim_1*dim_2,3))

        coord = np.arange(dim_1*dim_2).reshape((dim_1,dim_2))
        faces = np.zeros((dim_1*dim_2,4))
        for d in range(dim_1-1):
            faces[d*dim_2:d*dim_2 + (dim_2-1),:] = np.stack(
                [coord[d,:-1],coord[d,1:],coord[d+1,1:],coord[d+1,:-1]],
                axis=1
            )

        return [points,faces]

    def _slicer(self, deformation, dim, ind):
        """ Slice the deformation in the given index at the given dimension.
        It also subsample the deformation at the max_resolution set by the user
        in init, default value is 30.

        :param deformation: [1,D,H,W,3] numpy array
        :param dim: int \in \{0,1,2\}
        :param ind: int
        :return: sliced deformation
        """
        _,D,H,W,_ =deformation.shape
        if type(self.max_resolution) == np.int:
            self.max_resolution = (self.max_resolution,)*3
        if dim == 0 :
            h_sampler = np.linspace(0,H-1,min(H,self.max_resolution[1]),dtype=np.long)
            w_sampler = np.linspace(0,W-1,min(W,self.max_resolution[2]),dtype=np.long)
            return deformation[0,ind,h_sampler,:,:][:,w_sampler,:]
        elif dim == 1:
            d_sampler = np.linspace(0,D-1,min(D,self.max_resolution[0]),dtype=np.long)
            w_sampler = np.linspace(0,W-1,min(W,self.max_resolution[2]),dtype=np.long)
            return deformation[0,d_sampler,ind,:,:][:,w_sampler,:]
        elif dim == 2:
            d_sampler = np.linspace(0,D-1,min(D,self.max_resolution[0]),dtype=np.long)
            h_sampler = np.linspace(0,H-1,min(H,self.max_resolution[1]),dtype=np.long)
            return deformation[0,d_sampler,:,ind,:][:,h_sampler,:]
        else:
            raise IndexError("dim has to be {0,1,2} got "+str(dim))

    def forward(self,deformation,dim=0,n_surf=5):
        if len(deformation.shape) != 5 and deformation.shape[-1] !=3:
            raise TypeError("deformation shape must be of the form [1,D,H,W,3]",
                            "got array of dim "+ str(deformation.shape))
        if deformation.shape[0] > 1:
            deformation = deformation[0][None]
            print("Warning, only first Batch dimension will be considered")
        _,D,H,W,_ =deformation.shape

        surf_indexes = np.linspace(0,deformation.shape[dim+1]-1,n_surf,
                                   dtype=np.long)
        surfaces = [None] * n_surf
        for i,ind in enumerate(surf_indexes):
            mesh_pts = self._construct_surface_coord(
                self._slicer(deformation,dim,ind)
            )
            surfaces[i] = vedo.Mesh(mesh_pts)

        surf_meshes = vedo.merge(surfaces).lineWidth(1).alpha(self.alpha).computeNormals().c(self.color)
        # plots
        plot = vedo.Plotter()
        plt = plot.show(surf_meshes,
            # Points(surf_meshes.points(), c='red6', r=5),
            # deformed_grid.c('yellow'),
            # bg='k', bg2='lg',
            axes=1,
        ).addCutterTool(surf_meshes,'box') # comment this line for using the class in jupyter notebooks
        return plt


def deformation_grid3D_surfaces(deformation,dim =0,n_surf=5):
    """ function API for `deformation_grid3D_vedo` class

    :param deformation: [1,D,H,W,3] numpy array
    :param dim: Dimension you want to make the surface along.
    :param n_surf: Numbers of surfaces you want to plot.
    :return:
    """
    return deformation_grid3D_vedo().forward(deformation,dim=dim,n_surf=n_surf)

# ______________ TEST _________________________
# Generating a given 3D deformation
Tau,D,H,W = (1,50,150,200)  # we will use a temporal vector field at a given time
x,y,z = np.meshgrid(np.linspace(-1,1,H),
                    np.linspace(-1,1,D),
                    np.linspace(-1,1,W))
grid = np.stack([x,y,z],axis=-1)[None]
v = np.zeros((Tau,D,H,W,3)) # vector field initialization
coef_mat = np.exp(-5*((grid[0,...,0])**2 +
                         (grid[0,...,1])**2 +
                         (grid[0,...,2])**2))[...,None]
v[0] = coef_mat*np.stack(
     (.3*np.ones(grid[0,...,0].shape),
      np.sin(grid[0,...,2]),
     np.cos(grid[0,...,1])),
axis=3)

deformation = (grid + v[0][None]) # deformation shape is [1,D,H,W,3]


deformation_grid3D_surfaces(deformation,dim=1,n_surf=8)

[marcomusy]

Hi, the cutter tool is pretty experimental i'm afraid, and will not work in jupyter..
For sure it would be relatively easy to add a slider (see vedo -r sliders, vedo --search slider), but it also would not work in jupyter.

• you can easily add buttons too
• export screenshots to common format or full 3d interactive scene

import vedo
import numpy as np

# vedo.embedWindow('ipyvtk')
vedo.settings.useDepthPeeling = True # if you use transparency <1

class deformation_grid3D_vedo:

    def __init__(self, max_resolution=30, color='yellow', alpha=1):
        """
        When you initialize the class you have the opportunity
        to select some constants that will define the plot.

        :param max_resolution: int or 3sized tuple of int. Maximum number of
        mesh summits in each dimension. If type(max_resolution) == int then
        its value will be the maximum in each dimension.
        :param color: Color of the surfaces
        :param alpha: alpha value of the surfaces
        """
        self.max_resolution = max_resolution
        self.color = color
        self.alpha = alpha
        self.plot = None
        self.bu = None

    def _construct_surface_coord(self, deformation_slice):
        """ takes a deformation slice of the selected dimensions and
        prepare it to be plotted by vedo

        :param defomation_slice: [dim_1,dim_2,3]
        :return: an array with point coordinates and face links.
        """
        if deformation_slice.shape[-1] != 3:
            raise TypeError("deformation slice must have shape (dim_1,dim_2,3) got "+
                            str(deformation_slice.shape))
        dim_1,dim_2,_ = deformation_slice.shape

        points = deformation_slice.reshape((dim_1*dim_2,3))

        coord = np.arange(dim_1*dim_2).reshape((dim_1,dim_2))
        faces = np.zeros((dim_1*dim_2,4))
        for d in range(dim_1-1):
            faces[d*dim_2:d*dim_2 + (dim_2-1),:] = np.stack(
                [coord[d,:-1],coord[d,1:],coord[d+1,1:],coord[d+1,:-1]],
                axis=1
            )

        return [points,faces]

    def _slicer(self, deformation, dim, ind):
        """
        Slice the deformation in the given index at the given dimension.
        It also subsample the deformation at the max_resolution set by the user
        in init, default value is 30.

        :param deformation: [1,D,H,W,3] numpy array
        :param dim: int \in \{0,1,2\}
        :param ind: int
        :return: sliced deformation
        """
        _,D,H,W,_ =deformation.shape
        if type(self.max_resolution) == np.int:
            self.max_resolution = (self.max_resolution,)*3
        if dim == 0 :
            h_sampler = np.linspace(0,H-1,min(H,self.max_resolution[1]),dtype=np.long)
            w_sampler = np.linspace(0,W-1,min(W,self.max_resolution[2]),dtype=np.long)
            return deformation[0,ind,h_sampler,:,:][:,w_sampler,:]
        elif dim == 1:
            d_sampler = np.linspace(0,D-1,min(D,self.max_resolution[0]),dtype=np.long)
            w_sampler = np.linspace(0,W-1,min(W,self.max_resolution[2]),dtype=np.long)
            return deformation[0,d_sampler,ind,:,:][:,w_sampler,:]
        elif dim == 2:
            d_sampler = np.linspace(0,D-1,min(D,self.max_resolution[0]),dtype=np.long)
            h_sampler = np.linspace(0,H-1,min(H,self.max_resolution[1]),dtype=np.long)
            return deformation[0,d_sampler,:,ind,:][:,h_sampler,:]
        else:
            raise IndexError("dim has to be {0,1,2} got "+str(dim))
    
    def buttonfunc(self):
        self.plot.export("scene.npz") # vedo 3d format, use command line "vedo scene.npz"
        self.plot.screenshot("scene.png")
        vedo.printc("\save exported scene.npz and screenshot to scene.png", c='g')

    def forward(self, deformation,dim=0, n_surf=5):
        if len(deformation.shape) != 5 and deformation.shape[-1] !=3:
            raise TypeError("deformation shape must be of the form [1,D,H,W,3]",
                            "got array of dim "+ str(deformation.shape))
        if deformation.shape[0] > 1:
            deformation = deformation[0][None]
            print("Warning, only first Batch dimension will be considered")
        _,D,H,W,_ =deformation.shape

        surf_indexes = np.linspace(0,deformation.shape[dim+1]-1,n_surf,
                                   dtype=np.long)
        surfaces = [None] * n_surf
        for i,ind in enumerate(surf_indexes):
            mesh_pts = self._construct_surface_coord(
                self._slicer(deformation,dim,ind)
            )
            surfaces[i] = vedo.Mesh(mesh_pts)

        surf_meshes = vedo.merge(surfaces).computeNormals()
        surf_meshes.lineWidth(0.1).alpha(self.alpha).c(self.color).lighting('off')
        
        # plotter instance
        self.plot = vedo.Plotter()
        
        # add a button
        self.bu = self.plot.addButton(self.buttonfunc,
                                      pos=(0.15, 0.05),  # x,y fraction from bottom left corner
                                      states=["export/save"],
                                      c=["w"],
                                      bc=["dg"],  # colors of states
                                      font="times",   # arial, courier, times
                                      size=25,
                                      bold=True,
                                      italic=False,
        )
               
        # add some coloring
        elevations = surf_meshes.points()[:,0]
        surf_meshes.cmap("jet", elevations).addScalarBar("time")
        self.plot.show(surf_meshes, axes=8, bg2='lightblue', viewup='x', interactive=False)
        
        # self.plot.addCutterTool(surf_meshes,'box') # comment this line for using the class in jupyter notebooks
        vedo.interactive()  # stay interactive
        return

def deformation_grid3D_surfaces(deformation,dim =0,n_surf=5):
    """ function API for `deformation_grid3D_vedo` class

    :param deformation: [1,D,H,W,3] numpy array
    :param dim: Dimension you want to make the surface along.
    :param n_surf: Numbers of surfaces you want to plot.
    :return:
    """
    return deformation_grid3D_vedo().forward(deformation,dim=dim,n_surf=n_surf)

# ______________ TEST _________________________
# Generating a given 3D deformation
Tau,D,H,W = (1,50,150,200)  # we will use a temporal vector field at a given time
x,y,z = np.meshgrid(np.linspace(-1,1,H),
                    np.linspace(-1,1,D),
                    np.linspace(-1,1,W))
grid = np.stack([x,y,z],axis=-1)[None]
v = np.zeros((Tau,D,H,W,3)) # vector field initialization
coef_mat = np.exp(-5*((grid[0,...,0])**2 +
                         (grid[0,...,1])**2 +
                         (grid[0,...,2])**2))[...,None]
v[0] = coef_mat*np.stack(
     (.3*np.ones(grid[0,...,0].shape),
      np.sin(grid[0,...,2]),
     np.cos(grid[0,...,1])),
axis=3)

deformation = (grid + v[0][None]) # deformation shape is [1,D,H,W,3]

deformation_grid3D_surfaces(deformation, dim=1, n_surf=8)