Skip to content

Instantly share code, notes, and snippets.

@somada141
Last active May 1, 2021 05:24
Show Gist options
  • Save somada141/b125fd74916018ffe028 to your computer and use it in GitHub Desktop.
Save somada141/b125fd74916018ffe028 to your computer and use it in GitHub Desktop.
Volume rendering in Python using VTK-SimpleITK #python #visualization #vtk #simpleitk #itk
#!/usr/bin/python
import SimpleITK as sitk
import vtk
import numpy as np
from vtk.util.vtkConstants import *
def numpy2VTK(img,spacing=[1.0,1.0,1.0]):
# evolved from code from Stou S.,
# on http://www.siafoo.net/snippet/314
importer = vtk.vtkImageImport()
img_data = img.astype('uint8')
img_string = img_data.tostring() # type short
dim = img.shape
importer.CopyImportVoidPointer(img_string, len(img_string))
importer.SetDataScalarType(VTK_UNSIGNED_CHAR)
importer.SetNumberOfScalarComponents(1)
extent = importer.GetDataExtent()
importer.SetDataExtent(extent[0], extent[0] + dim[2] - 1,
extent[2], extent[2] + dim[1] - 1,
extent[4], extent[4] + dim[0] - 1)
importer.SetWholeExtent(extent[0], extent[0] + dim[2] - 1,
extent[2], extent[2] + dim[1] - 1,
extent[4], extent[4] + dim[0] - 1)
importer.SetDataSpacing( spacing[0], spacing[1], spacing[2])
importer.SetDataOrigin( 0,0,0 )
return importer
def volumeRender(img, tf=[],spacing=[1.0,1.0,1.0]):
importer = numpy2VTK(img,spacing)
# Transfer Functions
opacity_tf = vtk.vtkPiecewiseFunction()
color_tf = vtk.vtkColorTransferFunction()
if len(tf) == 0:
tf.append([img.min(),0,0,0,0])
tf.append([img.max(),1,1,1,1])
for p in tf:
color_tf.AddRGBPoint(p[0], p[1], p[2], p[3])
opacity_tf.AddPoint(p[0], p[4])
# working on the GPU
# volMapper = vtk.vtkGPUVolumeRayCastMapper()
# volMapper.SetInputConnection(importer.GetOutputPort())
# # The property describes how the data will look
# volProperty = vtk.vtkVolumeProperty()
# volProperty.SetColor(color_tf)
# volProperty.SetScalarOpacity(opacity_tf)
# volProperty.ShadeOn()
# volProperty.SetInterpolationTypeToLinear()
# working on the CPU
volMapper = vtk.vtkVolumeRayCastMapper()
compositeFunction = vtk.vtkVolumeRayCastCompositeFunction()
compositeFunction.SetCompositeMethodToInterpolateFirst()
volMapper.SetVolumeRayCastFunction(compositeFunction)
volMapper.SetInputConnection(importer.GetOutputPort())
# The property describes how the data will look
volProperty = vtk.vtkVolumeProperty()
volProperty.SetColor(color_tf)
volProperty.SetScalarOpacity(opacity_tf)
volProperty.ShadeOn()
volProperty.SetInterpolationTypeToLinear()
# Do the lines below speed things up?
# pix_diag = 5.0
# volMapper.SetSampleDistance(pix_diag / 5.0)
# volProperty.SetScalarOpacityUnitDistance(pix_diag)
vol = vtk.vtkVolume()
vol.SetMapper(volMapper)
vol.SetProperty(volProperty)
return [vol]
def vtk_basic( actors ):
"""
Create a window, renderer, interactor, add the actors and start the thing
Parameters
----------
actors : list of vtkActors
Returns
-------
nothing
"""
# create a rendering window and renderer
ren = vtk.vtkRenderer()
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(ren)
renWin.SetSize(600,600)
# ren.SetBackground( 1, 1, 1)
# create a renderwindowinteractor
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
for a in actors:
# assign actor to the renderer
ren.AddActor(a )
# render
renWin.Render()
# enable user interface interactor
iren.Initialize()
iren.Start()
#####
filename = 'toto.nii.gz'
img = sitk.ReadImage( filename ) # SimpleITK object
data = sitk.GetArrayFromImage( img ) # numpy array
from scipy.stats.mstats import mquantiles
q = mquantiles(data.flatten(),[0.7,0.98])
q[0]=max(q[0],1)
q[1] = max(q[1],1)
tf=[[0,0,0,0,0],[q[0],0,0,0,0],[q[1],1,1,1,0.5],[data.max(),1,1,1,1]]
actor_list = volumeRender(data, tf=tf, spacing=img.GetSpacing())
vtk_basic(actor_list)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment