[vtkusers] Intersection of plane with 3D surface using vtkCutter
Daniël Lacko
lacko.daniel at outlook.com
Thu Aug 1 11:46:53 EDT 2013
Hi all,
I'm trying to calculate the circumference of the human head in 3D. I'm using VTK 5.10.1 and Python 2.7 64-bit. I want to calculate the intersection of a plane passing through certain point on the head. I then want to extract the boundary and calculate the length of this polyline. First, I want to visualize the boundary to make sure the extraction is working as expected.
As far as I understand, I can use vtkCutter to cut a slice from the original image. When I modify the vtkCutter ClipCow.py example, I can use vtkFeatureExtract to extract the boundary and write it to a .vtk file for visualization:
# Here we are cutting the cow. Cutting creates lines where the cut
# function intersects the model. (Clipping removes a portion of the
# model but the dimension of the data does not change.)
#
# The reason we are cutting is to generate a closed polygon at the
# boundary of the clipping process. The cutter generates line
# segments, the stripper then puts them together into polylines. We
# then pull a trick and define polygons using the closed line
# segments that the stripper created.
cutEdges = vtk.vtkCutter()
cutEdges.SetInputConnection(cow.GetOutputPort())
cutEdges.SetCutFunction(plane)
cutEdges.GenerateCutScalarsOn()
cutEdges.SetValue(0, 0.5)
cutStrips = vtk.vtkStripper()
cutStrips.SetInputConnection(cutEdges.GetOutputPort())
cutStrips.Update()
cutPoly = vtk.vtkPolyData()
cutPoly.SetPoints(cutStrips.GetOutput().GetPoints())
cutPoly.SetPolys(cutStrips.GetOutput().GetLines())
# Extract boundary from cutPoly
cutBoundary = vtk.vtkFeatureEdges()
cutBoundary.SetInput(cutPoly)
cutBoundary.Update()
# Write cut line to .vtk file
polyWriter = vtk.vtkPolyDataWriter()
polyWriter.SetInput(cutBoundary.GetOutput())
polyWriter.SetFileName("cutPoly.vtk")
polyWriter.Write()
The code works as intended. However, if I try using it with my own code, it won't read the surface data correctly:
# Circumference plane (used to cut)
circPlaneSource = vtk.vtkPlaneSource()
circPlaneSource.SetOrigin(G -
150*self.midSagittalPlaneNormal)
circPlaneSource.SetPoint1(Op -
150*self.midSagittalPlaneNormal)
circPlaneSource.SetPoint2(G +
150*self.midSagittalPlaneNormal)
circPlaneSource.SetResolution(15, 15)
circPlaneSource.Update()
surface = vtk.vtkSTLReader()
surface.SetFileName(fileName)
# Cut a slice out of the surface and convert it to a PolyData object
cutEdges = vtk.vtkCutter()
cutEdges.SetInputConnection(surface.GetOutputPort())
cutEdges.SetCutFunction(circPlaneSource.GetOutput())
cutEdges.GenerateCutScalarsOn()
cutEdges.SetValue(0, 0.5)
cutStrips = vtk.vtkStripper()
cutStrips.SetInputConnection(cutEdges.GetOutputPort())
cutStrips.Update()
cutPoly = vtk.vtkPolyData()
cutPoly.SetPoints(cutStrips.GetOutput().GetPoints())
cutPoly.SetPolys(cutStrips.GetOutput().GetLines())
# Extract boundary from cutPoly
cutBoundary = vtk.vtkFeatureEdges()
cutBoundary.SetInput(cutPoly)
cutBoundary.Update()
# Write cut line to .vtk file
polyWriter = vtk.vtkPolyDataWriter()
polyWriter.SetInput(cutBoundary.GetOutput())
polyWriter.SetFileName("test.vtk")
polyWriter.Write()
I don't get an error message, the code just hangs at the cutEdges.SetInputConnection(...). The only significant difference as far as I can tell is that I'm using an .stl file instead of a .g file. Shouldn't the outputport data from vtkSTLReader and vtkBYUReader be the same (a vtkAlgorithm object in this case)?
Though I've been using bits and pieces of vtk code over the past year, I'm still very unexperienced in VTK (and image processing in general). Am I supposed to extract the data from vtkSTLReader in a different way? Any help would be greatly appreciated.
Kind regards,
Daniel
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.vtk.org/pipermail/vtkusers/attachments/20130801/bc7dd233/attachment.htm>
More information about the vtkusers
mailing list