Skip to content

Instantly share code, notes, and snippets.

@lassoan
Created March 16, 2018 02:40
Show Gist options
  • Save lassoan/74b5efb699b1c9d3a6ef76d98c77b697 to your computer and use it in GitHub Desktop.
Save lassoan/74b5efb699b1c9d3a6ef76d98c77b697 to your computer and use it in GitHub Desktop.
Reconstruct open surface from labelmap
# 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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment