Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- # 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()
Add Comment
Please, Sign In to add comment