Server side data processing in VTK

  • simple tools on client side (just meshes)
  • we use VTK to process volumetric mesh on server side
  • ipywidgets interactively show results as 3d surface

In [1]:
import k3d
import os
import vtk
from vtk.util import numpy_support

import numpy as np

In [2]:
plot = k3d.plot(camera_auto_fit=False)
plot.display()



In [3]:
filename = './cutter/output_fem.vtu'
reader = vtk.vtkXMLUnstructuredGridReader()
reader.SetFileName(filename)
reader.Update()

In [4]:
grid = reader.GetOutput()
bbox = np.array(grid.GetBounds()).reshape(3,2)
center = np.mean(bbox,axis=1)
center


Out[4]:
array([0.04984861, 0.03934663, 0.04888905])

In [5]:
plane = vtk.vtkPlane()
plane.SetOrigin(*center)
plane.SetNormal(1,0.3,0)

In [6]:
def vtk_ExtractSurface(vtk_outputport,vtk_o,vtk_n):

    plane.SetOrigin(*vtk_o)
    plane.SetNormal(*vtk_n)

    myExtractGeometry = vtk.vtkExtractGeometry()
    myExtractGeometry.SetInputConnection(vtk_outputport)
    myExtractGeometry.SetImplicitFunction(plane)
    myExtractGeometry.ExtractInsideOn()
    myExtractGeometry.SetExtractBoundaryCells(0);
    myExtractGeometry.Update()

    myExtractSurface = vtk.vtkDataSetSurfaceFilter()
    myExtractSurface.SetInputConnection(myExtractGeometry.GetOutputPort())
    myExtractSurface.Update()
    
    return myExtractSurface.GetOutput()

In [7]:
def update_from_cut(reader,vtk_o,vtk_n,plt_vtk):
    poly_data =  vtk_ExtractSurface(reader.GetOutputPort(),vtk_o,vtk_n)
    if poly_data.GetNumberOfCells() > 0:
        vertices,indices,attribute = get_mesh_data(poly_data)
        with plt_vtk.hold_sync():
            plt_vtk.vertices = vertices
            plt_vtk.indices = indices
            plt_vtk.attribute = attribute

In [8]:
def get_mesh_data(poly_data,color_attribute=('Umag', 0.0, 0.1)):

    if poly_data.GetPolys().GetMaxCellSize() > 3:
            cut_triangles = vtk.vtkTriangleFilter()
            cut_triangles.SetInputData(poly_data)
            cut_triangles.Update()
            poly_data = cut_triangles.GetOutput()

    if color_attribute is not None:
        attribute = numpy_support.vtk_to_numpy(poly_data.GetPointData().GetArray(color_attribute[0]))
        color_range = color_attribute[1:3]
    else:
        attribute = []
        color_range = []

    vertices = numpy_support.vtk_to_numpy(poly_data.GetPoints().GetData())
    indices = numpy_support.vtk_to_numpy(poly_data.GetPolys().GetData()).reshape(-1, 4)[:, 1:4]

    return (np.array(vertices, np.float32),np.array(indices, np.uint32),np.array(attribute, np.float32))

In [9]:
def clipping_plane_to_vtkPlane(clipping_plane):
    vtk_n = -np.array(clipping_plane[:3])
    vtk_o = clipping_plane[3]*vtk_n
    return (vtk_o,vtk_n)

In [10]:
plane = vtk.vtkPlane()
plane.SetOrigin(*center)
plane.SetNormal(1,0.3,0)

vtk_n = np.array([ 0. ,  .3,  0. ])
vtk_o = np.array([ 0.04984861,  20.03934663,  0.04888905])

In [11]:
plt_vtk = k3d.vtk_poly_data(
    vtk_ExtractSurface(
        reader.GetOutputPort(),
        vtk_o,vtk_n
    ), 
    color_attribute=('Umag', 0.0, 0.1), 
    color_map = k3d.basic_color_maps.Jet)
plot += plt_vtk

In [12]:
update_from_cut(reader,center+0.0,[1,0,0], plt_vtk)

In [13]:
from ipywidgets import FloatSlider, interact
@interact(s = FloatSlider(min=-0.1,max=0.1,step=0.00004))
def _(s):
    update_from_cut(reader,center+s,[1,0,0],plt_vtk)



In [14]:
plt_vtk.color_map = k3d.colormaps.paraview_color_maps.Coolwarm
plt_vtk.color_range = [0,0.32]
plt_vtk.flat_shading = True

In [15]:
plot.camera= [0.07, 0.03, 0.05, 0.050, 0.039, 0.046, -0.04, 0.063, 0.997]

In [16]:
plt_mesh = k3d.vtk_poly_data(vtk_ExtractSurface(reader.GetOutputPort(),vtk_o,vtk_n))

plt_mesh.wireframe = True
plt_mesh.color = 0xaaaaaa
plt_mesh.opacity = 0.2

plot += plt_mesh

In [ ]: