Wrong connected vertices in a non-triangulated mesh

[sergei9838]

Hi, Marco,
I need to find the vertices in a mesh connected to a given vertex (even more than that: the ones which are at most n edges away). Presumably, it can be done this way:

K=Cube().clean().tomesh()  # the mesh - cube
print(K.vertices[0])  # the source vertex [-0.5, -0.5, -0.5]
idx = K.connected_vertices(0)
print(idx)  # Output: [1, 2, 3, 4, 7, 5]

So the result is wrong: it should be [1,3,4] as seen from

print(K.edges)

I am bothered by the notice printed in red at the last command:

INFO| Executing edge extractor: points are renumbered

but still these are indeed the correct edge-neighbours:

print(K.vertices[[1,3,4]])
# array([[-0.5, -0.5,  0.5],
#       [-0.5,  0.5, -0.5],
#      [ 0.5, -0.5, -0.5]], dtype=float32)

It looks like Mesh.connected_vertices just takes the union of all the vertices of the cells containing the source vertex which would produce the correct result if the mesh is triangulated. But wrong, if it is not, like in the case of a cube.

Strangely, there is no Mesh.connected_edges method...

[sergei9838]

A fix could be:

src = 0
idx = K.connected_vertices(src)
eds = K.edges
neighbours_1 = [i for i in idx if [src, i] in eds or [i, src] in eds]  # Output: [1,3,4]

or simply

src = 0
eds = K.edges
neighbours_1 = [i for i in range(len(eds)) if [src, i] in eds or [i, src] in eds]  # Output: [1,3,4]

provided K.edges does not renumber the vertices

[sergei9838]

I wrote a code that definitely does not renumber the vertices:

def get_edges(msh: Mesh) -> list:
    """
    :param msh: a Mesh
    :return: the list of edges of `msh`
    """
    edges = []
    for cell in msh.cells:
        edges += [tuple(sorted([cell[i], cell[i+1]])) for i in range(-1, len(cell)-1)]
        # conversion to tuple is here for set() below to work (removing duplicates)
    return [[e[0], e[1]] for e in set(edges)]

but the result is the same as K.edges gives. Maybe vtkExtractEdges does not renumber them?

[sergei9838]

If this is for some use to the community, here is my solution:

def adjacency_list(msh: Mesh) -> list:
    """
    Computes the adjacency list for mesh `msh` edge-graph
    :param msh: a mesh
    :return: a list with i-th entry being the set if indices of vertices connected by an edge to i-th vertex
    """
    inc = [set()] * msh.nvertices
    for cell in msh.cells:
        if len(cell) > 1:
            for i in range(len(cell)-1):
                inc[cell[i]] = inc[cell[i]].union({cell[i-1], cell[i+1]})
    return inc


def graph_ball(index, n: int, msh: Mesh) -> set:
    """
    Computes the ball of radius `n` in the mesh' edge-graph metric centred in vertex `index`
    :param index: index
    :param n: radius
    :param msh: Mesh
    :return: the set of indices of the vertices which are at most `n` edges from vertex `index`
    """
    if n < 0:
        return set()
    if n == 0:
        return {index}
    else:
        al = adjacency_list(msh)
        ball = {index}
        i = 0
        while i < n and len(ball) < msh.nvertices:
            for v in ball:
                ball = ball.union(al[v])
            i += 1
        return ball

class Peeling:
    """
    Computes the spheres and balls of radii n for n <= `depth` centred in vertex `self.center`
        in the mesh' edge-graph metric
    
    Properties:
    :mesh: a Mesh
    :center: center vertex
    :depth: depth
    :balls: the list of sets. Its n-th element is the set of indices of the vertices
        reachable from `center` by at most n edges
    :spheres: the list of sets. Its n-th element is the set of indices of the vertices
        for which the minimal number of edges to reach the `center` vertex is n
    :adjacency: the adjacency list for the mesh' edge-graph: a list with i-th entry 
        being the set of vertex indices connected by an edge to i-th vertex
    """

    def __init__(self, msh: Mesh, center: int, depth: int = 1):
        self.center = center
        self.mesh = msh
        self.depth = depth
        self.balls = [{self.center}]
        self.spheres = [{self.center}]
        self.adjacency = adjacency_list(msh)

        n = 0
        last_ball = {self.center}
        while n < depth and len(last_ball) < self.mesh.nvertices:
            for v in self.spheres[-1]:
                last_ball = last_ball.union(self.adjacency[v])
            self.balls.append(last_ball)
            self.spheres.append(self.balls[-1] - self.balls[-2])
            n += 1

    def __repr__(self):
        out = f"center vertex: {self.center}\n balls: {self.balls}\n spheres: {self.spheres}"
        return out

[marcomusy]

Thanks Sergei! I'm adding this to the Mesh class!