[feat] update mesheuler, add test to v2s grayscale, finddisconnsurf

Author: fangq

iso2mesh/core.py

@@ -325,7 +325,7 @@ def vol2surf(img, ix, iy, iz, opt, dofix=0, method="cgalsurf", isovalues=None):
 
                 keepratio = (
                     opt.get("keepratio", 1)
-                    if len(opt) == 1
+                    if isinstance(opt, dict)
                     else opt[i].get("keepratio", 1)
                 )
                 print(f"Resampling surface mesh for level {i + 1}...")
@@ -337,21 +337,23 @@ def vol2surf(img, ix, iy, iz, opt, dofix=0, method="cgalsurf", isovalues=None):
             else:
                 radbound = (
                     opt.get("radbound", 1)
-                    if len(opt) == 1
+                    if isinstance(opt, dict)
                     else opt[i].get("radbound", 1)
                 )
                 distbound = (
                     opt.get("distbound", radbound)
-                    if len(opt) == 1
+                    if isinstance(opt, dict)
                     else opt[i].get("distbound", radbound)
                 )
                 maxsurfnode = (
                     opt.get("maxnode", 40000)
-                    if len(opt) == 1
+                    if isinstance(opt, dict)
                     else opt[i].get("maxnode", 40000)
                 )
                 surfside = (
-                    opt.get("side", "") if len(opt) == 1 else opt[i].get("side", "")
+                    opt.get("side", "")
+                    if isinstance(opt, dict)
+                    else opt[i].get("side", "")
                 )
 
                 if surfside == "upper":

iso2mesh/trait.py

@@ -193,6 +193,9 @@ def extractloops(edges):
     # Remove degenerate edges (edges where the start and end points are the same)
     edges = edges[edges[:, 0] != edges[:, 1], :]
 
+    if len(edges) == 0:
+        return loops
+
     # Initialize the loop with the first edge
     loops.extend(edges[0, :])
     loophead = edges[0, 0]
@@ -751,18 +754,21 @@ def flatsegment(node, edge):
 
 def mesheuler(face):
     """
-    X, V, E, F = mesheuler(face)
+    X, V, E, F, b, g, C = mesheuler(face)
 
     Compute Euler's characteristic of a mesh.
 
     Arguments:
     face : a closed surface mesh (Mx3 array where M is the number of faces and each row contains vertex indices)
 
     Output:
-    X : Euler's characteristic (X = 2 - 2 * g, where g is the genus)
+    X : Euler's characteristic (X = V - E + F - C)
     V : number of vertices
     E : number of edges
-    F : number of faces
+    F : number of triangles (if face is tetrahedral mesh, exterior surface)
+    b : number of boundary loops (for surfaces)
+    g : genus (holes)
+    C : number of tetrahedra
 
     Author: Qianqian Fang
     This function is part of the iso2mesh toolbox (http://iso2mesh.sf.net)
@@ -771,16 +777,30 @@ def mesheuler(face):
     # Number of vertices
     V = len(np.unique(face))
 
-    # Construct edges from faces
+    # Number of unique edges
     E = uniqedges(face)[0].shape[0]
 
-    # Number of faces
-    F = face.shape[0]
+    b = 0  # open-boundary loops
+    g = 0  # genus
+    C = 0  # tet cells
+
+    # Number of unique faces
+    if face.shape[1] == 4:
+        F = uniqfaces(face)[0].shape[0]
+        C = face.shape[0]
+    else:
+        ed = surfedge(face)[0]
+        loops = extractloops(ed)
+        b = np.sum(np.isnan(loops))
+        F = face.shape[0]
+
+    # Euler's formula, X = V - E + F - C - 2*g
+    X = V - E + F - C
 
-    # Euler's characteristic formula: X = V - E + F
-    X = V - E + F
+    if face.shape[1] == 3:
+        g = (X + b - 2) // 2
 
-    return X, V, E, F
+    return X, V, E, F, b, g, C
 
 
 def orderloopedge(edge):

test/run_test.py

@@ -441,10 +441,13 @@ def test_meshconn(self):
 
     def test_mesheuler(self):
         eu = mesheuler(self.el)
-        self.assertEqual(eu, (-9, 12, 33, 12))
+        self.assertEqual(eu, (1, 12, 33, 34, 0, 0, 12))
 
         eu = mesheuler(self.fc)
-        self.assertEqual(eu, (2, 12, 30, 20))
+        self.assertEqual(eu, (2, 12, 30, 20, 0, 0, 0))
+
+        eu = mesheuler(self.fc[1:-2, :])
+        self.assertEqual(eu, (0, 12, 29, 17, 2, 0, 0))
 
     def test_neighborelem(self):
         expected = [
@@ -936,8 +939,9 @@ def __init__(self, *args, **kwargs):
         xi, yi, zi = np.meshgrid(
             np.arange(0, 61), np.arange(0, 58), np.arange(0, 55), indexing="ij"
         )
-        dist = (xi - 30) ** 2 + (yi - 30) ** 2 + (zi - 30) ** 2
-        self.mask = dist < 400
+        self.dist = (xi - 30) ** 2 + (yi - 30) ** 2 + (zi - 30) ** 2
+        self.dist = self.dist.astype(float)
+        self.mask = self.dist < 400
         self.mask = self.mask.astype(float)
         self.mask[25:35, 25:35, :] = 2
 
@@ -1014,11 +1018,22 @@ def test_cgalv2m(self):
         self.assertAlmostEqual(sum(elemvolume(no[:, :3], el[:, :4])), 0.7455, 3)
 
     def test_v2s_label(self):
-        no, fc, _, _ = v2s(self.mask, 0.04, 1)
+        no, fc, _, _ = v2s(self.mask, 0.5, 1)
         self.assertAlmostEqual(
-            sum(elemvolume(no[:, :3], fc[:, :3])) * 0.001, 6.356848339883229, 3
+            sum(elemvolume(no[:, :3], fc[:, :3])) * 0.001, 5.802998130608866, 2
         )
 
+    def test_v2s_grayscale(self):
+        no, fc, _, _ = v2s(self.dist, [200, 400], 1)
+        self.assertAlmostEqual(
+            sum(elemvolume(no[:, :3], fc[:, :3])) * 0.0001, 4.816236473144449, 3
+        )
+
+    def test_finddisconnsurf(self):
+        no, fc, _, _ = v2s(self.dist, [90, 200, 400], 5)
+        fcs = finddisconnsurf(fc[:, :3])
+        self.assertEqual(len(fcs), 5)
+
     def test_v2m_cgalmesh(self):
         no, el, _ = v2m(
             self.mask.astype(np.uint8),
@@ -1036,7 +1051,7 @@ def test_v2m_simplify(self):
             self.mask.astype(np.uint8), 1.5, {"keepratio": 0.5}, 10, "simplify"
         )
         self.assertAlmostEqual(
-            sum(elemvolume(no[:, :3], el[:, :4])), 5456.390624999979, 0
+            sum(elemvolume(no[:, :3], el[:, :4])) * 0.001, 5.456390624999979, 2
         )