[test] add tests for m2v, mesh2vol, mesh2mask, mask has bias

the bias is due to matplotlib rendering under antialiased=False is sensitive to patch orders, reported at
matplotlib/matplotlib#26827

Author: fangq

README.md

@@ -202,21 +202,21 @@ tracked in https://github.com/NeuroJSON/pyiso2mesh/issues/1
 | ✅ `v2s.m` | ✅ tested |                    | ✅ `savedxf.m` | ⭕️ tested |
 | ✅ `s2m.m` | ✅ tested |                    | ✅ `savestl.m` | ⭕️ tested |
 | ✅ `s2v.m` | ✅ tested |                    | ✅ `savebinstl.m` | ⭕️ tested |
-| ✅ `m2v.m` | ⭕️ tested |                    | ✅ `saveinr.m` | ⭕️ tested |
+| ✅ `m2v.m` | ✅️ tested |                    | ✅ `saveinr.m` | ⭕️ tested |
 | ✅ `sms.m` | ✅ tested |                    | ✅ `saveoff.m` | ✅ tested |
 | > Streamlined mesh generation| |            | ⭕️ `savesmf.m` | ⭕️ tested |
 | ✅ `vol2mesh.m` | ✅ tested |               | ✅ `savesurfpoly.m` | ✅ tested |
 | ✅ `vol2surf.m` | ✅ tested |               | ⭕️ `savegts.m` | ⭕️ tested |
 | ✅ `surf2mesh.m` | ✅ tested |              | ⭕️ `readgts.m` | ⭕️ tested |
 | ✅ `surf2vol.m` | ✅ tested |               | ⭕️ `savemsh.m` | ⭕️ tested |
-| ✅ `mesh2vol.m` | ⭕️ tested |               | ⭕️ `savevrml.m` | ⭕️ tested |
+| ✅ `mesh2vol.m` | ✅️ tested |               | ⭕️ `savevrml.m` | ⭕️ tested |
 | > Iso2mesh main function backend| |         | ✅ `readasc.m` | ⭕️ tested |
 | ✅ `binsurface.m` | ✅ tested |             | ⭕️ `readinr.m` | ⭕️ tested |
 | ✅ `cgalv2m.m` | ✅ tested |                | ✅ `readmedit.m` | ⭕️ tested |
 | ✅ `cgals2m.m` | ✅ tested |                | ✅ `readoff.m` | ✅ tested |
 | ✅ `vol2restrictedtri.m` | ✅ tested |      | ⭕️ `readsmf.m` | ⭕️ tested |
 | ✅ `surf2volz.m` | ✅ tested |              | ✅ `readtetgen.m` | ✅ tested |
-| ✅ `mesh2mask.m` | ⭕️ tested |              | ✅ `deletemeshfile.m` | ✅ tested |
+| ✅ `mesh2mask.m` | ✅️ tested |              | ✅ `deletemeshfile.m` | ✅ tested |
 | > Iso2mesh primitive meshing|           |   | ✅ `mcpath.m` | ✅ tested |
 | ✅ `meshabox.m` | ✅ tested |               | ✅ `mwpath.m` | ✅ tested |
 | ✅ `meshasphere.m` | ✅ tested |            | ✅ `savemedit.m` | ✅ tested |

iso2mesh/modify.py

@@ -253,7 +253,7 @@ def qmeshcut(elem, node, value, cutat):
     cutweight = np.abs(cutweight / totalweight[:, np.newaxis])
 
     nodeid = edges[cutedges] - 1
-    nodeid = np.column_stack([nodeid, cutweight[:, 1]])
+    nodeid = np.column_stack([nodeid, cutweight[:, 1]]) + 1
 
     cutpos = (
         node[edges[cutedges, 0] - 1] * cutweight[:, [1]]
@@ -283,7 +283,7 @@ def qmeshcut(elem, node, value, cutat):
         elemid = linecut
         if value.shape[0] == elem.shape[0] and "cutvalue" not in locals():
             cutvalue = value[elemid]
-        return cutpos, cutvalue, facedata, elemid, nodeid
+        return cutpos, cutvalue, facedata, elemid + 1, nodeid
 
     tricut = np.where(etag == 3)[0]
     quadcut = np.where(etag == 4)[0]
@@ -300,7 +300,7 @@ def qmeshcut(elem, node, value, cutat):
 
     facedata = np.vstack([tripatch[:, [0, 1, 2, 2]], quadpatch[:, [0, 1, 3, 2]]])
 
-    return cutpos, cutvalue, facedata, elemid, nodeid
+    return cutpos, cutvalue, facedata, elemid + 1, nodeid
 
 
 def meshcheckrepair(node, elem, opt=None, **kwargs):

iso2mesh/raster.py

@@ -19,7 +19,7 @@
 ##====================================================================================
 
 
-def m2v(*args):
+def m2v(*args, **kwargs):
     """
     Shortcut for mesh2vol, rasterizing a tetrahedral mesh to a volume.
 
@@ -29,10 +29,10 @@ def m2v(*args):
     Returns:
     Volumetric representation of the mesh.
     """
-    return mesh2vol(*args)
+    return mesh2vol(*args, **kwargs)
 
 
-def mesh2vol(node, elem, xi, yi=None, zi=None):
+def mesh2vol(node, elem, xi, yi=None, zi=None, **kwargs):
     """
     mesh2vol(node, elem, xi, yi=None, zi=None)
 
@@ -108,10 +108,10 @@ def mesh2vol(node, elem, xi, yi=None, zi=None):
             continue
 
         if weight is not None:
-            maskz, weightz = mesh2mask(cutpos, facedata, xi, yi)
+            maskz, weightz = mesh2mask(cutpos, facedata, xi, yi, **kwargs)
             weight[:, :, :, i] = weightz
         else:
-            maskz = mesh2mask(cutpos, facedata, xi, yi)[0]
+            maskz = mesh2mask(cutpos, facedata, xi, yi, **kwargs)[0]
 
         idx = ~np.isnan(maskz)
         if nodeval is not None:
@@ -131,7 +131,7 @@ def mesh2vol(node, elem, xi, yi=None, zi=None):
     return mask, weight
 
 
-def mesh2mask(node, face, xi, yi=None, hf=None):
+def mesh2mask(node, face, xi, yi=None, hf=None, **kwargs):
     """
     Fast rasterization of a 2D mesh to an image with triangle index labels.
 
@@ -173,9 +173,7 @@ def mesh2mask(node, face, xi, yi=None, hf=None):
         )
 
     fig = (
-        plt.figure(
-            figsize=(xi.size * 0.01, yi.size * 0.01), dpi=100, layout="compressed"
-        )
+        plt.figure(figsize=(xi.size * 0.01, yi.size * 0.01), dpi=100)
         if hf is None
         else hf
     )
@@ -185,30 +183,28 @@ def mesh2mask(node, face, xi, yi=None, hf=None):
     ax.set_ylim(mn[1], mx[1])
     ax.set_axis_off()
 
-    colors = cm.gray(np.linspace(0, 1, len(face)))
+    colors = cm.jet(np.linspace(0, 1, len(face)))
 
     patches = []
     for i, f in enumerate(face[:, :3]):
-        polygon = Polygon(
-            node[f - 1, :2],
-            closed=True,
-            edgecolor="none",
-            linewidth=0,
-            linestyle="none",
-        )
+        polygon = Polygon(node[f - 1, :2], closed=True, zorder=1)
         patches.append(polygon)
 
     collection = PatchCollection(
-        patches, facecolors=colors, linewidths=0.01, edgecolors="none", edgecolor="face"
+        patches,
+        facecolors=colors,
+        linewidths=0,
+        edgecolor="face",
+        antialiased=(not kwargs.get("edge", True)),
     )
     ax.add_collection(collection)
 
     plt.draw()
     fig.canvas.draw()
     img = np.array(fig.canvas.renderer.buffer_rgba())
-    mask_raw = img[:, :, 0]
-    mask = np.zeros(mask_raw.shape, dtype=np.int32)
-    color_vals = (colors[:, :3] * 255).astype(np.uint8)
+
+    mask = np.zeros(img.shape[:2], dtype=np.int32) * np.nan
+    color_vals = np.floor(colors[:, :3] * 255 + 0.5).astype(np.uint8)
 
     for idx, cval in enumerate(color_vals):
         match = np.all(img[:, :, :3] == cval, axis=-1)

test/run_test.py

@@ -1278,6 +1278,49 @@ 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)
 
+    def test_mesh2mask(self):
+        node = np.array([[0, 0], [2, 0], [2, 1], [0, 1], [1, 0.5]])
+        face = np.array([[1, 2, 5], [2, 3, 5], [3, 4, 5], [4, 1, 5]])
+        xi = np.arange(-0.1, 2.2, 0.1)
+        yi = np.arange(-0.1, 1.2, 0.1)
+        mask, _ = mesh2mask(node, face, xi, yi)
+        counts = np.histogram(mask, bins=np.arange(1, 6) - 0.5)[0]
+
+        self.assertEqual(mask.shape, (xi.size, yi.size))
+        self.assertEqual(np.count_nonzero(np.isnan(mask)), (xi.size + yi.size) * 2 - 4)
+        self.assertEqual(counts.tolist(), [45, 55, 60, 71])
+
+        mask, _ = mesh2mask(node, face, xi, yi, edge=False)
+        self.assertEqual(np.count_nonzero(np.isnan(mask)), 179)
+
+    def test_mesh2vol(self):
+        no1, el1 = meshgrid5(np.arange(1, 4), np.arange(1, 3), np.arange(1, 3))
+        mask, _ = mesh2vol(
+            no1,
+            el1,
+            np.arange(1, 3.1, 0.05),
+            np.arange(1, 2.1, 0.05),
+            np.arange(1, 2, 0.05),
+        )
+        counts = np.histogram(mask, bins=np.arange(0, el1.shape[0] + 2) - 0.5)[0]
+
+        self.assertEqual(np.count_nonzero(np.isnan(mask)), 2657)
+        self.assertTrue(np.all(counts > 100))
+
+    def test_m2v(self):
+        no1, el1 = meshgrid5(np.arange(1, 3), np.arange(1, 3), np.arange(1, 3))
+        mask, _ = m2v(
+            no1,
+            el1,
+            np.arange(1, 3.1, 0.05),
+            np.arange(1, 2.1, 0.05),
+            np.arange(1, 2, 0.05),
+        )
+        counts = np.histogram(mask, bins=np.arange(0, el1.shape[0] + 2) - 0.5)[0]
+
+        self.assertEqual(np.count_nonzero(np.isnan(mask)), 9928)
+        self.assertTrue(np.all(counts > 50))
+
 
 @unittest.skipIf(
     (int(mpl_ver[0]), int(mpl_ver[1])) < (3, 6), "Requires Matplotlib 3.6 or higher"