Exploring Volumetric data with VTK

Based on this example



In [9]:

    
import os
os.chdir("D:/vtk_data/VTKData/Data")
#os.chdir(r"C:\Users\Diego\Programas\vtk\VTKData\Data")
import vtk



In [10]:

    
pl3d = vtk.vtkMultiBlockPLOT3DReader()

xyx_file = "combxyz.bin"
q_file = "combq.bin"
pl3d.SetXYZFileName(xyx_file)
pl3d.SetQFileName(q_file)
pl3d.SetScalarFunctionNumber(100)
pl3d.SetVectorFunctionNumber(202)
pl3d.Update()

From the example code:

Here we read data from a annular combustor. A combustor burns fuel and air in a gas turbine (e.g., a jet engine) and the hot gas eventually makes its way to the turbine section.



In [11]:

    
blocks = pl3d.GetOutput()
b0 = blocks.GetBlock(0)

Loot at what is inside the data set



In [12]:

    
# uncomment to get long output
# print b0

Set up vtk environment



In [13]:

    
renderer = vtk.vtkRenderer()
render_window = vtk.vtkRenderWindow()
render_window.AddRenderer(renderer)
interactor = vtk.vtkRenderWindowInteractor()
interactor.SetInteractorStyle(vtk.vtkInteractorStyleTrackballCamera())
render_window.SetInteractor(interactor)
renderer.SetBackground(0.2,0.2,0.2)
interactor.Initialize()
render_window.Render()

Look at an outline of the data for context



In [14]:

    
outline = vtk.vtkStructuredGridOutlineFilter()
outline.SetInputData(b0)
outline_mapper = vtk.vtkPolyDataMapper()
outline_mapper.SetInputConnection(outline.GetOutputPort())
outline_actor = vtk.vtkActor()
outline_actor.SetMapper(outline_mapper)
outline_actor.GetProperty().SetColor(1,1,1)
renderer.AddActor(outline_actor)
renderer.ResetCamera()



In [15]:

    
render_window.Render()
interactor.Start()

Scalars

The data contains two scalar arrays and one vector array. We are going to start exploring the data set by looking at the first scalar array "Density". We will use several different methods to look at this data

Points



In [59]:

    
points = vtk.vtkVertexGlyphFilter()
points.SetInputData(b0)
scalars_mapper = vtk.vtkPolyDataMapper()
scalars_actor = vtk.vtkActor()
scalars_actor.SetMapper(scalars_mapper)
scalars_mapper.SetInputConnection(points.GetOutputPort())



In [60]:

    
scalars_mapper.SetScalarModeToUsePointData()
#print scalars_mapper.GetLookupTable()



In [61]:

    
renderer.AddActor(scalars_actor)
render_window.Render()
interactor.Start()
renderer.RemoveActor(scalars_actor)
render_window.Render()

Solid



In [62]:

    
poly = vtk.vtkGeometryFilter()
poly.SetInputData(b0)
#scalars_mapper = vtk.vtkPolyDataMapper()
#scalars_actor = vtk.vtkActor()
#scalars_actor.SetMapper(scalars_mapper)
scalars_mapper.SetInputConnection(poly.GetOutputPort())



In [63]:

    
renderer.AddActor(scalars_actor)
render_window.Render()
interactor.Start()
renderer.RemoveActor(scalars_actor)
render_window.Render()

Sample with plane



In [64]:

    
samples = vtk.vtkPlaneWidget()
samples.SetResolution(24)
samples.SetOrigin(2,-2,26)
samples.SetPoint1(2, 2,26)
samples.SetPoint2(2,-2,32)
samples_pd = samples.GetPolyDataAlgorithm()
samples.SetInteractor(interactor)
samples.SetRepresentationToOutline()
samples.On()



In [65]:

    
probe = vtk.vtkProbeFilter()
probe.SetSourceData(b0)
probe.SetInputConnection(samples_pd.GetOutputPort())
probed_mapper = vtk.vtkPolyDataMapper()
probed_mapper.SetInputConnection(probe.GetOutputPort())
probed_actor = vtk.vtkActor()
probed_actor.SetMapper(probed_mapper)



In [66]:

    
samples.On()
renderer.AddActor(probed_actor)
render_window.Render()
interactor.Start()
samples.Off()
renderer.RemoveActor(probed_actor)

Iso-surfaces



In [67]:

    
min_s, max_s = b0.GetScalarRange()
print min_s, max_s









    



0.197813093662 0.710419237614



In [68]:

    
contours = vtk.vtkContourFilter()
contours.SetInputData(b0)
contours.GenerateValues(5, (min_s,max_s))
contours_mapper = vtk.vtkPolyDataMapper()
contours_mapper.SetInputConnection(contours.GetOutputPort())
contours_actor = vtk.vtkActor()
contours_actor.SetMapper(contours_mapper)



In [69]:

    
renderer.AddActor(contours_actor)
render_window.Render()
interactor.Start()
renderer.RemoveActor(contours_actor)
render_window.Render()

We will add some transparency to be able to look inside



In [70]:

    
contours_actor.GetProperty().SetOpacity(0.5)



In [71]:

    
renderer.AddActor(contours_actor)
render_window.Render()
interactor.Start()
renderer.RemoveActor(contours_actor)
render_window.Render()

Volume Rendering



In [72]:

    
volume_mapper =  vtk.vtkUnstructuredGridVolumeRayCastMapper()
#volume_mapper = vtk.vtkProjectedTetrahedraMapper()
#volume_mapper = vtk.vtkUnstructuredGridVolumeZSweepMapper()

#volume_mapper.SetBlendModeToComposite()  # try others
#volume_mapper.SetBlendModeToMinimumIntensity()
tri_filter = vtk.vtkDataSetTriangleFilter()
tri_filter.SetInputData(b0)
volume_mapper.SetInputConnection(tri_filter.GetOutputPort())
volume = vtk.vtkVolume()
volume.SetMapper(volume_mapper)
volume_property = volume.GetProperty()



In [73]:

    
intensity = vtk.vtkPiecewiseFunction()
intensity.AddPoint(min_s,0.1)
intensity.AddPoint((max_s+min_s)/2.0,  1.0)
intensity.AddPoint(max_s, 1.0)
color = vtk.vtkColorTransferFunction()
color.AddRGBPoint(min_s,0,1,0)
color.AddRGBPoint(max_s,1,0,0)
volume_property.SetScalarOpacity(intensity)
volume_property.SetInterpolationTypeToLinear()
volume_property.SetColor(color)
volume_property.ShadeOff()
volume.VisibilityOn()



In [74]:

    
renderer.AddViewProp(volume)
render_window.Render()
interactor.Start()
renderer.RemoveViewProp(volume)

Vectors

Now we are going to explore the vector data inside the dataset. We will use two methods

Glyphs



In [75]:

    
arrow = vtk.vtkArrowSource()
glyphs = vtk.vtkGlyph3D()
glyphs.SetInputData(b0)
glyphs.SetSourceConnection(arrow.GetOutputPort())



In [76]:

    
glyph_mapper =  vtk.vtkPolyDataMapper()
glyph_mapper.SetInputConnection(glyphs.GetOutputPort())
glyph_actor = vtk.vtkActor()
glyph_actor.SetMapper(glyph_mapper)



In [77]:

    
glyphs.SetVectorModeToUseVector()
glyphs.SetScaleModeToScaleByScalar()
#glyphs.SetScaleModeToDataScalingOff()

glyph_mapper.UseLookupTableScalarRangeOn()



In [78]:

    
renderer.AddActor(glyph_actor)
render_window.Render()
interactor.Start()
renderer.RemoveActor(glyph_actor)
render_window.Render()

Now lets scale by the norm of the vectors



In [79]:

    
glyphs.SetScaleModeToScaleByVector()
glyphs.SetScaleFactor(0.001)
glyphs.SetColorModeToColorByVector()



In [80]:

    
glyphs.Update()
s0,sf = glyphs.GetOutput().GetScalarRange()
lut = vtk.vtkColorTransferFunction()
lut.AddRGBPoint(s0, 1,0,0)
lut.AddRGBPoint(sf, 0,1,0)
glyph_mapper.SetLookupTable(lut)



In [81]:

    
renderer.AddActor(glyph_actor)
render_window.Render()
interactor.Start()
renderer.RemoveActor(glyph_actor)
render_window.Render()

The above graph is very crowded, we can simplify it by showing a random sampleof vectors. We can now also make them bigger.



In [82]:

    
mask = vtk.vtkMaskPoints()
mask.SetInputData(b0)
mask.GenerateVerticesOn()
mask.SetMaximumNumberOfPoints(10000)
mask.SetRandomModeType(1)
mask.RandomModeOn()
glyphs.SetInputConnection(mask.GetOutputPort())
glyphs.SetScaleFactor(0.005)



In [83]:

    
renderer.AddActor(glyph_actor)
render_window.Render()
interactor.Start()
renderer.RemoveActor(glyph_actor)
render_window.Render()

Here we risk loosing important features in the dataset. Another approach is showing only the vectors above a threshold



In [84]:

    
threshold = vtk.vtkThresholdPoints()
threshold.SetInputData(b0)
print b0.GetScalarRange()
threshold.ThresholdByUpper(0.5)
glyphs.SetInputConnection(threshold.GetOutputPort())









    



(0.19781309366226196, 0.710419237613678)



In [85]:

    
renderer.AddActor(glyph_actor)
render_window.Render()
interactor.Start()
renderer.RemoveActor(glyph_actor)
render_window.Render()

Streamlines



In [86]:

    
seeds = vtk.vtkPlaneWidget()
seeds.SetResolution(8)
seeds.SetOrigin(2,-2,26)
seeds.SetPoint1(2, 2,26)
seeds.SetPoint2(2,-2,32)
seeds_pd = seeds.GetPolyDataAlgorithm()
seeds.SetInteractor(interactor)
seeds.On()



In [87]:

    
streamline = vtk.vtkStreamLine()
streamline.SetInputData(pl3d.GetOutput().GetBlock(0))
streamline.SetSourceConnection(seeds_pd.GetOutputPort())
streamline.SetMaximumPropagationTime(200)
streamline.SetIntegrationStepLength(.2)
streamline.SetStepLength(0.001)
streamline.SetNumberOfThreads(1)
streamline.SetIntegrationDirectionToForward()
streamline.VorticityOn()

streamline_mapper = vtk.vtkPolyDataMapper()
streamline_mapper.SetInputConnection(streamline.GetOutputPort())
streamline_actor = vtk.vtkActor()
streamline_actor.SetMapper(streamline_mapper)
streamline_actor.VisibilityOn()



In [88]:

    
renderer.AddActor(streamline_actor)
seeds.On()
render_window.Render()
interactor.Start()
renderer.RemoveActor(streamline_actor)
render_window.Render()
seeds.Off()

Stream Tubes



In [89]:

    
tubes =vtk.vtkTubeFilter()
tubes.SetInputConnection(streamline.GetOutputPort())
tubes.SetRadius(0.1)
tubes.SetNumberOfSides(8)
tubes.SetVaryRadiusToVaryRadiusByScalar()
tubes.SetRadiusFactor(3)

tubes_mapper = vtk.vtkPolyDataMapper()
tubes_mapper.SetInputConnection(tubes.GetOutputPort())
tubes_actor = vtk.vtkActor()
tubes_actor.SetMapper(tubes_mapper)



In [90]:

    
renderer.AddActor(tubes_actor)
seeds.On()
render_window.Render()
interactor.Start()
renderer.RemoveActor(tubes_actor)
render_window.Render()
seeds.Off()

We have looked at several techniques one at a time. Now you could try mixing them up to show a complete visualization. Also up to now we have ignored the third scalar in the data...

Additional sources



In [90]: