`GC_SAFETY_ASSERT FAILURE` for manifold mesh

[jo-mueller]

Hi @nmwsharp ,

first of all: Great library.

No to the actual problem: I am trying to use the EdgeFlipGeodesicSolver to find a path on a relatively simple mesh (vertices and faces below) and am receiving an GC_SAFETY_ASSERT FAILURE error like this:

Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "c:\Users\my_user\AppData\Local\mambaforge\envs\stress\lib\site-packages\potpourri3d\mesh.py", line 64, in __init__
    self.bound_solver = pp3db.EdgeFlipGeodesicsManager(V, F)
RuntimeError: GC_SAFETY_ASSERT FAILURE from D:\a\potpourri3d\potpourri3d\deps\geometry-central\src\surface\manifold_surface_mesh.cpp:110 - duplicate edge in list 1 -- 2

Now, I saw in #18 that this is most likely due to a non-manifold edge. So I did some cleaning (I am using vedo for most mesh operations) with mesh.clean() and also sorted and searched all faces for duplicate triangles. When I check the mesh for whether it is a manifold-mesh using vedo like this, I get:

>>> mesh.is_manifold()
True

Here are the vertices and faces, if it helps to trace down the issue. Here is the mesh in question, in case it helps to debug the issue:

test.zip

Best & thanks for any help, Johannes

[ElisabethBrockhausQC]

Hi @jo-mueller,

I've had the same problem in the MeshVectorHeatSolver. I ended up manually reordering the mesh. It's just a quick a dirty solution, so not very efficient or anything, but maybe it helps you as well. reorder_edges_of_faces() does the trick for me:

Python code

from sys import setrecursionlimit
from typing import Annotated, Literal, TypeVar

import numpy as np
import numpy.typing as npt

DType = TypeVar("DType", bound=np.generic)
ArrayNx3 = Annotated[npt.NDArray[DType], Literal["N", "3"]]


def _get_duplicate_edges(F: ArrayNx3) -> list[list[int]]:
  """Get list of directed edges that belong to more than one face.

  Parameters
  ----------
  F : ArrayNx3
      Array of vertex indices that form the faces of the tessellation.

  Returns
  -------
  list
      List of duplicate edges.
  """
  edges = np.array([[[a, b], [b, c], [c, a]] for a, b, c in F]).reshape(
      F.shape[0] * 3, 2
  )
  e, c = np.unique(edges, axis=0, return_counts=True)
  return e[c > 1].tolist()


def _reorder_neighbours(
  F: ArrayNx3,
  current_face: int,
  faces_to_neighbours_map: dict,
  reordered_faces: list[int],
) -> ArrayNx3:
  """Recursively reorder neighbouring faces until there are no duplicate edges.

  Parameters
  ----------
  F : ArrayNx3
      Array of vertex indices that form the faces of the tessellation.
  current_face : int
      Index of the correctly ordered reference face whose neighbors are reordered.
  faces_to_neighbours_map : dict
      Dictonary mapping face indices to the indices of all not yet ordered neighbours.
  reordered_faces : list[int]
      List of faces that are already in correct order.

  Returns
  -------
  ArrayNx3
      Array of faces where the neighbours of `current_face` are ordered.
  """
  neighbours = faces_to_neighbours_map[current_face]
  for neighbour in neighbours:
      violation = _get_duplicate_edges(F[[current_face] + [neighbour]])
      if violation:
          F[neighbour] = [F[neighbour][2], F[neighbour][1], F[neighbour][0]]
      reordered_faces.append(neighbour)
      faces_to_neighbours_map[current_face].remove(neighbour)
      faces_to_neighbours_map[neighbour].remove(current_face)

  for i in reordered_faces:
      if faces_to_neighbours_map[i]:
          F = _reorder_neighbours(F, i, faces_to_neighbours_map, reordered_faces)
          break

  if not any(faces_to_neighbours_map.values()):
      return F
  else:
      raise RecursionError("Could not reach all vertices.")


def reorder_edges_of_faces(F: ArrayNx3) -> ArrayNx3:
  """Reorder edges of faces for `potpourri3d.MeshVectorHeatSolver`.

  Parameters
  ----------
  F : ArrayNx3
      Array of unordered faces

  Returns
  -------
  ArrayNx3
      Array of ordered faces
  """
  setrecursionlimit(5 * F.shape[0])

  empty_list: list[int] = []

  edges = np.array([[[a, b], [b, c], [c, a]] for a, b, c in F]).reshape(
      F.shape[0] * 3, 2
  )
  edges_to_faces_map = {(min(e), max(e)): empty_list.copy() for e in edges}
  for i, f in enumerate(F):
      edges_to_faces_map[(min(f[0], f[1]), max(f[0], f[1]))].append(i)
      edges_to_faces_map[(min(f[1], f[2]), max(f[1], f[2]))].append(i)
      edges_to_faces_map[(min(f[2], f[0]), max(f[2], f[0]))].append(i)

  faces_to_neighbours_map = {i: empty_list.copy() for i in range(F.shape[0])}
  for neighbours in edges_to_faces_map.values():
      if len(neighbours) > 1:
          faces_to_neighbours_map[neighbours[0]].append(neighbours[1])
          faces_to_neighbours_map[neighbours[1]].append(neighbours[0])
      if len(neighbours) > 2:
          faces_to_neighbours_map[neighbours[0]].append(neighbours[2])
          faces_to_neighbours_map[neighbours[2]].append(neighbours[0])
          faces_to_neighbours_map[neighbours[1]].append(neighbours[2])
          faces_to_neighbours_map[neighbours[2]].append(neighbours[1])
      if len(neighbours) > 3:
          raise ValueError("Something weird is happening.")

  return _reorder_neighbours(F, 0, faces_to_neighbours_map, [0])

[jo-mueller]

Hey there @ElisabethBrockhausQC ,

sorry for the late reply. That absolutely solves the problem, although I must have admit I don't understand why. Is there a hidden convention for how faces are supposed to be sorted for this library to work?

PS: Ok, it only solves my problem in some cases, in others I get a recursion error or still the same GC_SAFETY_ASSERT FAILURE error :/

[ElisabethBrockhausQC]

Is there a hidden convention for how faces are supposed to be sorted for this library to work?

It has been a while since I've worked with this, but IIRC the problem arises if the same edge is part of two faces and the corresponding to vertices are in the same order in the definitions of both faces.

it only solves my problem in some cases, in others I get a recursion error

The RecursionError is explicitly raised in my code snippet:

  if not any(faces_to_neighbours_map.values()):
      return F
  else:
      raise RecursionError("Could not reach all vertices.")

You can try whether increasing the recursion limit here solves the problem.

or still the same GC_SAFETY_ASSERT FAILURE error

😢 I think I cannot help with this, sorry.

[jo-mueller]

Hi @ElisabethBrockhausQC , thanks for coming back to this so quickly. I found the error myself. It's not enough for the mesh to be non-manifold, but it also needs to be oriented, so that all normals either face out or in. For simple geometries like the ones I am working with that is fairy easy to achieve but I guess it could be a bit harder for more intricate structures.

I'll close this, hopefully it helps others that stumble upon the same problem.