The following code creates image data with two regions, where one region (ID 2) is surrounded by background values (ID 0) on all sides except for one, where there is a single shared boundary with another region (ID 5).
Using vtkSurfaceNets3D, it’s possible to extract this internal boundary using a combination of SetOutputStyleToSelected, SetLabel, and AddSelectedLabel.
With VTK 9.6.2, the code below generates this output (using contour_old(image)):
where the 1 component of the two-component ‘BoundaryLabels’ array is shown as colored cells. The blue cells show the internal boundary [2, 5] (plotted as value 5), whereas the red cells show the external boundary with the background [2, 0] (plotted as value 0).
With the latest VTK built from master, we still get the same result using contour_old(image), but there is a deprecation warning for SetOutputStyleToSelected and AddSelectedLabel, suggesting to use vtkSurfaceNetsAtlas instead. Notably SetLabel is not deprecated, but maybe it should be? I tried recreating the same result using the new API (and still using the old SetLabel on vtkSurfaceNets3D), but I am unable to reproduce the result using the contour_new(image) code below. This is the result I get instead, with the internal boundary between regions missing, such that we now see “inside” the contour.
Am I missing something, or was this perhaps an oversight for the deprecation, and the deprecation for SetLabel was missed?
Maybe @spyridon97 can comment on this?
import vtk
def create_labeled_image():
# Create 4x3x3 image with two adjacent labels
# First label (ID 2):
# has a single point near center of image,
# is adjacent to second label,
# is otherwise surrounded by background,
# Second label (ID 5):
# has two points near center of image,
# is adjacent to first label,
# has one side touching image boundary,
# is otherwise surrounded by background
dim = (4, 3, 3)
image = vtk.vtkImageData()
image.SetDimensions(dim)
n_points = dim[0] * dim[1] * dim[2]
labels = vtk.vtkIntArray()
labels.SetName("labels")
labels.SetNumberOfComponents(1)
labels.SetNumberOfTuples(n_points)
for i in range(n_points):
labels.SetValue(i, 0)
labels.SetValue(17, 2)
labels.SetValue(18, 5)
labels.SetValue(19, 5)
image.GetPointData().AddArray(labels)
image.GetPointData().SetActiveScalars("labels")
return image
def contour_old(image):
surface_nets = vtk.vtkSurfaceNets3D()
surface_nets.SetInputData(image)
surface_nets.SmoothingOff()
surface_nets.AddSelectedLabel(2) # Moved to atlas in 9.7
surface_nets.SetOutputStyleToSelected() # Moved to atlas in 9.7
surface_nets.SetLabel(2, 2) # Missing from atlas in 9.7 (?)
surface_nets.SetLabel(5, 5)
surface_nets.Update()
return surface_nets.GetOutput()
def contour_new(image):
surface_nets = vtk.vtkSurfaceNets3D()
surface_nets.SetInputData(image)
surface_nets.SmoothingOff()
surface_nets.SetLabel(2, 2)
surface_nets.SetLabel(5, 5)
surface_nets.Update()
atlas = vtk.vtkSurfaceNetsAtlas()
atlas.SetInputDataObject(surface_nets.GetOutput())
atlas.SetExtractionModeToLabelSet()
atlas.AddSelectedLabel(2)
atlas.SetOutputStyleToBoundary()
atlas.Update()
pdc = atlas.GetOutput()
blocks = [
pds.GetPartition(j)
for i in range(pdc.GetNumberOfPartitionedDataSets())
for pds in [pdc.GetPartitionedDataSet(i)]
for j in range(pds.GetNumberOfPartitions())
if pds.GetPartition(j) is not None
]
append = vtk.vtkAppendPolyData()
for block in blocks:
append.AddInputData(block)
append.Update()
return append.GetOutput()
def plot(poly):
# Mapper
mapper = vtk.vtkPolyDataMapper()
mapper.SetInputData(poly)
mapper.SetScalarModeToUseCellData()
mapper.SelectColorArray("BoundaryLabels")
mapper.ScalarVisibilityOn()
mapper.SetColorModeToMapScalars()
mapper.SetArrayComponent(1)
# Actor
actor = vtk.vtkActor()
actor.SetMapper(mapper)
# Renderer
renderer = vtk.vtkRenderer()
renderer.AddActor(actor)
renderer.SetBackground(0.2, 0.3, 0.4)
# Render window
render_window = vtk.vtkRenderWindow()
render_window.AddRenderer(renderer)
render_window.SetSize(800, 600)
# Interactor
interactor = vtk.vtkRenderWindowInteractor()
interactor.SetRenderWindow(render_window)
renderer.ResetCamera()
camera = renderer.GetActiveCamera()
camera.SetPosition(1, 1, 1)
camera.SetFocalPoint(0, 0, 0)
camera.SetViewUp(0, 0, 1)
renderer.ResetCamera()
render_window.Render()
# Start rendering
render_window.Render()
interactor.Start()
image = create_labeled_image()
# poly = contour_old(image)
poly = contour_new(image)
plot(poly)
EDIT: Fix colors that are referenced in plot
