[brain2mesh] accelerate brain2mesh, binsurface avoid meshcheckrepair

Author: fangq

README.md

@@ -4,7 +4,7 @@
 
 * **Copyright**: (C) Qianqian Fang (2024-2025) <q.fang at neu.edu>
 * **License**: GNU Public License V3 or later
-* **Version**: 0.4.0
+* **Version**: 0.4.1
 * **URL**: [https://pypi.org/project/iso2mesh/](https://pypi.org/project/iso2mesh/)
 * **Homepage**: [https://iso2mesh.sf.net](https://iso2mesh.sf.net)
 * **Github**: [https://github.com/NeuroJSON/pyiso2mesh](https://github.com/NeuroJSON/pyiso2mesh)

iso2mesh/__init__.py

@@ -177,7 +177,7 @@
     thinbinvol,
 )
 
-__version__ = "0.4.0"
+__version__ = "0.4.1"
 __all__ = [
     "advancefront",
     "barycentricgrid",

iso2mesh/brain2mesh.py

@@ -673,8 +673,8 @@ def normalizer(x):
 
     # Generate brain_f if requested
     brain_f = np.array([])
-    #    if len(sys._getframe().f_back.f_code.co_names) > 2:  # Check if brain_f is requested
-    #        brain_f = layersurf(brain_el)
+    if cfg.get("face", False):  # Check if brain_f is requested
+        brain_f = layersurf(brain_el)[0]
 
     return brain_n, brain_el, brain_f
 
@@ -1387,7 +1387,7 @@ def brain1020(
     # At this point, initpoints contains {nz, iz, lpa, rpa, cz0}
     # Plot the head mesh
     if showplot:
-        hh = plotmesh(node, face, alpha=0.6, color="wheat", linewidth=0)
+        hh = plotmesh(node, face, alpha=0, color="wheat", linewidth=0.1)
         ax = hh["ax"][0]
         ax.set_xlabel("X")
         ax.set_ylabel("Y")

iso2mesh/core.py

@@ -331,7 +331,7 @@ def vol2surf(img, ix, iy, iz, opt, dofix=0, method="cgalsurf", isovalues=None):
             else:
                 levelmask = newimg >= isovalues[i]
 
-            levelno, levelel = binsurface(levelmask)
+            levelno, levelel = binsurface(levelmask, clean=False)
 
             if levelel.size > 0:
                 if opt.get("autoregion", 0):
@@ -694,7 +694,7 @@ def surf2vol(node, face, xi, yi, zi, **kwargs):
     return img.astype(np.uint8)
 
 
-def binsurface(img, nface=3):
+def binsurface(img, nface=3, **kwargs):
     """
     node, elem = binsurface(img, nface)
 
@@ -787,7 +787,7 @@ def binsurface(img, nface=3):
     xi, yi, zi = np.unravel_index(id, newdim, order="F")
     node = np.column_stack([xi, yi, zi])
 
-    if nface == 3:
+    if nface == 3 and kwargs.get("clean", True):
         node, elem = meshcheckrepair(node, elem)
 
     return node, elem

iso2mesh/modify.py

@@ -378,6 +378,7 @@ def meshcheckrepair(node, elem, opt=None, **kwargs):
         if status:
             raise RuntimeError(f"jmeshlib command failed: {output}")
         node, elem = readoff(mwpath("post_sclean.off"))
+        elem = np.flipud(elem)
 
     if opt == "meshfix":
         exesuff = fallbackexeext(getexeext(), "meshfix")

iso2mesh/plot.py

@@ -428,7 +428,15 @@ def plotmesh(node, *args, **kwargs):
         if getattr(idx, "size", None) == 0:
             print("Warning: nothing to plot")
             return None
-        handles = plottetra(node, elem[idx, :], opt, *args, **kwargs)
+        if "alpha" not in kwargs and selector:
+            faceselect, elemid = volface(elem[idx, :4])
+            if elem.shape[1] > 4:
+                faceselect = np.hstack(
+                    (faceselect, elem[idx[elemid - 1], 4].reshape(-1, 1))
+                )
+            handles = plotsurf(node, faceselect, opt, *args, **kwargs)
+        else:
+            handles = plottetra(node, elem[idx, :], opt, *args, **kwargs)
 
     if not "hold" in extraarg or not extraarg["hold"] or extraarg["hold"] == "off":
         plt.draw()

iso2mesh/trait.py

@@ -142,10 +142,9 @@ def surfedge(f, junction=None):
 
     openedge = edges[ix[qx], :]
 
-    elemid = None
-    elemid, _ = np.unravel_index(ix[qx], f.shape)
+    elemid, _ = np.unravel_index(ix[qx], f.shape, order="F")
 
-    return openedge, elemid
+    return openedge, elemid + 1
 
 
 # _________________________________________________________________________________________________________
@@ -495,11 +494,11 @@ def layersurf(elem, **kwargs):
     # Process each label
     for i in range(len(labels)):
         if outsideislower == ">=" and labels[i] not in innermost:
-            newface = volface(elem[elem[:, 4] >= labels[i], :4])
+            newface = volface(elem[elem[:, 4] >= labels[i], :4])[0]
         elif outsideislower == "<=" and labels[i] not in innermost:
-            newface = volface(elem[elem[:, 4] <= labels[i], :4])
+            newface = volface(elem[elem[:, 4] <= labels[i], :4])[0]
         else:
-            newface = volface(elem[elem[:, 4] == labels[i], :4])
+            newface = volface(elem[elem[:, 4] == labels[i], :4])[0]
 
         # Add label to faces
         newface = np.hstack((newface, np.full((newface.shape[0], 1), labels[i])))

setup.py

@@ -6,7 +6,7 @@
 setup(
     name="iso2mesh",
     packages=["iso2mesh"],
-    version="0.4.0",
+    version="0.4.1",
     license="GPLv3+",
     description="One-liner 3D Surface and Tetrahedral Mesh Generation Toolbox",
     long_description=readme,

test/run_test.py

@@ -1112,24 +1112,24 @@ def __init__(self, *args, **kwargs):
     def test_binsurface(self):
         no, fc = binsurface(self.im)
         expected_fc = [
-            [10, 4, 1],
-            [7, 10, 1],
-            [12, 6, 5],
-            [11, 5, 4],
-            [2, 3, 9],
-            [1, 2, 8],
-            [11, 12, 5],
-            [10, 11, 4],
-            [2, 9, 8],
-            [1, 8, 7],
-            [8, 9, 12],
-            [7, 8, 11],
-            [6, 3, 2],
-            [5, 2, 1],
-            [8, 12, 11],
-            [7, 11, 10],
-            [5, 6, 2],
             [4, 5, 1],
+            [5, 6, 2],
+            [7, 11, 10],
+            [8, 12, 11],
+            [5, 2, 1],
+            [6, 3, 2],
+            [7, 8, 11],
+            [8, 9, 12],
+            [1, 8, 7],
+            [2, 9, 8],
+            [10, 11, 4],
+            [11, 12, 5],
+            [1, 2, 8],
+            [2, 3, 9],
+            [11, 5, 4],
+            [12, 6, 5],
+            [7, 10, 1],
+            [10, 4, 1],
         ]
         self.assertEqual(fc.tolist(), expected_fc)
 
@@ -1168,7 +1168,7 @@ def test_v2s(self):
     def test_v2m(self):
         no, el, fc = v2m(self.im, 0.5, 0.03, 10)
         self.assertAlmostEqual(sum(elemvolume(no, fc[:, :3])), 5.01, 2)
-        self.assertAlmostEqual(sum(elemvolume(no, el[:, :4])), 0.8786, 4)
+        self.assertAlmostEqual(sum(elemvolume(no, el[:, :4])), 0.8786654361973504, 4)
 
     def test_cgalv2m(self):
         no, el, fc = v2m(
@@ -1184,11 +1184,12 @@ def test_cgalv2m(self):
     @unittest.skipIf(sys.platform.startswith("darwin"), "skip cgals2m")
     def test_cgals2m(self):
         node, elem, face = cgals2m(self.no[:, :3], self.fc[:, :3], 1, 50)
+
         self.assertAlmostEqual(
-            sum(elemvolume(node[:, :3], face[:, :3])) * 0.001, 2.464594888865948, 3
+            sum(elemvolume(node[:, :3], face[:, :3])) * 0.001, 2.456809397314205, 3
         )
         self.assertAlmostEqual(
-            sum(elemvolume(node[:, :3], elem[:, :4])) * 0.001, 3.9983115207040685, 2
+            sum(elemvolume(node[:, :3], elem[:, :4])) * 0.001, 3.983804447661502, 2
         )
 
     def test_v2s_label(self):
@@ -1252,7 +1253,7 @@ def test_surf2volz(self):
             np.arange(19, 42, 1),
             np.arange(19, 42, 1),
         )
-        self.assertAlmostEqual(np.sum(vol.astype(np.float32)) * 0.001, 1.099, 2)
+        self.assertAlmostEqual(np.sum(vol.astype(np.float32)) * 0.001, 1.125, 2)
 
     def test_surf2vol(self):
         vol = surf2vol(
@@ -1262,7 +1263,7 @@ def test_surf2vol(self):
             np.arange(19, 42, 1),
             np.arange(19, 42, 1),
         )
-        self.assertAlmostEqual(np.sum(vol.astype(np.float32)) * 0.001, 1.704, 2)
+        self.assertAlmostEqual(np.sum(vol.astype(np.float32)) * 0.001, 1.72, 2)
 
     def test_surf2vol_fill(self):
         vol = surf2vol(
@@ -1273,11 +1274,11 @@ def test_surf2vol_fill(self):
             np.arange(19, 42, 0.5),
             fill=1,
         )
-        self.assertAlmostEqual(np.sum(vol.astype(np.float32)) * 0.0001, 3.5536, 2)
+        self.assertAlmostEqual(np.sum(vol.astype(np.float32)) * 0.0001, 3.5391, 2)
 
     def test_s2v(self):
         vol = s2v(self.no, self.fc, 100, fill=1)
-        self.assertAlmostEqual(np.sum(vol.astype(np.float32)) * 0.00001, 6.0827, 2)
+        self.assertAlmostEqual(np.sum(vol.astype(np.float32)) * 0.00001, 6.0598097, 2)
 
     def test_mesh2mask(self):
         node = np.array([[0, 0], [2, 0], [2, 1], [0, 1], [1, 0.5]])