How to extract common points between two surface meshes

[scarpma]

Hello everyone, thanks for the beautiful ecosystem. PyVista is so tremendously useful !!

My question is how can I obtain the indexes of points that are common between two pyvista.PolyData objects.

I have a surface that is split into two parts. Every time that I want the whole surface I just do surface1.merge(surface2). In this way the common points (i.e. points that are along the cut which splits the surface) are just merged because are recognized to be in the same exact location.

Now my problem is that I want to know the indexes of these points that will be merged (indexes relative to one specific pyvista.PolyData object). I would like to select those points by having their indexes idx and by doing surface1.points[idx]. This way I can build an array with specific values along the cut, for example by doing

example_array = np.zeros(surface1.n_points)
example_array[idx] = np.ones(idx.shape[0])
surface1['exampleArray'] = example_array

Do you know how to do this ?

Thank you very much in advance

[scarpma]

I found a possible solution.

basically I first extract the boundary edges of both surfaces. Then I compare these edges (each surface can have multiple separate boundaries) until I find a pair that is equal. Once the common edge is found, one can obtain the original indexes through the 'vtkOriginalPointIds'.

I wrote a little function as an example:

def getCommonPointsIdx(s1, s2):
    '''
    returns the indeces of points inside s1 (surface one)
    which are in common with s2 (surface two).
    '''
    e1 = s1.extract_feature_edges(
        manifold_edges=False,
        boundary_edges=True,
        non_manifold_edges=False,
        feature_edges=False,)

    e1 = list(e1.split_bodies())
    e1c = np.array([a.points.mean(axis=0) for a in e1])

    e2 = s2.extract_feature_edges(
        manifold_edges=False,
        boundary_edges=True,
        non_manifold_edges=False,
        feature_edges=False,)

    e2 = list(e2.split_bodies())
    e2c = np.array([a.points.mean(axis=0) for a in e2])

    for ii, i in enumerate(e1c):
        for jj, j in enumerate(e2c):
            diff = ((i-j)**2).sum()
            if diff < 1e-7:
                edge = e1[ii]
                break

    return edge['vtkOriginalPointIds']

I do the comparison by computing the cartesian coordinates of the center of each boundary.