Reconstruct open surface from labelmap

ReconstructSurface.py

# Important: the input volume must have isotropic spacing.
# If surface is thick, use Extract skeleton module with Skeleton type = 2D, Do not prune branches = enabled.
# Increase radius parameter to fill more holes.
# Increase dimension to preserve more details.

inputLabelmap = getNode('Input labelmap')

ici = vtk.vtkImageChangeInformation()
ici.SetInputData(inputLabelmap.GetImageData())
ici.SetOutputSpacing(inputLabelmap.GetSpacing())
ici.SetOrigin(inputLabelmap.GetOrigin())

imageToPoints=vtk.vtkImageDataToPointSet()
imageToPoints.SetInputConnection(ici.GetOutputPort())

geom = vtk.vtkStructuredGridGeometryFilter()
geom.SetInputConnection(imageToPoints.GetOutputPort())

threshold = vtk.vtkThreshold()
threshold.ThresholdByUpper(1)
threshold.SetInputConnection(geom.GetOutputPort())

geom2 = vtk.vtkGeometryFilter()
geom2.SetInputConnection(threshold.GetOutputPort())
geom2.Update()
polyData=geom2.GetOutput()

# Estimate normals
normals = vtk.vtkPCANormalEstimation()
normals.SetInputData(polyData)
sampleSize = polyData.GetNumberOfPoints() * .00005
if sampleSize < 10:
  sampleSize = 50;
normals.SetSampleSize(int(sampleSize))
normals.SetNormalOrientationToGraphTraversal()
normals.FlipNormalsOn()

distance = vtk.vtkSignedDistance()
distance.SetInputConnection(normals.GetOutputPort())

dimension = 256
bounds = [0,-1,0,-1,0,-1]
polyData.GetBounds(bounds);
range = [bounds[1]-bounds[0], bounds[3]-bounds[2], bounds[5]-bounds[4]]
radius = max(range) / float(dimension) * 4  # approximately 4 voxels

distance.SetRadius(radius)
distance.SetDimensions(dimension, dimension, dimension)
distance.SetBounds(
  bounds[0] - range[0] * .1,
  bounds[1] + range[0] * .1,
  bounds[2] - range[1] * .1,
  bounds[3] + range[1] * .1,
  bounds[4] - range[2] * .1,
  bounds[5] + range[2] * .1)

surface = vtk.vtkExtractSurface()
surface.SetInputConnection(distance.GetOutputPort())
surface.SetRadius(radius * .99)
surface.Update()

resultSurface = slicer.modules.models.logic().AddModel(surface.GetOutput())
resultSurface.GetDisplayNode().BackfaceCullingOff()