[doc] update README, accept subplot/parent/hold in plotmesh, 0.3.4

Author: fangq

Makefile

@@ -0,0 +1,18 @@
+
+PY=python3
+
+all: pretty test
+pretty:
+	$(PY) -m black test/*.py iso2mesh/*.py setup.py
+
+test:
+	$(PY) -m unittest test.run_test
+
+report:
+	@echo '====== all imported functions ======'
+	@grep '^\s*[a-z]*,' iso2mesh/__init__.py | sed -e 's/^\s*//g' -e 's/,//g' | sort | uniq -c
+	@echo '====== all tested functions ======'
+	@grep 'def\s*test' test/run_test.py | sed -e 's/\s*def test_//g' -e 's/(self)://g' -e 's/_.*//g' | sort | uniq -c
+
+.DEFAULT_GOAL=all
+.PHONY: all pretty test report

README.md

@@ -1,12 +1,17 @@
-![](https://neurojson.org/wiki/upload/neurojson_banner_long.png)
+![](https://iso2mesh.sourceforge.net/images/iso2mesh_2015_banner.png)
 
 # pyiso2mesh - One-liner 3D Surface and Tetrahedral Mesh Generation Toolbox
 
-* Copyright: (C) Qianqian Fang (2024-2025) <q.fang at neu.edu>, Edward Xu (2024) <xu.ed at northeastern.edu>
-* License: GNU Public License V3 or later
-* Version: 0.3.3
-* URL: [https://pypi.org/project/iso2mesh/](https://pypi.org/project/iso2mesh/)
-* Github: [https://github.com/NeuroJSON/pyiso2mesh](https://github.com/NeuroJSON/pyiso2mesh)
+* **Copyright**: (C) Qianqian Fang (2024-2025) <q.fang at neu.edu>
+* **License**: GNU Public License V3 or later
+* **Version**: 0.3.4
+* **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)
+* **Acknowledgement**: This project is supported by the US National Institute of Health (NIH)
+  grants [U24-NS124027](https://reporter.nih.gov/project-details/10308329) and 
+  [R01-CA204443](https://reporter.nih.gov/project-details/10982160)
+
 
 ![Python Module](https://github.com/NeuroJSON/pyiso2mesh/actions/workflows/build_all.yml/badge.svg)
 
@@ -170,3 +175,140 @@ no, fc, _, _ = v2s(img, 0.5, {'distbound': 0.2})
 ax = plotmesh(no, fc, 'y < 30', alpha=0.5,  edgecolor='none')
 plotmesh(no, fc, 'y > 30', parent = ax)
 ```
+
+## Iso2Mesh function port status
+
+The progress of converting MATLAB-based Iso2Mesh functions to Python is
+tracked in https://github.com/NeuroJSON/pyiso2mesh/issues/1
+
+| Ported | Tested |     | Ported | Tested |
+| ------ | ------ | --- | ------ | ------ |
+| > Streamlined mesh generation - shortcuts| |  | > File I/O | |
+| ✅ `v2m.m` | ✅ tested |                    | ✅ `saveasc.m` | ⭕️ tested |
+| ✅ `v2s.m` | ✅ tested |                    | ✅ `savedxf.m` | ⭕️ tested |
+| ✅ `s2m.m` | ✅ tested |                    | ✅ `savestl.m` | ⭕️ tested |
+| ✅ `s2v.m` | ⭕️ tested |                    | ✅ `savebinstl.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 |
+| > 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 |
+| > Iso2mesh primitive meshing functions| |   | ✅ `mcpath.m` | ✅ tested |
+| ✅ `meshabox.m` | ✅ tested |               | ✅ `mwpath.m` | ✅ tested |
+| ✅ `meshasphere.m` | ✅ tested |            | ✅ `savemedit.m` | ✅ tested |
+| ✅ `meshanellip.m` | ✅ tested |            | ⭕️ `savejson.m` | ⭕️ tested |
+| ✅ `meshunitsphere.m` | ✅ tested |         | ⭕️ `loadjson.m` | ⭕️ tested |
+| ✅ `meshacylinder.m` | ✅ tested |          | ⭕️ `saveubjson.m` | ⭕️ tested |
+| ✅ `meshgrid5.m` | ✅ tested |              | ⭕️ `loadubjson.m` | ⭕️ tested |
+| ✅ `meshgrid6.m` | ✅ tested |              | ⭕️ `loadmsgpack.m` | ⭕️ tested |
+| ✅ `latticegrid.m` | ✅ tested |            | ⭕️ `savemsgpack.m` | ⭕️ tested |
+| ✅ `extrudecurve.m` | ⭕️ tested |           | ⭕️ `savebj.m` | ⭕️ tested |
+| ✅ `meshcylinders.m` | ✅ tested |          | ⭕️ `loadbj.m` | ⭕️ tested |
+| > Mesh decomposition and query| |           | ⭕️ `savemphtxt.m` | ⭕️ tested |
+| ✅ `finddisconnsurf.m` | ✅ tested |        | ⭕️ `savetetgenele.m` | ⭕️ tested |
+| ✅ `surfedge.m` | ✅ tested |               | ⭕️ `savetetgennode.m` | ⭕️ tested |
+| ✅ `volface.m` | ✅ tested |                | ⭕️ `saveabaqus.m` | ⭕️ tested |
+| ✅ `extractloops.m` | ✅ tested |           | ⭕️ `savenirfast.m` | ⭕️ tested |
+| ✅ `meshconn.m` | ✅ tested |               | ⭕️ `readnirfast.m` | ⭕️ tested |
+| ✅ `meshcentroid.m` | ✅ tested |           | ⭕️ `readnifti.m` | ⭕️ tested |
+| ✅ `nodevolume.m` | ✅ tested |             | ⭕️ `readmptiff.m` | ⭕️ tested |
+| ✅ `elemvolume.m` | ✅ tested |             | ⭕️ `loadjsnirf.m` | ⭕️ tested |
+| ✅ `neighborelem.m` | ✅ tested |           | ⭕️ `savejsnirf.m` | ⭕️ tested |
+| ✅ `layersurf.m` | ⭕️ tested |              | ⭕️ `loadsnirf.m` | ⭕️ tested |
+| ✅ `faceneighbors.m` | ✅ tested |          | ⭕️ `savesnirf.m` | ⭕️ tested |
+| ✅ `edgeneighbors.m` | ✅ tested |          | ⭕️ `readobjmesh.m` | ⭕️ tested |
+| ✅ `maxsurf.m` | ⭕️ tested |                | ⭕️ `loadjmesh.m` | ⭕️ tested |
+| ✅ `flatsegment.m` | ⭕️ tested |            | ⭕️ `readobjmesh.m` | ⭕️ tested |
+| ✅ `orderloopedge.m` | ⭕️ tested |          | > Volumetric image pre-processing| |
+| ✅ `mesheuler.m` | ✅ tested |              | ⭕️ `bwislands.m` | ⭕️ tested |
+| ✅ `bbxflatsegment.m` | ⭕️ tested |         | ⭕️ `fillholes3d.m` | ⭕️ tested |
+| ✅ `surfplane.m` | ⭕️ tested |              | ⭕️ `deislands2d.m` | ⭕️ tested |
+| ✅ `surfinterior.m` | ⭕️ tested |           | ⭕️ `deislands3d.m` | ⭕️ tested |
+| ✅ `surfpart.m` | ⭕️ tested |               | ⭕️ `ndgaussian.m` | ⭕️ tested |
+| ✅ `surfseeds.m` | ⭕️ tested |              | ⭕️ `ndimfilter.m` | ⭕️ tested |
+| ✅ `meshquality.m` | ✅ tested |            | ⭕️ `imedge3d.m` | ⭕️ tested |
+| ✅ `meshedge.m` | ✅ tested |               | ⭕️ `internalpoint.m` | ⭕️ tested |
+| ✅ `meshface.m` | ✅ tested |               | ⭕️ `smoothbinvol.m` | ⭕️ tested |
+| ✅ `surfacenorm.m` | ✅ tested |            | ⭕️ `thickenbinvol.m` | ⭕️ tested |
+| ✅ `nodesurfnorm.m` | ✅ tested |           | ⭕️ `thinbinvol.m` | ⭕️ tested |
+| ✅ `uniqedges.m` | ✅ tested |              | ⭕️ `maskdist.m` | ⭕️ tested |
+| ✅ `uniqfaces.m` | ✅ tested |              | > Mesh plotting| |
+| ✅ `advancefront.m` | ⭕️ tested |           | ✅ `plotmesh.m` | ✅ tested |
+| ✅ `innersurf.m` | ⭕️ tested |              | ✅ `plotsurf.m` | ✅ tested |
+| ✅ `outersurf.m` | ⭕️ tested |              | ✅ `plottetra.m` | ✅ tested |
+| ✅ `surfvolume.m` | ✅ tested |             | ✅ `plotedges.m` | ✅ tested |
+| ✅ `insurface.m` | ✅ tested |              | ✅ `qmeshcut.m` | ✅ tested |
+| > Mesh processing and reparing| |           | > Miscellaneous functions| |
+| ✅ `meshcheckrepair.m` | ✅ tested |        | ⭕️ `surfdiffuse.m` | ⭕️ tested |
+| ✅ `meshreorient.m` | ✅ tested |           | ⭕️ `volmap2mesh.m` | ⭕️ tested |
+| ✅ `removedupelem.m` | ✅ tested |          | ⭕️ `isoctavemesh.m` | ⭕️ tested |
+| ✅ `removedupnodes.m` | ✅ tested |         | ⭕️ `getvarfrom.m` | ⭕️ tested |
+| ✅ `removeisolatednode.m` | ✅ tested |     | ✅ `raytrace.m` | ✅ tested |
+| ✅ `removeisolatedsurf.m` | ⭕️ tested |     | ⭕️ `linextriangle.m` | ⭕️ tested |
+| ✅ `surfaceclean.m` | ⭕️ tested |           | ⭕️ `getplanefrom3pt.m` | ⭕️ tested |
+| ✅ `getintersecttri.m` | ⭕️ tested |        | ✅ `getexeext.m` | ✅ tested |
+| ✅ `delendelem.m` | ⭕️ tested |             | ✅ `fallbackexeext.m` | ✅ tested |
+| ✅ `surfreorient.m` | ✅ tested |           | ⭕️ `iso2meshver.m` | ⭕️ tested |
+| > Mesh registration - Metch Toolbox| |      | ⭕️ `raysurf.m` | ⭕️ tested |
+| ✅ `proj2mesh.m` | ⭕️ tested |              | ⭕️ `getoptkey.m` | ⭕️ tested |
+| ✅ `dist2surf.m` | ⭕️ tested |              | ✅ `rotatevec3d.m` | ⭕️ tested |
+| ✅ `regpt2surf.m` | ⭕️ tested |             | ✅ `rotmat2vec.m` | ⭕️ tested |
+| ✅ `affinemap.m` | ⭕️ tested |              | ✅ `varargin2struct.m` | ⭕️ tested |
+| > Polyline handling| |                      | ✅ `jsonopt.m` | ⭕️ tested |
+| ⭕️ `slicesurf.m` | ⭕️ tested |              | ⭕️ `mergestruct.m` | ⭕️ tested |
+| ⭕️ `slicesurf3.m` | ⭕️ tested |             | ⭕️ `orthdisk.m` | ⭕️ tested |
+| ⭕️ `polylinelen.m` | ⭕️ tested |            | ⭕️ `nestbracket2dim.m` | ⭕️ tested |
+| ⭕️ `polylinesimplify.m` | ⭕️ tested |       | ⭕️ `memmapstream.m` | ⭕️ tested |
+| ⭕️ `polylineinterp.m` | ⭕️ tested |         | ⭕️ `aos2soa.m` | ⭕️ tested |
+| ⭕️ `closestnode.m` | ⭕️ tested |            | ⭕️ `soa2aos.m` | ⭕️ tested |
+| > Mesh resampling and optimization| |
+| ✅ `meshresample.m` | ✅ tested |
+| ✅ `remeshsurf.m` | ⭕️ tested |
+| ✅ `smoothsurf.m` | ✅ tested |
+| ✅ `sortmesh.m` | ⭕️ tested |
+| ✅ `mergemesh.m` | ⭕️ tested |
+| ✅ `meshrefine.m` | ⭕️ tested |
+| ✅ `mergesurf.m` | ⭕️ tested |
+| ✅ `surfboolean.m` | ✅ tested |
+| ✅ `fillsurf.m` | ⭕️ tested |
+| ✅ `highordertet.m` | ⭕️ tested |
+| ✅ `elemfacecenter.m` | ⭕️ tested |
+| ✅ `barydualmesh.m` | ⭕️ tested |
+| ✅ `meshinterp.m` | ⭕️ tested |
+| ✅ `meshremap.m` | ⭕️ tested |
+| ✅ `extrudesurf.m` | ⭕️ tested |
+
+
+
+## Acknowledgement
+
+The `pyiso2mesh` module was converted from the MATLAB/Octave version of
+of Iso2Mesh (https://github.com/fangq/iso2mesh) written by the same author.
+
+We utilized large-language-model (LLM) and AI chatbot in the initial MATLAB-to-Python
+conversion, with specific instructions to **faithfully replicate** the algorithms
+developed in the MATLAB code to avoid introducing external copyrighted materials
+into this toolbox. With the assumption that the AI chatbot functions primarily
+as an automated syntax translator without altering the originality of the code,
+the upstream author of the original MATLAB-based Iso2Mesh retains the full copyright
+of this derived Python module.
+
+The initial translation was further manually tested, adjusted and restructured
+to produce matching outputs as the original MATLAB toolbox. Dr. Edward Xu
+<xu.ed at northeastern.edu> had contributed to the initial conversion and
+testing of a subset of functions in the `geometry.py` and `trait.py` units.
+
+The development of this software is supported by the US National Institute of Health (NIH)
+under research awards [U24-NS124027](https://reporter.nih.gov/project-details/10308329) and 
+[R01-CA204443](https://reporter.nih.gov/project-details/10982160), with the author,
+Dr. Qianqian Fang, serving as the principle investigator on both awards.

iso2mesh/__init__.py

@@ -137,7 +137,7 @@
     barycentricgrid,
 )
 
-__version__ = "0.3.3"
+__version__ = "0.3.4"
 __all__ = [
     "advancefront",
     "barycentricgrid",

iso2mesh/plot.py

@@ -38,7 +38,7 @@ def plotsurf(node, face, *args, **kwargs):
 
     sc = np.random.rand(10, 3)
 
-    ax = _createaxis(*args, **kwargs)
+    ax = plt.gca()
 
     h = {"fig": [], "ax": [], "obj": []}
     h["fig"].append(plt.gcf())
@@ -180,7 +180,7 @@ def plottetra(node, elem, *args, **kwargs):
 
     np.random.seed(randseed)
 
-    ax = _createaxis(*args, **kwargs)
+    ax = plt.gca()
 
     h = {"fig": [], "ax": [], "obj": []}
     h["fig"].append(plt.gcf())
@@ -264,7 +264,7 @@ def plotedges(node, edges, *args, **kwargs):
     edlen = edges.shape[0]
     rng_state = np.random.get_state()
 
-    ax = _createaxis(*args, **kwargs)
+    ax = plt.gca()
 
     hh = {"fig": [], "ax": [], "obj": []}
     hh["fig"].append(plt.gcf())
@@ -307,14 +307,6 @@ def plotedges(node, edges, *args, **kwargs):
         edges = edges.astype(int) - 1  # 1-based to 0-based
 
         if node.shape[1] >= 3:
-            ax = plt.gca()
-
-            if ax.name != "3d":
-                plt.figure()  # Create a new figure
-                ax = plt.gcf().add_subplot(
-                    projection="3d"
-                )  # Add 3D axes to the current figure
-
             segments = [[node[start], node[end]] for start, end in edges]
             h = Line3DCollection(segments, **kwargs)
             ax.add_collection3d(h)
@@ -404,9 +396,6 @@ def plotmesh(node, *args, **kwargs):
         (h,) = ax.plot(x[idx], y[idx], z[idx], *opt, **kwargs)
         handles["obj"].append(h)
         _autoscale_3d(ax, node)
-        if not "hold" in extraarg or not extraarg["hold"] or extraarg["hold"] == "off":
-            plt.show(block=False)
-        return handles
 
     # Plot surface mesh
     if face is not None:
@@ -440,6 +429,7 @@ def plotmesh(node, *args, **kwargs):
         handles = plottetra(node, elem[idx, :], opt, *args, **kwargs)
 
     if not "hold" in extraarg or not extraarg["hold"] or extraarg["hold"] == "off":
+        plt.draw()
         plt.show(block=False)
 
     return handles
@@ -456,23 +446,29 @@ def _autoscale_3d(ax, points):
 
 def _createaxis(*args, **kwargs):
     subplotid = kwargs.get("subplot", 111)
-    docreate = False if len(args) == 0 else args[0]
+    docreate = True
+    fig = None
 
     if "parent" in kwargs:
-        ax = kwargs["parent"]
-        if isinstance(ax, dict):
-            ax = ax["ax"][-1]
-        elif isinstance(ax, list):
-            ax = ax[-1]
+        hh = kwargs["parent"]
+        if isinstance(hh, dict):
+            fig = hh["fig"][0]
+            ax = hh["ax"][-1]
+        elif isinstance(hh, list):
+            ax = hh[-1]
+        if "subplot" in kwargs and fig:
+            ax = fig.add_subplot(subplotid, projection="3d")
     elif not docreate and len(plt.get_fignums()) > 0 and len(plt.gcf().axes) > 0:
-        ax = plt.gcf().axes[-1]
+        if not fig:
+            fig = plt.gcf()
+        ax = fig.axes[-1]
     else:
-        if docreate:
-            plt.figure()
-        ax = plt.gcf().add_subplot(subplotid, projection="3d")
+        if docreate and not fig:
+            fig = plt.figure()
+        ax = fig.add_subplot(subplotid, projection="3d")
 
     if ax.name != "3d":
-        plt.figure()  # Create a new figure
-        ax = plt.gcf().add_subplot(subplotid, projection="3d")
+        fig = plt.figure()  # Create a new figure
+        ax = fig.add_subplot(subplotid, projection="3d")
 
     return ax

iso2mesh/trait.py

@@ -347,7 +347,7 @@ def nodevolume(node, elem, evol=None):
 
     # Loop through each element and accumulate the volumes
     for i in range(elemnum):
-        nodevol[elem[i, :dim]] += evol[i]
+        nodevol[elem[i, :dim] - 1] += evol[i]
 
     # Divide by the dimensionality to get the final node volumes
     nodevol /= dim

setup.py

@@ -6,9 +6,9 @@
 setup(
     name="iso2mesh",
     packages=["iso2mesh"],
-    version="0.3.3",
-    license='GPLv3+',
-    description="Image-based 3D Surface and Volumetric Mesh Generator",
+    version="0.3.4",
+    license="GPLv3+",
+    description="One-liner 3D Surface and Tetrahedral Mesh Generation Toolbox",
     long_description=readme,
     long_description_content_type="text/markdown",
     author="Qianqian Fang",
@@ -18,12 +18,20 @@
     keywords=[
         "Iso2Mesh",
         "Mesh generation",
+        "Finite element method",
         "FEM",
+        "FEA",
         "MATLAB",
         "Mesh-based Monte Carlo",
         "CGAL",
         "Tetgen",
         "CSG",
+        "Boolean",
+        "Surface mesh",
+        "Tetrahedral mesh",
+        "Modeling",
+        "Multiphysics",
+        "Image processing",
     ],
     platforms="any",
     install_requires=["numpy>=1.8.0", "matplotlib"],

test/run_test.py

@@ -577,9 +577,58 @@ def test_elemvolume(self):
         result = np.unique(elemvolume(self.no, self.el).round(decimals=6)).tolist()
         self.assertEqual(result, expected)
 
+    def test_nodevolume(self):
+        expected = np.sum(elemvolume(self.no, self.el))
+        result = np.sum(nodevolume(self.no, self.el))
+        self.assertAlmostEqual(result, expected, 7)
+
     def test_surfvolume(self):
         self.assertEqual(surfvolume(self.no, self.fc), 1)
 
+    def test_surfacenorm(self):
+        snorm = surfacenorm(self.no, self.fc)
+        expected = [
+            [0.0, 0.0, -1.0],
+            [0.0, -1.0, 0.0],
+            [0.0, 0.0, -1.0],
+            [-1.0, 0.0, 0.0],
+            [0.0, -1.0, 0.0],
+            [-1.0, 0.0, 0.0],
+            [1.0, 0.0, 0.0],
+            [1.0, 0.0, 0.0],
+            [0.0, 1.0, 0.0],
+            [0.0, 1.0, 0.0],
+            [0.0, -1.0, 0.0],
+            [-1.0, 0.0, 0.0],
+            [0.0, -1.0, 0.0],
+            [-1.0, 0.0, 0.0],
+            [1.0, 0.0, 0.0],
+            [1.0, 0.0, 0.0],
+            [0.0, 1.0, 0.0],
+            [0.0, 1.0, 0.0],
+            [0.0, 0.0, 1.0],
+            [0.0, 0.0, 1.0],
+        ]
+        self.assertEqual(snorm.tolist(), expected)
+
+    def test_nodesurfnorm(self):
+        nnorm = nodesurfnorm(self.no, self.fc)
+        expected = [
+            [-0.57735, -0.57735, -0.57735],
+            [0.816497, -0.408248, -0.408248],
+            [-0.408248, 0.816497, -0.408248],
+            [0.408248, 0.408248, -0.816497],
+            [-0.707107, -0.707107, 0.0],
+            [0.707107, -0.707107, 0.0],
+            [-0.707107, 0.707107, 0.0],
+            [0.707107, 0.707107, 0.0],
+            [-0.408248, -0.408248, 0.816497],
+            [0.408248, -0.816497, 0.408248],
+            [-0.816497, 0.408248, 0.408248],
+            [0.57735, 0.57735, 0.57735],
+        ]
+        self.assertEqual(nnorm.round(6).tolist(), expected)
+
     def test_insurface(self):
         pts = np.array([[1.5, -0.9, 2.1], [1, 0, 2], [-1, 0, 2], [1.2, 0, 2.5]])
         expected = [1, 1, 0, 1]