How to add edges to join two detached meshes?

[JeffreyWardman]

I have a a = mesh.boundaries() and a a_smooth = mesh.boundaries().smooth_mls_1d().shift(dz=-1). I want to be able to join them together via new edges determined by nearest neighbours.

My approach was to calculate the nearest neighbour of a in a_smooth and to then append them to a copy of a. I have calculated the Euclidean distance between the points. I thought at this stage I could append the lines and vertices and it would work as intended, however, it only shows the original boundary points. When I run .clean() it removes the smoothed boundary.

# Add points to mesh to join it with smoothed border, that the boundaries can join to after extruding downwards.
vertices = boundaries.vertices.copy()
smoothed_vertices = boundaries_smoothed.vertices.copy()
nearest_neighbour = torch.cdist(torch.as_tensor(vertices), torch.as_tensor(smoothed_vertices)).topk(
    k=1, largest=False, sorted=True
)

merged = boundaries.copy()
merged.vertices = np.vstack([merged.vertices, boundaries_smoothed.vertices.copy()])
for line in boundaries_smoothed.lines:
    merged.lines.append(line)

# Smoothed vertices indices will be appended to non-smoothed so will start after them.
initial_val = len(boundaries.vertices.copy())  # No need to index + 1
for i, neighbour in enumerate(nearest_neighbour.indices.squeeze(dim=-1).numpy()):
    merged.lines.append([i, initial_val + neighbour])

[marcomusy]

You maybe forgot to create a copy?

from vedo import *
mesh = Mesh(dataurl+'bunny.obj').cut_with_plane([0.02,0,0]).c("blue7")
boundaries = mesh.boundaries()
boundaries_smoothed = boundaries.clone().smooth_mls_1d(f=2.0)
lines = Lines(boundaries, boundaries_smoothed).lw(2)
show(mesh, boundaries, boundaries_smoothed, lines, axes=1)

[JeffreyWardman]

Ah apologies! I forgot to add in the shifting downward into the larger code section. It should be boundaries_smoothed = boundaries.copy().smooth_mls_1d(f=0.2).shift(dz=-1), not boundaries_smoothed = boundaries.copy().smooth_mls_1d(f=0.2). So the two boundaries are therefore completely separate from each other and without overlap.

[marcomusy]

not sure if I understand what your're trying to do... extruding the boundary?

from vedo import *

mesh = Mesh(dataurl+'bunny.obj').cut_with_plane([0.02,0,0]).c("blue7")

boundaries = mesh.boundaries().join(reset=True)
boundaries_smoothed = boundaries.copy().smooth_mls_1d(f=2.0).shift(dx=-.1)

ribbon = Ribbon(boundaries, boundaries_smoothed).c("red7").alpha(0.5)

show(mesh, boundaries, boundaries_smoothed, ribbon, axes=1)

[JeffreyWardman]

boundaries.zip

I'm trying to connect these two boundaries. Using a ribbon doesn't seem to work unfortunately.

This also relates to #968. The last step in replicating the madcad algorithm is to join the inflated and deflated meshes with mkquad (https://pymadcad.readthedocs.io/en/latest/_modules/madcad/mesh/mesh.html)

[marcomusy]

It doesn't work because the points are aligned but not ordered:

from vedo import *

b0 = Mesh("boundaries.vtk")
b1 = Mesh("boundaries_smoothed.vtk")

show(
    b0, b1,
    Lines(b0,b1),
    b0.labels("id", scale=0.05),
    b1.labels("id", scale=0.05),
    axes=1,
).close()

# THIS DOES NOT WORK YET BUT IT'S THE WAY TO GO...:
b0.add_ids()
b1.add_ids()

ids0 = b0.pointdata["PointID"]
ids1 = b1.pointdata["PointID"]

points = b0.coordinates.tolist() + b1.coordinates.tolist()
ids = ids0.tolist() + ids1.tolist()

n = b1.npoints
faces = []
for k in range(b0.npoints-1):
    i0 = ids[k]
    i1 = ids[k+1]
    j0 = ids[k+n]+n
    faces.append([i0, i1, j0])
# print(faces)
m = Mesh([points, faces], c="k6")
show(m, axes=1).close()

[marcomusy]

...not sure if this can help:

from vedo import *
pts = Points("data/boundaries.vtk").c("red5").ps(5)
line = pts.clone().subsample(0.005).generate_segments().clean().join()
print(line)
show(pts, line, axes=1)

[JeffreyWardman]

I can't figure out a method to reorder the points unfortunately. Do you have a suggestion? Unfortunately I want to ensure the original boundary is unaffected so your last comment won't be appropriate.

I can't set lines (apart from deleting a value one by one with del) but maybe swapping lines with this:

import networkx as nx

lines = boundaries.lines.copy()
G = nx.DiGraph()
G.add_edges_from(lines)
traversal_path = list(nx.eulerian_circuit(G))

[JeffreyWardman]

@marcomusy how can I overwrite the existing boundary.lines with new indices? I should be able to resolve it after doing this (or creating a new Mesh with just vertices and lines).

[marcomusy]

Have you tried Mesh([vertices, [], lines]) ?

[JeffreyWardman]

Got it to work! Thanks again for your help.

I reordered points with:

def reorder_points(boundaries):
    lines = boundaries.lines.copy()
    vertices = boundaries.vertices.copy()

    G = nx.DiGraph()
    G.add_edges_from(lines)
    traversal_path = list(nx.eulerian_circuit(G))
    eulerian_nodes = [edge[0] for edge in traversal_path] + [traversal_path[-1][-1]]

    new_pts = np.array([vertices[i] for i in eulerian_nodes])
    new_boundary = vedo.Mesh([new_pts, [], [[i, i + 1] for i in range(len(lines))]])
    return new_boundary

Lines were connected with two triangles via a slight modification to your recommendation (the second append to faces):

def connect_lines_with_faces(b0: vedo.Mesh, b1: vedo.Mesh) -> vedo.Mesh:
    b0.add_ids()
    b1.add_ids()

    ids0 = b0.pointdata["PointID"]
    ids1 = b1.pointdata["PointID"]

    points = b0.coordinates.tolist() + b1.coordinates.tolist()
    ids = ids0.tolist() + ids1.tolist()

    n = b1.npoints
    faces = []
    for k in range(b0.npoints - 1):
        i0 = ids[k]
        i1 = ids[k + 1]
        j0 = ids[k + n] + n
        faces.append([i0, i1, j0])

        j1 = ids[k + 1 + n] + n
        faces.append([i1, j0, j1])

    m = vedo.Mesh([points, faces], c="k6")
    return m

[marcomusy]

Thanks Jeffrey, that's a clever solution! As vedo supports triangle strips too, your script gave me the idea to try out this:

from vedo import *
from vedo import vtkclasses as vtki

def to_strips(b0, b1, closed=True):

    b0 = b0.clone().join()
    b1 = b1.clone().join()

    vertices0 = b0.vertices.tolist()
    vertices1 = b1.vertices.tolist()
    lines0 = b0.lines
    lines1 = b1.lines
    m =  len(lines0) #TODO generalize to multiple lines in b0 and b1
    assert m == len(lines1)

    for j in range(m):
        n = len(lines0[j])
        assert n == len(lines1[j])

        ids0j = list(lines0[j])
        ids1j = list(lines1[j])

        if closed:
            ids0j.append(ids0j[0])
            ids1j.append(ids1j[0])
            vertices0.append(vertices0[ids0j[0]])
            vertices1.append(vertices1[ids1j[0]])
            n = n + 1

        # create a triangle strip
        # https://vtk.org/doc/nightly/html/classvtkTriangleStrip.html
        tstrip = vtki.new("TriangleStrip")
        tstrip_ids = tstrip.GetPointIds()
        tstrip_ids.SetNumberOfIds(2*n)

        points = vtki.new("Points")
        for ipt in range(n):
            points.InsertNextPoint(vertices0[ids0j[ipt]])
            points.InsertNextPoint(vertices1[ids1j[ipt]])

        for i in range(2*n):
            tstrip_ids.SetId(i, i)

        cells = vtki.new("CellArray")
        cells.InsertNextCell(tstrip)

    polydata = vtki.vtkPolyData()
    polydata.SetPoints(points)
    polydata.SetStrips(cells)

    m = Mesh(polydata, c="k6")
    return m

b0 = Mesh("data/boundaries/boundaries.vtk")
b1 = Mesh("data/boundaries/boundaries_smoothed.vtk").shift(dz=10)

# #TODO multiple lines, (doesnt work right now)
# b0c = b0.clone().shift([90,0,0])
# b1c = b1.clone().shift([90,0,0])
# b0 = merge(b0, b0c)
# b1 = merge(b1, b1c)

m = to_strips(b0, b1).triangulate()
print(m)

show(
    b0, b1,
    # Lines(b0,b1),
    m,
    axes=1,
).close()

the advantage is that we don't need networkx + it would be generalizable to Mesh objs that contain multiple lines.

[marcomusy]

Just pushed a version to add a slot in constructor of Mesh() so now the above simplifies to

from vedo import *

def to_strips(b0, b1, closed=True):

    b0 = b0.clone().join()
    b1 = b1.clone().join()

    vertices0 = b0.vertices.tolist()
    vertices1 = b1.vertices.tolist()

    lines0 = b0.lines
    lines1 = b1.lines
    m =  len(lines0)
    assert m == len(lines1)

    strips = []
    points = []

    for j in range(m):

        ids0j = list(lines0[j])
        ids1j = list(lines1[j])

        n = len(ids0j)
        assert n == len(ids1j)

        if closed:
            ids0j.append(ids0j[0])
            ids1j.append(ids1j[0])
            vertices0.append(vertices0[ids0j[0]])
            vertices1.append(vertices1[ids1j[0]])
            n = n + 1

        strip = []  # create a triangle strip
        npt = len(points)
        for ipt in range(n):
            points.append(vertices0[ids0j[ipt]])
            points.append(vertices1[ids1j[ipt]])

        strip = list(range(npt, npt+2*n))
        strips.append(strip)

    return Mesh([points, [], [], strips], c="k6")

b0 = Mesh("data/boundaries/boundaries.vtk")
b1 = Mesh("data/boundaries/boundaries_smoothed.vtk").shift(dz=10)

b0c = b0.clone().shift([90,0,0])
b1c = b1.clone().shift([90,0,0])
b0 = merge(b0, b0c)
b1 = merge(b1, b1c)

m = to_strips(b0, b1)#.triangulate()
print(m)

show(b0, b1, m, axes=1).close()

...and works for multiple lines too :)

[JeffreyWardman]

Awesome! Nice stuff!