streamlines filter

[banesullivan]

Resolve #198

This adds a .streamlines() filter that will integrate a vector field to generate streamlines.

Todo:

• Document the filter
• Add parameter descriptions

@JiaweiZhuang: would you look over these changes and provide feedback?

Example A:

import vtki
from vtki import examples
mesh = examples.download_carotid()

streamlines = mesh.streamlines(max_time=100.0,
                           initial_step_length=2., terminal_speed=0.1,
                           n_points=25, source_radius=2.0,
                           source_center=(133.1, 116.3, 5.0) )

p = vtki.Plotter(notebook=1)
p.add_mesh(mesh.outline(), color='k')
p.add_mesh(streamlines.tube(radius=0.15))
p.add_mesh(mesh.contour([160]).wireframe(),
           color='grey', opacity=0.25)
p.camera_position = [(182., 177., 50),
                     (139, 105, 19),
                     (-0.2, -0.2, 1)]
p.show()

Example B:

import vtki
from vtki import examples

mesh = examples.download_blood_vessels().cell_data_to_point_data()
mesh.set_active_scalar('velocity')

streamlines, src = mesh.streamlines(return_source=True, source_radius=10,
                           source_center=(92.46, 74.37, 135.5) ,)

boundary = mesh.decimate_boundary()

p = vtki.Plotter(notebook=1)
p.add_mesh(streamlines.tube(radius=0.2), ligthing=False)
p.add_mesh(src)
p.add_mesh(boundary, color='grey', opacity=0.25)
p.camera_position = [(10, 9.5, -43),
                     (87.0, 73.5, 123.0),
                     (-0.5, -0.7, 0.5)]
p.show()

[banesullivan]

Note that I had to convert the blood vessel dataset to point data:

mesh = examples.download_blood_vessels().cell_data_to_point_data()

I have no idea why, but this implementation currently does not handle vector field on the cells of a mesh... ParaView has no issue with it so I wonder what's different? Maybe they're running a cell data to point data filter behind the scenes?

[JiaweiZhuang]

Thanks for such fast work!

would you look over these changes and provide feedback?

I am able to reproduce your plot at test_streamlines.ipynb Looks great. No major things to add.

I had to convert the blood vessel dataset to point data

Yeah that's weird. mesh.cell_centers().streamlines() also gives empty data. "Cell data" vs "point data" is something that confuses me a lot in VTK. While mesh.points is the 3D coordinate of grid points, mesh.cells is an array of integer indices, referring to the point array. It's size is 9 * n_cells which I am not sure how to interpret (see the below screenshot from the notebook). Maybe such information is insufficient for computing the streamlines.

[banesullivan]

I want to point out that this is an awesome demonstration of vtki's benefit - "Example A" above using the carotid dataset is a replica of this example from VTK. Note that with just the VTK package, this takes about 80 lines of code while with vtki, this takes 10 lines of code to make the same exact 3D scene (not including comments or line breaks for code style/readability)

[banesullivan]

@JiaweiZhuang:

mesh.cell_centers().streamlines() also gives empty data ... Maybe such information is insufficient for computing the streamlines.

That should be expected because the cell_centers() filter returns vertex cells of all the cell center locations - a vertex cell is just a point in 3D space which doesn't have an area or volume associated with it - thus it cannot be integrated by the vtkStreamTracer filter

Cell data" vs "point data" is something that confuses me a lot in VTK. While mesh.points is the 3D coordinate of grid points, mesh.cells is an array of integer indices, referring to the point array. It's size is 9 * n_cells which I am not sure how to interpret (see the below screenshot from the notebook).

This is fair - we haven't documented the difference between point and cell data very well in vtki and it's not very intuitive to new VTK users (hence the reason for #161). In brief, the difference comes down to where in space your scalar/vector values live - cell data means the entire cell (2D triangle, 3D voxel, 3D hexahedron, etc.) is associated with your data (scalar array or vector), while point data means the scalars or vectors live on the nodes of the mesh and not the cells (nodes are the points that make up a cell's extent). When rendering, if plotting the cell data, the entire cell's color is mapped to its value and when plotting point data, cells colors are interpolations of the values at all the nodes of that cell:

import vtki
from vtki import examples

mesh = examples.load_uniform()

p = vtki.Plotter(shape=(1,2))
p.add_mesh(mesh, scalars='Spatial Cell Data', show_edges=True)
p.subplot(0,1)
p.add_mesh(mesh, scalars='Spatial Point Data', show_edges=True)
p.show()

In vtki, we've mostly kept this hidden from users to avoid confusion but the difference is really important for VTK and when trying to visualize data resolution.

As for the points and cells properties on the meshes, this is something that needs updating in the docs... the points property is simply the nodes of the mesh in 3D space - the cells simply connect the points in a specific order to create a volume or surface, etc.

As for why the cells array is 9 times the length of your number of cells, this is because you have a mesh where all of the cells are of the same type and defined by 8 node locations (the points of the voxel) - the 9th element comes from having to separate the cell indices to make individual cells.

The mesh.cells array comes directly from VTK and has a structure where the first value (8 in your case) tells you how many indices make that first cell (8 indices) then the following n (8 again) elements in that array are the indices of the points/nodes array that make that cell.

note that the cells property is only available on vtki.UnstructuredGrid class - why might want to make a more universal implementation

# using init snippet above
mesh = mesh.cast_to_unstructured_grid()

print(mesh.cells[0:9])

[  8   0   1  10  11 100 101 110 111]

last_cells_nodes = mesh.points[mesh.cells[-8:]]

p = vtki.Plotter()
p.add_mesh(mesh, show_edges=True)
p.add_mesh(last_cells_nodes, point_size=10)
p.show()

[JiaweiZhuang]

@banesullivan Thanks for this detailed explanation; this makes it much clearer!

cell data means the entire cell (2D triangle, 3D voxel, 3D hexahedron, etc.) is associated with your data (scalar array or vector), while point data means the scalars or vectors live on the nodes of the mesh and not the cells (nodes are the points that make up a cell's extent).

Yeah this is also what I imagined. Point data are like finite difference grids and cell data are like finite volume grids.

As for why the cells array is 9 times the length of your number of cells, this is because you have a mesh where all of the cells are of the same type and defined by 8 node locations (the points of the voxel) - the 9th element comes from having to separate the cell indices to make individual cells.

Thanks for the clarification. I naively expected that mesh.cells would also return some sort of coordinate values, like mesh.points does.