Project points/surface to a plane

[S0Phon]

I have an stl file and would like to get the indices of all the cells visible from a chosen viewpoint. I know this is possible in vtk using the vtkHardwareSelector but am curious if it is possible in pyvista?

In addition I would also like to know if there is an existing fuction for projecting 3D point data onto a plane?

[banesullivan]

Hi @S0Phon, these are interesting questions!

an stl file and would like to get the indices of all the cells visible from a chosen viewpoint. I know this is possible in vtk using the vtkHardwareSelector but am curious if it is possible in pyvista?

I'll make a new feature request out of this in the main PyVista repo and keep you posted.

I would also like to know if there is an existing fuction for projecting 3D point data onto a plane?

This is something I've done before actually - let me dig through some old code and try to make an example for you

[banesullivan]

Here's an example of projecting the points of PolyData to a plane which I will implement as a new filter for PolyData object in pyvista/pyvista#276:

import numpy as np
import pyvista as pv
from pyvista import examples

def make_example_data():
    surface = examples.download_saddle_surface()
    points = examples.download_sparse_points()
    poly = surface.interpolate(points, radius=12.0)
    return poly

poly = make_example_data()
poly.plot()

Now let's project that data to a plane defined by some origin and some normal:

def project_points_to_plane(mesh, origin=None, normal=(0,0,1), inplace=False):
    """Project points of this mesh to a plane"""
    if not isinstance(mesh, (pv.PolyData)):
        raise TypeError('Please use surface meshes only.')
    import vtk
    if origin is None:
        origin = mesh.center
    if not inplace:
        mesh = mesh.copy()
    # Make plane
    normal = normal / np.linalg.norm(normal) # MUST HAVE MAGNITUDE OF 1
    plane = vtk.vtkPlane()
    plane.SetOrigin(origin)
    plane.SetNormal(normal)
    # Perform projection in place on the copied mesh
    f = lambda p: plane.ProjectPoint(p, p)
    np.apply_along_axis(f, 1, mesh.points)
    if not inplace:
        return mesh
    return

# Project that surface to a plane
og = poly.center
og[-1] -= poly.length / 3.
projected = project_points_to_plane(poly, origin=og)

p = pv.Plotter()
p.add_mesh(poly,)
p.add_mesh(projected, )
p.show()

[S0Phon]

Thank you very much for your help @banesullivan! For reorganising my questions into the appropriate repositories and for the above method. I have tried the code you provided above using the example data and then again with data from an stl file and it works very well.

[akaszynski]

We probably want to save all of these as examples on the docs.

[S0Phon]

I have reopened this as I noticed an issue with the code when the origin is not defined and the plane normal vector is changed significantly. In this case I used a normal vector of normal=(-1,0,1) (this was the only thing I changed).

As you can see the new projected plane just interects the old one. This may not be an issue in some cases (and is not for my specific case) but if needing both the original data and the projected data in the same image then it is an issue.
Below is an example using a test stl file for the same normal vector as above where the intersection is more apparent.

The zip folder containing the stl file I used for the image above.
double_sphere.zip

[banesullivan]

@S0Phon - this seems like an issue with the default parameter for how a plane is chosen. Do you have any suggestions for how to choose a default plane location (the origin parameter) to make this a bit better?

Right now, the function just uses the center of the input mesh for the plane location... so it will always intersect:

if origin is None:
    origin = mesh.center

Maybe we should do something more like:

if origin is None:
    origin = mesh.center
    origin[-1] -= mesh.length

but even that is not guaranteed to never overlap

[banesullivan]

@akaszynski

We probably want to save all of these as examples on the docs.

I included this example in pyvista/pyvista#276 so it'll be in the gallery!

Also, my intention with using a GitHub repo for a support forum like this was that all of these discussions/examples will be here permanently so that we can search through them and add good examples to the docs when we have time in the future without trying to dig through emails/slack discussions

[S0Phon]

I am not very good at writing neat code or pythonic code or whatever but I think something along these lines should work:

if origin is None:
        original_origin = (mesh.center)
        origin = [original_origin[0] - normal[0]*mesh.length/2, original_origin[1] - normal[1]*mesh.length/2, original_origin[2] - normal[2]*mesh.length/2]

What I am doing here is creating a new point to be the origin where the new point is a distance of mesh.length/2 along the normal vector.
If I have understood the length of mesh.length correctly then this will have a point intersect only where an object occupies the boundary box such as for a cuboid. To avoid this the distance along the normal axis could be mesh.length*0.6.
I am sure that this could be written much better than I have done here. Perhaps using numpy arrays and writing as:

origin = np.array(mesh.center)-np.array(normal)*(mesh.length/2)

or if using the slightly longer distance mesh.length*0.6:

origin = np.array(mesh.center)-np.array(normal)*(mesh.length*0.6)

This may benefit from some rounding as in the example I tried I had 19 significant figures for the new origin.

It is a little unsatisfying as the distance of the projection from the edge of the object is not the same for all possible normal vectors but I am not sure how to get the dimensions of the boundary box so cannot suggest a fix for this at the moment.

Possible Issue
I had python crash on me fairly regularly while testing this. Perhaps the vtk renderer is remaining open and then causing python to crash on a subsequent run? (This is just a guess from issues I have read about when using vtk and could be totally wrong). I am not sure what information is needed to try to discover the problem so will wait until someone responds and says what is needed.

[banesullivan]

Thanks for the suggestion, @S0Phon! I went ahead and used 2. as the factor and implemented it:

if origin is None:
    origin = np.array(self.center) - np.array(normal)*self.length/2.

It is a little unsatisfying as the distance of the projection from the edge of the object is not the same for all possible normal vectors but I am not sure how to get the dimensions of the boundary box so cannot suggest a fix for this at the moment.

I agree... to get the dimensions of the bounding box, you could use the .bounds property:

>>> mesh.bounds
[-20.012727737426758, 20.0, -0.647991418838501, 40.23965072631836, -0.6093360781669617, 15.127287864685059]

And you could get the ranges (distance along each axis of the bounding box) by:

>>> np.array(mesh.bounds).reshape(-1, 2).ptp(axis=1)
[40.01272774, 40.88764215, 15.73662394]

If you have a better solution using that, let me know!

---

On your issue: what you describe is very strange for this example... I haven't experienced this.

I am not sure what information is needed to try to discover the problem so will wait until someone responds and says what is needed.

Could you install a new package I've made called scooby which will be a dependency of PyVista on our next release for environment reporting during bug reports?

From command line:

pip install --upgrade scooby

then copy paste the output of:

import scooby
core = ['pyvista', 'vtk', 'numpy', 'imageio', 'appdirs', 'scooby']
optional = ['matplotlib', 'PyQt5', 'IPython', 'ipywidgets', 'colorcet',
                'cmocean']
# This will output a report in an interactive Python shell
scooby.Report(core=core, optional=optional)

If that doesn't show anything/ if you are using script files, use the print function to show the report:

print(scooby.Report(core=core, optional=optional))

[S0Phon]

Output is as follows

  Date: Mon Jul 01 14:31:53 2019 GMT Daylight Time

           Windows : OS
                 8 : CPU(s)
             AMD64 : Machine
             64bit : Architecture
           15.8 GB : RAM
            Python : Environment

  Python 3.7.3 (default, Mar 27 2019, 17:13:21) [MSC v.1915 64 bit (AMD64)]

            0.20.4 : pyvista
             8.2.0 : vtk
            1.16.2 : numpy
             2.5.0 : imageio
             1.4.3 : appdirs
             0.3.0 : scooby
             3.1.0 : matplotlib
             5.9.2 : PyQt5
             7.4.0 : IPython
             7.4.2 : ipywidgets

  Intel(R) Math Kernel Library Version 2019.0.3 Product Build 20190125 for
  Intel(R) 64 architecture applications

I hope this is what is needed as I had to change one of the lines as said in the note below. I ran this at the beginning of the script after the module imports.

Note: I think scooby.investigte should be replaced with scooby.Report as it threw an attribute error saying that module 'scooby' has no attribute 'investigate' and when checking the issues section for scooby I found one that said that investigate has been deleted and to use Report.

A more full description of the problem is that roughly every second or third time I run the script, I get python crashing when the second window opens to show the projected image. When it crashes the window opens black and stays that way for several seconds before the windows says not responding.
It has not crashed on opening the first window which shows the imported geometry.

[banesullivan]

Thank you for sharing all the above info, @S0Phon!

To make sure I try reproducing this corectly, can you outline your steps exactly? I'm assuming they are:

1. Place the following in a .py file:

import numpy as np
import pyvista as pv
from pyvista import examples

def make_example_data():
    surface = examples.download_saddle_surface()
    points = examples.download_sparse_points()
    poly = surface.interpolate(points, radius=12.0)
    return poly

poly = make_example_data()
poly.plot()

def project_points_to_plane(mesh, origin=None, normal=(0,0,1), inplace=False):
    """Project points of this mesh to a plane"""
    if not isinstance(mesh, (pv.PolyData)):
        raise TypeError('Please use surface meshes only.')
    import vtk
    if origin is None:
        origin = np.array(mesh.center) - np.array(normal)*mesh.length/2.
    if not inplace:
        mesh = mesh.copy()
    # Make plane
    normal = normal / np.linalg.norm(normal) # MUST HAVE MAGNITUDE OF 1
    plane = vtk.vtkPlane()
    plane.SetOrigin(origin)
    plane.SetNormal(normal)
    # Perform projection in place on the copied mesh
    f = lambda p: plane.ProjectPoint(p, p)
    np.apply_along_axis(f, 1, mesh.points)
    if not inplace:
        return mesh
    return

# Project that surface to a plane
projected = project_points_to_plane(poly)

p = pv.Plotter()
p.add_mesh(poly,)
p.add_mesh(projected, )
p.show()

2. Run the file from the command line:

$ python my_file.py

3. Close the renderers that appear and repeat 2.

Or are you doing something different?

[S0Phon]

I have been doing something different. I copied the above into visual studio code, saved it as a python file and ran it using an anaconda distribution of python using the debug feature of VS Code. This is the only way I have been using python so far.

When I tried the above using the anaconda prompt, I get an issue 'No module named 'pyvista'. I cannot understand this as I have checked in the environment library list for it and it shows up and when I try to reintall it, it says that I already have it. I then made sure was in the correct enviroment and still the same issue.

To further confuse matters, when I tried using a seperate installation of python and run the file I get the error 'NameError: name 'self' not defined.

I appreciate your patience and help with this. I am quite new to programming and have yet to figure out a lot of the management of installations and libraries.

[banesullivan]

Oh whoops the self name error is my bad! I didn't actually proof that script I posted above. I just edited and I think it should be good now.

[S0Phon]

When trying with the self replaced with mesh , it has worked 10 times in both anaconda prompt and in command prompt with no failures. I can't see why that would fix the issue with anaconda but I did restarted my computer so maybe that was it,

Visual Studio Code still has issues but as it seems to be running properly from terminal I think that is a seperate issue.

Thanks you for all your help!

[banesullivan]

Strange indeed - if it is working in anaconda prompt and command prompt then all seems good.

I have a feeling that Visual Studio might not be cleaning up the session between runs (which to me seems like a VS Python issue in general), but I've never used VS before so I can't hypothesize further.

I'm going to close this issue as we implemented this filter in the last release of PyVista 🚀

[Paulschneider007]

then, if I have two surfaces, how could I project the points in the first surface to the second one, and display the color map?

[bluetyson]

Different to doingbthe above twice?

[Paulschneider007]

I find a example in pyvista to calculate the surface distance.