marcomusy / vedo

A python module for scientific analysis of 3D data based on VTK and Numpy
https://vedo.embl.es
MIT License
2.05k stars 265 forks source link

Suggestion: cut with plane #738

Open saurabhpre opened 1 year ago

saurabhpre commented 1 year ago

Hello, Thank you for your effort. I like it how easy to use this is. I was working with cut_with_plane and notice that it only returns one mesh and discard the other although vtk has the option to return either. I made some changes to the code so that the user can control which mesh to return. Hopefully this improves on already a great repo. `

def cut_with_plane2(mesh, origin=(0, 0, 0), normal=(1, 0, 0), direction=0):
    """
    Cut the mesh with the plane defined by a point and a normal.
    Parameters
    ----------
    origin : array
        the cutting plane goes through this point
    normal : array
        normal of the cutting plane
    direction: bool
        0,1 which cut to keep
    Example:
        .. code-block:: python
            from vedo import Cube
            cube = Cube().cut_with_plane(normal=(1,1,1))
            cube.bc('pink').show()
    .. hint:: examples/simulations/trail.py
        .. image:: https://vedo.embl.es/images/simulations/trail.gif
    Check out also:
        ``crop()``, ``cut_with_box()``, ``cut_with_cylinder()``, ``cut_with_sphere()``
    """
    s = str(normal)
    if "x" in s:
        normal = (1, 0, 0)
        if "-" in s:
            normal = -np.array(normal)
    elif "y" in s:
        normal = (0, 1, 0)
        if "-" in s:
            normal = -np.array(normal)
    elif "z" in s:
        normal = (0, 0, 1)
        if "-" in s:
            normal = -np.array(normal)
    plane = vtk.vtkPlane()
    plane.SetOrigin(origin)
    plane.SetNormal(normal)

    clipper = vtk.vtkClipPolyData()
    clipper.SetInputData(mesh.polydata(True))  # must be True
    clipper.SetClipFunction(plane)
    clipper.GenerateClippedOutputOn()
    clipper.SetValue(0)
    clipper.Update()

    cpoly = clipper.GetOutput(direction)

    if mesh.GetIsIdentity() or cpoly.GetNumberOfPoints() == 0:
        mesh._update(cpoly)
    else:
        # bring the underlying polydata to where _data is
        M = vtk.vtkMatrix4x4()
        M.DeepCopy(self.GetMatrix())
        M.Invert()
        tr = vtk.vtkTransform()
        tr.SetMatrix(M)
        tf = vtk.vtkTransformPolyDataFilter()
        tf.SetTransform(tr)
        tf.SetInputData(cpoly)
        tf.Update()
        mesh._update(tf.GetOutput())

    return mesh

`

marcomusy commented 1 year ago

Hi @saurabhpre Thanks! In principle the user can decide which one to keep by flipping the signs of normal, I understand that maybe it is something that one does not think of , in the sense that it may not be intuitive enough.. That is why other similar methods have an invert=True keyword .. maybe cut_with_plane() should have the same thing even if it's redundant. You are welcome to make a PR if you wish!

Another interesting option i'm thinking of right now is to return the same (uncut) mesh with an additional boolean flag in mesh.pointdata.