<html>
<head>
<style><!--
.hmmessage P
{
margin:0px;
padding:0px
}
body.hmmessage
{
font-size: 12pt;
font-family:Calibri
}
--></style></head>
<body class='hmmessage'><div dir='ltr'>That figures. I was so focused on the problem being the input surface that I didn't bother checking the type and output of vtkPlaneSource. Thanks for the swift reply, changing the code to the following worked perfectly:<br><br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; # Circumference plane (used to cut)<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; circPlane = vtk.vtkPlane()<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; circPlane.SetOrigin(G - 150*self.midSagittalPlaneNormal)<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; circPlane.SetNormal(self.circPlaneNormal)<br><br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; surface = vtk.vtkSTLReader()<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; surface.SetFileName(fileName)<br><br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; # Cut a slice out of the surface and convert it to a PolyData object<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; cutEdges = vtk.vtkCutter()<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; cutEdges.SetInputConnection(surface.GetOutputPort())<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; cutEdges.SetCutFunction(circPlane)<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; cutEdges.GenerateCutScalarsOn()<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; cutEdges.SetValue(0, 0.5)<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; cutStrips = vtk.vtkStripper()<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; cutStrips.SetInputConnection(cutEdges.GetOutputPort())<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; cutStrips.Update()<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; cutPoly = vtk.vtkPolyData()<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; cutPoly.SetPoints(cutStrips.GetOutput().GetPoints())<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; cutPoly.SetPolys(cutStrips.GetOutput().GetLines())<br><br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; # Extract boundary from cutPoly<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; cutBoundary = vtk.vtkFeatureEdges()<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; cutBoundary.SetInput(cutPoly)<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; cutBoundary.Update()<br><br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; # Write cut line to .vtk file<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; polyWriter = vtk.vtkPolyDataWriter()<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; polyWriter.SetInput(cutBoundary.GetOutput())<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; polyWriter.SetFileName("test.vtk")<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; polyWriter.Write()<br><br>A follow-up question: what if I want to visualize the plane to check whether it is at the right place? I suspect an implicit function can't be visualized as such and I would first need to convert it to a polydata object. What is the best way to do this? If I do:<br><br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; circPlane = vtk.vtkPlane()<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; circPlane.SetOrigin(G - 150*self.midSagittalPlaneNormal)<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; circPlane.SetNormal(circPlaneNormal)<br><br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; circPlaneSource = vtk.vtkPlaneSource()<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; circPlaneSource.SetOrigin(circPlane.GetOrigin())<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; circPlaneSource.SetNormal(circPlane.GetNormal())<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; circPlaneSource.SetPoint1(Op -<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 150*self.midSagittalPlaneNormal)<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; circPlaneSource.SetPoint2(G +<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 150*self.midSagittalPlaneNormal)<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; circPlaneSource.SetResolution(15,15)<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; circPlaneSource.Update()<br><br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; self.circPlaneMapper = vtk.vtkPolyDataMapper()<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; self.circPlaneMapper.SetInput(circPlaneSource.GetOutput())<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; self.circPlaneMapper.SetScalarVisibility(False)<br><br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; self.circPlaneActor = vtk.vtkActor()<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; self.circPlaneActor.GetProperty().SetColor( 0.0 , 1.0 , 0.0 )<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; self.circPlaneActor.GetProperty().SetOpacity( 1 )<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; self.circPlaneActor.SetMapper(self.circPlaneMapper)<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; self.circPlaneActor.GetProperty().SetRepresentationToWireframe()<br><br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; self.renderer.AddActor(self.circPlaneActor)<br><br>I do see a plane, but it's the inclination is wrong (different from the intersection made by the vtkCutter, because when I visualize the extracted boundary on the surface it is at exactly the right place). If I comment the circPlaneSource.SetPoint functions, I don't even get a plane, just a line. (Which is probably to be expected, since the renderer needs to know the boundaries within which the plane should be drawn, right?) <br><br><div><hr id="stopSpelling">Date: Thu, 1 Aug 2013 12:10:45 -0400<br>Subject: Re: [vtkusers] Intersection of plane with 3D surface using vtkCutter<br>From: bill.lorensen@gmail.com<br>To: lacko.daniel@outlook.com<br>CC: vtkusers@vtk.org<br><br><div dir="ltr">You need a vtkPlane, not a vtkPlaneSource to do the cutting. The vtkPlaceSource produces polydata. The vtkPlane is an implicit function. vtkCutter uses implicit functions to cut the data.<div><br></div><div>
<br></div></div><div class="ecxgmail_extra"><br><br><div class="ecxgmail_quote">On Thu, Aug 1, 2013 at 11:46 AM, Daniël Lacko <span dir="ltr">&lt;<a href="mailto:lacko.daniel@outlook.com" target="_blank">lacko.daniel@outlook.com</a>&gt;</span> wrote:<br>
<blockquote class="ecxgmail_quote" style="border-left:1px #ccc solid;padding-left:1ex;">


<div><div dir="ltr">Hi all,<br><br>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.<br>
<br>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:<br>
<br>&nbsp;&nbsp;&nbsp; # Here we are cutting the cow. Cutting creates lines where the cut<br>&nbsp;&nbsp;&nbsp; # function intersects the model. (Clipping removes a portion of the<br>&nbsp;&nbsp;&nbsp; # model but the dimension of the data does not change.)<br>&nbsp;&nbsp;&nbsp; #<br>
&nbsp;&nbsp;&nbsp; # The reason we are cutting is to generate a closed polygon at the<br>&nbsp;&nbsp;&nbsp; # boundary of the clipping process. The cutter generates line<br>&nbsp;&nbsp;&nbsp; # segments, the stripper then puts them together into polylines. We<br>&nbsp;&nbsp;&nbsp; # then pull a trick and define polygons using the closed line<br>
&nbsp;&nbsp;&nbsp; # segments that the stripper created.<br>&nbsp;&nbsp;&nbsp; cutEdges = vtk.vtkCutter()<br>&nbsp;&nbsp;&nbsp; cutEdges.SetInputConnection(cow.GetOutputPort())<br>&nbsp;&nbsp;&nbsp; cutEdges.SetCutFunction(plane)<br>&nbsp;&nbsp;&nbsp; cutEdges.GenerateCutScalarsOn()<br>&nbsp;&nbsp;&nbsp; cutEdges.SetValue(0, 0.5)<br>
&nbsp;&nbsp;&nbsp; cutStrips = vtk.vtkStripper()<br>&nbsp;&nbsp;&nbsp; cutStrips.SetInputConnection(cutEdges.GetOutputPort())<br>&nbsp;&nbsp;&nbsp; cutStrips.Update()<br>&nbsp;&nbsp;&nbsp; cutPoly = vtk.vtkPolyData()<br>&nbsp;&nbsp;&nbsp; cutPoly.SetPoints(cutStrips.GetOutput().GetPoints())<br>&nbsp;&nbsp;&nbsp; cutPoly.SetPolys(cutStrips.GetOutput().GetLines())<br>
<br>&nbsp;&nbsp;&nbsp; # Extract boundary from cutPoly<br>&nbsp;&nbsp;&nbsp; cutBoundary = vtk.vtkFeatureEdges()<br>&nbsp;&nbsp;&nbsp; cutBoundary.SetInput(cutPoly)<br>&nbsp;&nbsp;&nbsp; cutBoundary.Update()<br>&nbsp;&nbsp;&nbsp; # Write cut line to .vtk file<br>&nbsp;&nbsp;&nbsp; polyWriter = vtk.vtkPolyDataWriter()<br>
&nbsp;&nbsp;&nbsp; polyWriter.SetInput(cutBoundary.GetOutput())<br>&nbsp;&nbsp;&nbsp; polyWriter.SetFileName("cutPoly.vtk")<br>&nbsp;&nbsp;&nbsp; polyWriter.Write()<br><br>The code works as intended. However, if I try using it with my own code, it won't read the surface data correctly:<br>
<br>&nbsp;&nbsp;&nbsp; # Circumference plane (used to cut)<br>&nbsp;&nbsp;&nbsp; circPlaneSource = vtk.vtkPlaneSource()<br>&nbsp;&nbsp;&nbsp; circPlaneSource.SetOrigin(G - <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 150*self.midSagittalPlaneNormal)<br>&nbsp;&nbsp;&nbsp; circPlaneSource.SetPoint1(Op -<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 150*self.midSagittalPlaneNormal)<br>
&nbsp;&nbsp;&nbsp; circPlaneSource.SetPoint2(G +<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 150*self.midSagittalPlaneNormal)<br>&nbsp;&nbsp;&nbsp; circPlaneSource.SetResolution(15, 15)<br>&nbsp;&nbsp;&nbsp; circPlaneSource.Update()<br><br>&nbsp;&nbsp;&nbsp; surface = vtk.vtkSTLReader()<br>&nbsp;&nbsp;&nbsp; surface.SetFileName(fileName)<br>
<br>&nbsp;&nbsp;&nbsp; # Cut a slice out of the surface and convert it to a PolyData object<br>&nbsp;&nbsp;&nbsp; cutEdges = vtk.vtkCutter()<br>&nbsp;&nbsp;&nbsp; cutEdges.SetInputConnection(surface.GetOutputPort())<br>&nbsp;&nbsp;&nbsp; cutEdges.SetCutFunction(circPlaneSource.GetOutput())<br>
&nbsp;&nbsp;&nbsp; cutEdges.GenerateCutScalarsOn()<br>&nbsp;&nbsp;&nbsp; cutEdges.SetValue(0, 0.5)<br>&nbsp;&nbsp;&nbsp; cutStrips = vtk.vtkStripper()<br>&nbsp;&nbsp;&nbsp; cutStrips.SetInputConnection(cutEdges.GetOutputPort())<br>&nbsp;&nbsp;&nbsp; cutStrips.Update()<br>&nbsp;&nbsp;&nbsp; cutPoly = vtk.vtkPolyData()<br>
&nbsp;&nbsp;&nbsp; cutPoly.SetPoints(cutStrips.GetOutput().GetPoints())<br>&nbsp;&nbsp;&nbsp; cutPoly.SetPolys(cutStrips.GetOutput().GetLines())<br>&nbsp;&nbsp;&nbsp; <br>&nbsp;&nbsp;&nbsp; # Extract boundary from cutPoly<br>&nbsp;&nbsp;&nbsp; cutBoundary = vtk.vtkFeatureEdges()<br>&nbsp;&nbsp;&nbsp; cutBoundary.SetInput(cutPoly)<br>
&nbsp;&nbsp;&nbsp; cutBoundary.Update()<br><br>&nbsp;&nbsp;&nbsp; # Write cut line to .vtk file<br>&nbsp;&nbsp;&nbsp; polyWriter = vtk.vtkPolyDataWriter()<br>&nbsp;&nbsp;&nbsp; polyWriter.SetInput(cutBoundary.GetOutput())<br>&nbsp;&nbsp;&nbsp; polyWriter.SetFileName("test.vtk")<br>&nbsp;&nbsp;&nbsp; polyWriter.Write()<br>
<br>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)? <br>
<br>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.<br>
<br>Kind regards,<br><br>Daniel<br>                                               </div></div>
<br>_______________________________________________<br>
Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br>
<br>
Visit other Kitware open-source projects at <a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br>
<br>
Please keep messages on-topic and check the VTK FAQ at: <a href="http://www.vtk.org/Wiki/VTK_FAQ" target="_blank">http://www.vtk.org/Wiki/VTK_FAQ</a><br>
<br>
Follow this link to subscribe/unsubscribe:<br>
<a href="http://www.vtk.org/mailman/listinfo/vtkusers" target="_blank">http://www.vtk.org/mailman/listinfo/vtkusers</a><br>
<br></blockquote></div><br><br clear="all"><div><br></div>-- <br>Unpaid intern in BillsBasement at noware dot com<br>
</div></div>                                               </div></body>
</html>