How to access point data for a cell?

[S0Phon]

I would like to know how to extract point data for polydata on a per cell basis. Ideally this would be in a format of a nx3x3 numpy array, where n is the number of triangles in a trianglular mesh and the x,y,z coordinates of each vertex are given in each row of the 3x3 matrix. However just being able to get an array of the points for each cell would be enough as I could combine and reshape this afterwards.

My reason for doing this is to input these points into another program (which requires this format) after using pyvista to manipulate the original STL import.

I feel like this is probably quite an easy thing to do but so far have been unable to figure out how. I appologise if this has been covered somewhere before but I have searched for some time and cannot find examples of getting the triangle point data, only the mesh point data by using mesh.points.

To create a similar situation using pyvista example files, I have found an example @banesullivan posted on a discourse question.

import vtk
import pyvista as pv
from pyvista import examples

# get some poly data object - use PyVista eample with cell data present
poly_data = examples.download_doorman()

# Make some implicit function to clip with:
plane = vtk.vtkPlane()
plane.SetNormal([1,0,0])
plane.SetOrigin(poly_data.GetCenter())

alg = vtk.vtkClipPolyData()
alg.SetInputDataObject(poly_data)
alg.SetClipFunction(plane) # the the cutter to use the plane we made
alg.SetInsideOut(False) # invert the clip if needed
alg.Update() # Perfrom the Cut
clipped = alg.GetOutput()

# Plot it up... use PyVista because it's easy
pv.plot(clipped)

I would then need to extract the point data for each triangle in clipped as described at the start of the issue.

[GuillaumeFavelier]

I assume the variable clipped is a vtkPolyData so you can wrap it with:
clipped = pv.PolyData(clipped)
By default, I think it's only a shallow copy but more details can be found in the pyvista.core.PolyData class.

And then you have access to clipped.faces which is a one-dimensional array of cells following VTK convention [cell0_nverts, cell0_v0, cell0_v1, cell0_v2, cell1_nverts, ...]. You can obtain the number of faces by using the clipped.n_faces attribute.

Is it helpful @S0Phon?

btw, I see in the doc it's written that faces returns a pointer to the points as a numpy object, we will correct this (points for points and faces for cells).

[S0Phon]

Thank you for responding. If I understand correctly then I think this is very helpful and has answered my question.
Just to check, are the cell0_v0, cell0_v1 and cell0_v2 etc, indices to the points in the points array given by clipped.points?

[GuillaumeFavelier]

Just to check, are the cell0_v0, cell0_v1 and cell0_v2 etc, indices to the points in the points array given by clipped.points?

Exactly!

[S0Phon]

Great! Thank you again @GuillaumeFavelier . I shall close this issue

[banesullivan]

FYI, a true PyVista version of this would be:

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

poly_data = examples.download_doorman()

clipped = poly_data.clip(normal=[1,0,0], 
                         origin=poly_data.center, 
                         invert=False)

clipped.plot()

Then, you could grab the cell nodes like so (if I understand this correctly... if not, @S0Phon, would you please post your solution):

cell_indices = clipped.faces
# make sure all triangles
assert not cell_indices.size % 4 and np.all(cell_indices.reshape(-1, 4)[:,0] == 3)
# Grab ids for all cell nodes
cell_node_ids = cell_indices.reshape(-1, 4)[:,1:4].ravel()

cell_nodes = clipped.points[cell_node_ids]
cell_nodes

array([[1.666226, 3.261832, 0.646217],
       [1.658398, 3.293886, 0.599486],
       [1.615146, 3.251326, 0.649149],
       ...,
       [0.11784 , 5.763801, 0.25697 ],
       [0.126947, 5.757092, 0.256631],
       [0.121833, 5.77197 , 0.242532]], dtype=float32)

[S0Phon]

This is great. Thank you @banesullivan, it is much more concise than what I had done.