CODEWITHSUNDEEP
pythonintersectionvtk
Giojoe89
Giojoe89

Reputation: 11

How to compute the intersection area between a plane (vtkCutter object) and a 3D geometry in Python?

I have been (unsuccessfully) trying to use the vtk toolkit in Python (v3.8.8) to extract the intersection area between a (cutting) plane and a given 3D shape (STL file). I get to the point of defining a vtkCutter, perform the clipping (with a plane) and extract the polyline and points (perimeter) of the intersection using a vtkStripper and vtkPolydata (see snapshot below).

def findIntersection(mySTL,Pt1,Pt2,midPt):
    # create cut plane
    lineVec = unitVec(Pt1,Pt2)
    plane = vtk.vtkPlane()
    plane.SetOrigin(midPt[0][0], midPt[0][1], midPt[0][2])
    plane.SetNormal(lineVec[0],lineVec[1],lineVec[2])
    # create cutter
    cutEdges = vtk.vtkCutter()
    cutEdges.SetInputConnection(mySTL)
    cutEdges.SetCutFunction(plane)
    cutEdges.GenerateCutScalarsOn()
    cutEdges.SetValue(0, 0.5)
    cutStrips = vtk.vtkStripper()
    cutStrips.SetInputConnection(cutEdges.GetOutputPort())
    cutStrips.Update()
    # get polydata intersection (plane vs stl)
    cutPoly = vtk.vtkPolyData()
    cutPoly.SetPoints(cutStrips.GetOutput().GetPoints())
    cutPoly.SetPolys(cutStrips.GetOutput().GetLines())

    return cutEdges, cutPoly, cutStrips

I tried various methods to extract the intersection area (all failing but one – reported below - that returns values in the expected ranges). However, I am getting values far from independent measures performed in other software, which I trust.

def getIntersectArea(contour):
    CoM = vtk.vtkCenterOfMass()
    CoM.SetInputData(contour) #contour is the vtk.PolyData (e.g cutPoly)
    CoM.SetUseScalarsAsWeights(False)
    CoM.Update()
    Center = CoM.GetCenter()
    area = float(0)
    points = contour.GetCell(0).GetPoints()
    numberOfPoints = points.GetNumberOfPoints()
    for i in range(numberOfPoints):
        point0 = points.GetPoint(i)
        point1 = points.GetPoint((i+1)%numberOfPoints)
        pArea = vtk.vtkTriangle.TriangleArea(Center,point0,point1)
        area = area+float(pArea)

    return area

I have not been able to find an appropriate solution, as of yet. Does anyone have a solution or can point me to other discussions (which I may have missed)?

Thanks!

PS: I am exploiting the vtk library to perform a bunch of other operations on STL files, all working.

Upvotes: 1

Views: 1288

Answers (1)

mmusy
mmusy

Reputation: 1337

Using vedo:

from vedo import *

msh = Mesh(dataurl+'bunny.obj').scale(3).shift(0,-0.5,0.01)
plane = Grid(resx=100, resy=100).wireframe(False)

cutplane = plane.clone().cutWithMesh(msh).triangulate()
area = cutplane.area()

show([[msh, plane], [cutplane, f"area: {area}"]], N=2, axes=1)

enter image description here

Upvotes: 1

Related Questions