I need help with vtkImageStencil.
The attached program is doing the following:
- loading a 3D image (CT-scan data)
- applying a threshold
- creating polydata via vtkFlyingEdges3D (MarchingCube alg.)
- using vtkPolyDataToImageStencil to create a stencil
- applying the stencil (vtkImageStencil) algorithm to create a 3D mask (input is an empty vtkImageGridSource)
Finally I plot everyting with appropriate mappers to a render scene along with the original image (for which the threshold is applied)
Unfortunately the plot shows that mask (black) at exactly half the size as the original image (red). Although it has the exact same settings regarding the size, spacing etc.
Does anybody have an idea where I went wrong. I guess it must have something to do how I set up the input image with vtkImageGridSource?
Any help is greatly appreciated!
import vtk
from vtkmodules.vtkCommonDataModel import vtkSphere
colors = vtk.vtkNamedColors()
reader = vtk.vtkNrrdReader()
# Start by loading some data.
filename = "data\manually matched\Drydisk07.6.nrrd"
reader.SetFileName(filename)
reader.Update()
# An outline is shown for context.
original_outline_cube = vtk.vtkOutlineFilter()
original_outline_cube.SetInputConnection(reader.GetOutputPort())
original_outline_cube_mapper = vtk.vtkPolyDataMapper()
original_outline_cube_mapper.SetInputConnection(original_outline_cube.GetOutputPort())
original_outline_actor = vtk.vtkActor()
original_outline_actor.SetMapper(original_outline_cube_mapper)
# create a threshed image
image_threshed = vtk.vtkImageThreshold()
image_threshed.SetInputConnection(reader.GetOutputPort()) # Set your input vtkImageData
image_threshed.ThresholdByUpper(450) # Set the threshold value
image_threshed.ReplaceInOn() # Set the operation to replace in values
image_threshed.SetInValue(1) # Set the value for inside the threshold
image_threshed.ReplaceOutOn() # Set the operation to replace out values
image_threshed.SetOutValue(0) # Set the value for outside the threshold
# renderer
ren = vtk.vtkRenderer()
# render_window
ren_win = vtk.vtkRenderWindow()
ren_win.AddRenderer(ren)
# render_interactor
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(ren_win)
# Convert the image data to a surface mesh using vtkMarchingCubes
# marchingCubes = vtk.vtkMarchingCubes()
marching_cubes = vtk.vtkFlyingEdges3D()
marching_cubes.SetInputConnection(image_threshed.GetOutputPort())
marching_cubes.SetValue(0, 1) # Threshold value
# use clipped polydata as stencil
polydata_to_image_stencil = vtk.vtkPolyDataToImageStencil()
polydata_to_image_stencil.SetInputConnection(marching_cubes.GetOutputPort())
# Create an image stencil
imageStencil = vtk.vtkImageStencil()
image_grid_source = vtk.vtkImageGridSource()
imageStencil.SetInputConnection(image_grid_source.GetOutputPort())
imageStencil.SetStencilConnection(polydata_to_image_stencil.GetOutputPort())
imageStencil.SetBackgroundValue(1000)
imageStencil.ReverseStencilOn() # Set to cut out the data inside the stencil
# Visualize the cut-out image data
cutout_image_mapper = vtk.vtkFixedPointVolumeRayCastMapper()
# Get the output of the stencil
cutout_image_mapper.SetInputConnection(imageStencil.GetOutputPort())
original_mapper = vtk.vtkFixedPointVolumeRayCastMapper()
original_mapper.SetInputConnection(image_threshed.GetOutputPort())
# Create transfer mapping scalar value to opacity for visualisation
original_opacity_transfer_function = vtk.vtkPiecewiseFunction()
original_opacity_transfer_function.AddPoint(0, 0.0)
original_opacity_transfer_function.AddPoint(450, 1)
original_opacity_transfer_function.AddPoint(451, 1)
opacityThreshedTransferFunction = vtk.vtkPiecewiseFunction()
opacityThreshedTransferFunction.AddPoint(0, 0)
opacityThreshedTransferFunction.AddPoint(1, 1)
opacityThreshedTransferFunction.AddPoint(2500, 1)
# Create transfer mapping scalar value to color.
original_color_transfer_function = vtk.vtkColorTransferFunction()
original_color_transfer_function.AddRGBPoint(0.0, 0.0, 0.0, 0.0)
original_color_transfer_function.AddRGBPoint(255.0, 1.0, 1.0, 1.0)
threshed_color_transfer_function = vtk.vtkColorTransferFunction()
threshed_color_transfer_function.AddRGBPoint(0.0, 1.0, 0.0, 0.0)
threshed_color_transfer_function.AddRGBPoint(255.0, 1.0, 0.0, 0.0)
# # The property describes how the data will look.
original_volume_property = vtk.vtkVolumeProperty()
original_volume_property.SetColor(original_color_transfer_function)
original_volume_property.SetScalarOpacity(original_opacity_transfer_function)
original_volume_property.ShadeOn()
original_volume_property.SetInterpolationTypeToLinear()
cutOutImageActor = vtk.vtkVolume()
cutOutImageActor.SetMapper(cutout_image_mapper)
cutOutImageActor.SetProperty(original_volume_property)
originalVolumeProperty = vtk.vtkVolumeProperty()
originalVolumeProperty.SetColor(threshed_color_transfer_function)
originalVolumeProperty.SetScalarOpacity(opacityThreshedTransferFunction)
originalVolumeProperty.ShadeOn()
originalVolumeProperty.SetInterpolationTypeToLinear()
originalActor = vtk.vtkVolume()
originalActor.SetMapper(original_mapper)
originalActor.SetProperty(originalVolumeProperty)
# add the actors to the scene
ren.AddVolume(cutOutImageActor)
ren.AddVolume(originalActor)
ren.AddActor(original_outline_actor)
ren.SetBackground(colors.GetColor3d('Wheat'))
ren.GetActiveCamera().Azimuth(45)
ren.GetActiveCamera().Elevation(30)
interactorStyle3D_2 = vtk.vtkInteractorStyleTrackballCamera()
iren.SetInteractorStyle(interactorStyle3D_2)
image_threshed.Update()
extent = reader.GetExecutive().GetWholeExtent(reader.GetOutputInformation(0))
image_grid_source.SetDataOrigin(reader.GetDataOrigin())
image_grid_source.SetDataSpacing(reader.GetDataSpacing())
image_grid_source.SetDataExtent(reader.GetDataExtent())
image_grid_source.Modified()
image_grid_source.Update()
imageStencil.GetOutput().SetSpacing(reader.GetDataSpacing())
imageStencil.GetOutput().SetExtent(reader.GetDataExtent())
imageStencil.GetOutput().SetOrigin(reader.GetDataOrigin())
# Set the scalar type and number of scalar components
imageStencil.Update()
ren.ResetCameraClippingRange()
ren.ResetCamera()
ren_win.Render()
iren.Start()
