[vtkusers] Intersection of plane with 3D surface using vtkCutter

Fri Aug 2 08:17:36 EDT 2013

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:

        # Circumference plane (used to cut)
        circPlane = vtk.vtkPlane()
        circPlane.SetOrigin(G - 150*self.midSagittalPlaneNormal)
        circPlane.SetNormal(self.circPlaneNormal)

        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(circPlane)
        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()

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:

        circPlane = vtk.vtkPlane()
        circPlane.SetOrigin(G - 150*self.midSagittalPlaneNormal)
        circPlane.SetNormal(circPlaneNormal)

        circPlaneSource = vtk.vtkPlaneSource()
        circPlaneSource.SetOrigin(circPlane.GetOrigin())
        circPlaneSource.SetNormal(circPlane.GetNormal())
        circPlaneSource.SetPoint1(Op -
                150*self.midSagittalPlaneNormal)
        circPlaneSource.SetPoint2(G +
                150*self.midSagittalPlaneNormal)
        circPlaneSource.SetResolution(15,15)
        circPlaneSource.Update()

        self.circPlaneMapper = vtk.vtkPolyDataMapper()
        self.circPlaneMapper.SetInput(circPlaneSource.GetOutput())
        self.circPlaneMapper.SetScalarVisibility(False)

        self.circPlaneActor = vtk.vtkActor()
        self.circPlaneActor.GetProperty().SetColor( 0.0 , 1.0 , 0.0 )
        self.circPlaneActor.GetProperty().SetOpacity( 1 )
        self.circPlaneActor.SetMapper(self.circPlaneMapper)
        self.circPlaneActor.GetProperty().SetRepresentationToWireframe()

        self.renderer.AddActor(self.circPlaneActor)

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?) 

Date: Thu, 1 Aug 2013 12:10:45 -0400
Subject: Re: [vtkusers] Intersection of plane with 3D surface using vtkCutter
From: bill.lorensen at gmail.com
To: lacko.daniel at outlook.com
CC: vtkusers at vtk.org

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.

On Thu, Aug 1, 2013 at 11:46 AM, Daniël Lacko <lacko.daniel at outlook.com> wrote:

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

_______________________________________________

Powered by www.kitware.com

Visit other Kitware open-source projects at http://www.kitware.com/opensource/opensource.html

Please keep messages on-topic and check the VTK FAQ at: http://www.vtk.org/Wiki/VTK_FAQ

Follow this link to subscribe/unsubscribe:

http://www.vtk.org/mailman/listinfo/vtkusers

-- 
Unpaid intern in BillsBasement at noware dot com

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.vtk.org/pipermail/vtkusers/attachments/20130802/e2658adf/attachment.htm>