[release] improve brain1020 to be more robust,

porting fangq/brain2mesh@2e2760e to handle neck/shoulder and multi-loop cuts

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.5.1
+* **Version**: 0.5.2
 * **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

@@ -205,6 +205,7 @@
     closestnode,
     polylineinterp,
     polylinesimplify,
+    maxloop,
 )
 
 from .raster import (
@@ -241,7 +242,7 @@
     thinbinvol,
 )
 
-__version__ = "0.5.1"
+__version__ = "0.5.2"
 __all__ = [
     "advancefront",
     "barycentricgrid",
@@ -362,6 +363,7 @@
     "closestnode",
     "polylineinterp",
     "polylinesimplify",
+    "maxloop",
     "affinemap",
     "meshinterp",
     "meshremap",

iso2mesh/brain2mesh.py

@@ -37,7 +37,13 @@
 )
 from iso2mesh.geometry import meshabox
 from iso2mesh.volume import volgrow, fillholes3d
-from iso2mesh.line import polylineinterp, polylinelen, polylinesimplify, closestnode
+from iso2mesh.line import (
+    polylineinterp,
+    polylinelen,
+    polylinesimplify,
+    closestnode,
+    maxloop,
+)
 from iso2mesh.plot import plotmesh
 
 from typing import Tuple, Dict, Union, Optional, List
@@ -1367,11 +1373,54 @@ def brain1020(
                 landmarks["cz"],
             ]
         )
-
+    else:
+        landmarks = {
+            "nz": initpoints[0, :],
+            "iz": initpoints[1, :],
+            "lpa": initpoints[2, :],
+            "rpa": initpoints[3, :],
+            "cz": initpoints[4, :],
+        }
     # Convert tetrahedral mesh into a surface mesh
     if face.shape[1] >= 4:
         face = volface(face[:, :4])[0]  # Use first 4 columns for tetrahedral faces
 
+    if kwargs.get("clean", 1):
+        # Find the bounding box of the top part of the head, and remove all other triangles
+        p0 = landmarks.copy()
+
+        v_ni = p0["nz"] - p0["iz"]
+        v_lr = p0["lpa"] - p0["rpa"]
+        v_cz0 = np.cross(v_ni, v_lr)
+        v_cz0 = v_cz0 / np.linalg.norm(v_cz0)
+        v_cn = p0["cz"] - p0["nz"]
+        d_czlpa = np.dot(p0["cz"] - p0["lpa"], v_cz0)
+        d_cznz = np.dot(v_cn, v_cz0)
+
+        if abs(d_czlpa) > abs(d_cznz):  # if lpa is further away from cz than nz
+            # Move nz to the same level as lpa, can also add rpa
+            p0["nz"] = p0["nz"] - v_cz0 * (abs(d_czlpa) - abs(d_cznz))
+            # Move iz to the same level as lpa, can also add rpa
+            p0["iz"] = p0["iz"] - v_cz0 * (abs(d_czlpa) - abs(d_cznz))
+
+        v_cz = d_cznz * v_cz0
+        bbx0 = p0["nz"] - 0.6 * v_lr - 0.1 * v_cz + 0.1 * v_ni
+
+        # Calculate mesh centroids
+        c0 = meshcentroid(node, face)
+
+        # Calculate distance from bounding box plane
+        v_cz0_rep = np.tile(
+            v_cz0, (face.shape[0], 1)
+        )  # repmat(v_cz0, size(face, 1), 1)
+        bbx0_rep = np.tile(bbx0, (face.shape[0], 1))  # repmat(bbx0, size(face, 1), 1)
+        dz = np.sum(v_cz0_rep * (c0 - bbx0_rep), axis=1)
+
+        # Filter faces - keep only those with dz > 0
+        face = face[dz > 0, :]
+
+        del p0, bbx0, c0, dz
+
     # Remove nodes not located in the surface
     node, face, _ = removeisolatednode(node, face)
 
@@ -1387,16 +1436,6 @@ def brain1020(
         print("Initial points:")
         print(initpoints)
 
-    # Save input initpoints to landmarks output, cz is not finalized
-    if initpoints.shape[0] >= 5:
-        landmarks = {
-            "nz": initpoints[0, :],
-            "iz": initpoints[1, :],
-            "lpa": initpoints[2, :],
-            "rpa": initpoints[3, :],
-            "cz": initpoints[4, :],
-        }
-
     # At this point, initpoints contains {nz, iz, lpa, rpa, cz0}
     # Plot the head mesh
     if showplot:
@@ -1417,7 +1456,8 @@ def brain1020(
     while np.linalg.norm(initpoints[4, :] - lastcz) > tol and cziter < maxcziter:
 
         # Step 1: nz, iz and cz0 to determine saggital reference curve
-        nsagg = slicesurf(node, face, initpoints[[0, 1, 4], :])
+        nsagg, curveloop, _ = slicesurf(node, face, initpoints[[0, 1, 4], :], full=True)
+        nsagg = nsagg[maxloop(curveloop) - 1, :]
 
         # Step 1.1: get cz1 as the mid-point between iz and nz
         slen, nsagg, _ = polylinelen(
@@ -1429,7 +1469,11 @@ def brain1020(
         initpoints[4, :] = cz[0, :]
 
         # Step 1.2: lpa, rpa and cz1 to determine coronal reference curve, update cz1
-        curves["cm"] = slicesurf(node, face, initpoints[[2, 3, 4], :])
+        curves["cm"], curveloop, _ = slicesurf(
+            node, face, initpoints[[2, 3, 4], :], full=True
+        )
+        curves["cm"] = curves["cm"][maxloop(curveloop) - 1, :]
+
         len_cm, curves["cm"], _ = polylinelen(
             curves["cm"], initpoints[2, :], initpoints[3, :], initpoints[4, :]
         )
@@ -1458,7 +1502,11 @@ def brain1020(
     )
     landmarks["cm"] = coro  # t7, c3, cz, c4, t8
 
-    curves["sm"] = slicesurf(node, face, initpoints[[0, 1, 4], :])
+    curves["sm"], curveloop, _ = slicesurf(
+        node, face, initpoints[[0, 1, 4], :], full=True
+    )
+    curves["sm"] = curves["sm"][maxloop(curveloop) - 1, :]
+
     slen, curves["sm"], _ = polylinelen(
         curves["sm"], initpoints[0, :], initpoints[1, :], initpoints[4, :]
     )
@@ -1475,6 +1523,8 @@ def brain1020(
         landmarks["cm"][0, :],
         landmarks["sm"][-1, :],
         perc2 * 2,
+        0,
+        maxloop=1,
     )
 
     # Step 4: fpz, t8 and oz to determine right 10% axial reference curve
@@ -1485,6 +1535,8 @@ def brain1020(
         landmarks["cm"][-1, :],
         landmarks["sm"][-1, :],
         perc2 * 2,
+        0,
+        maxloop=1,
     )
 
     # Show plots of the landmarks
@@ -1611,6 +1663,8 @@ def brain1020(
             landmarks["sm"][idxcz - 1 - i, :],
             landmarks["aar"][i - 1, :],
             step,
+            0,
+            maxloop=1,
         )
 
         landmarks[f"cal_{i}"] = cal_landmarks
@@ -1661,6 +1715,8 @@ def brain1020(
             landmarks["sm"][idxcz - 1 + i, :],
             landmarks["apr"][i - 1, :],
             step,
+            0,
+            maxloop=1,
         )
 
         landmarks[f"cpl_{i}"] = cpl_landmarks
@@ -1696,15 +1752,29 @@ def brain1020(
             landmarks["papl"],
             curves["papl"],
         ) = slicesurf3(
-            node, face, landmarks["nz"], landmarks["lpa"], landmarks["iz"], perc2 * 2
+            node,
+            face,
+            landmarks["nz"],
+            landmarks["lpa"],
+            landmarks["iz"],
+            perc2 * 2,
+            0,
+            maxloop=1,
         )
         (
             landmarks["paar"],
             curves["paar"],
             landmarks["papr"],
             curves["papr"],
         ) = slicesurf3(
-            node, face, landmarks["nz"], landmarks["rpa"], landmarks["iz"], perc2 * 2
+            node,
+            face,
+            landmarks["nz"],
+            landmarks["rpa"],
+            landmarks["iz"],
+            perc2 * 2,
+            0,
+            maxloop=1,
         )
 
         if showplot:

iso2mesh/line.py

@@ -11,6 +11,7 @@
     "closestnode",
     "polylineinterp",
     "polylinesimplify",
+    "maxloop",
 ]
 ##====================================================================================
 ## dependent libraries
@@ -349,3 +350,69 @@ def segvec(n1, n2):
         length = np.array([])
 
     return newnodes, length
+
+
+def maxloop(curveloop: np.ndarray) -> np.ndarray:
+    """
+    Return the curve segment that has the largest number of nodes
+
+    Author: Qianqian Fang, <q.fang at neu.edu>
+    Python conversion: Preserves exact algorithm from MATLAB version
+
+    Parameters:
+    -----------
+    curveloop : ndarray (1D)
+        Curves defined by 1-based node indices, separated by nan.
+        The values in this array are node indices (1-based from MATLAB),
+        not Python array positions. NaN values serve as segment separators.
+
+    Returns:
+    --------
+    newloop : ndarray (1D)
+        The 1-based node indices defining the longest segment.
+        These are the actual node index values, maintaining 1-based indexing.
+
+    Notes:
+    ------
+    This function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
+
+    Important: curveloop and newloop contain 1-based node indices as data values.
+    The indexing for array access (loopend positions) uses Python's 0-based indexing,
+    but the actual node index values remain 1-based throughout.
+    """
+
+    newloop = curveloop.copy()
+
+    # Find array positions where nan occurs - equivalent to find(isnan(curveloop))
+    # loopend contains Python 0-based positions in the curveloop array where NaN appears
+    loopend = np.where(np.isnan(curveloop))[0]
+
+    if len(loopend) > 1:  # length(loopend) > 1
+        # Calculate segment lengths - equivalent to [loopend(1), diff(loopend)]
+        # MATLAB: seglen = [loopend(1), diff(loopend)];
+        # seglen[0] = number of elements before first NaN
+        # seglen[i] = number of elements between consecutive NaNs
+        seglen = np.concatenate(([loopend[0]], np.diff(loopend)))
+
+        # Find maximum segment length and its location
+        # MATLAB: [maxlen, maxloc] = max(seglen);
+        maxloc = np.argmax(seglen)
+        maxlen = seglen[maxloc]
+
+        # Prepend 0 to loopend - equivalent to loopend = [0 loopend];
+        # This represents a virtual NaN at position -1 for easier indexing
+        loopend = np.concatenate(([-1], loopend))
+
+        # Extract the longest segment
+        # MATLAB: newloop = curveloop((loopend(maxloc)+1):(loopend(maxloc+1)-maxloc));
+        # We're extracting array elements, so we use Python 0-based array indexing
+        # But the VALUES we extract are 1-based node indices that we keep as-is
+        start_idx = loopend[maxloc] + 1  # Position after the NaN (or start)
+        end_idx = loopend[maxloc + 1]  # Position of the next NaN
+        newloop = curveloop[start_idx:end_idx]
+
+    # Remove any remaining nan values - equivalent to newloop(isnan(newloop)) = [];
+    # The node index values in newloop remain 1-based
+    newloop = newloop[~np.isnan(newloop)]
+
+    return newloop.astype(int)

iso2mesh/modify.py

@@ -44,7 +44,13 @@
     surfedge,
     extractloops,
 )
-from iso2mesh.line import getplanefrom3pt, polylinesimplify, polylinelen, polylineinterp
+from iso2mesh.line import (
+    getplanefrom3pt,
+    polylinesimplify,
+    polylinelen,
+    polylineinterp,
+    maxloop,
+)
 
 ##====================================================================================
 ## implementations
@@ -1156,7 +1162,7 @@ def slicesurf(node, face, *args, **kwargs):
     return bcutpos, bcutloop, bcutvalue
 
 
-def slicesurf3(node, elem, p1, p2, p3, step=None, minangle=None):
+def slicesurf3(node, elem, p1, p2, p3, step=None, minangle=None, **kwargs):
     """
     slicesurf3(node, elem, p1, p2, p3, step=None, minangle=None)
 
@@ -1187,7 +1193,9 @@ def slicesurf3(node, elem, p1, p2, p3, step=None, minangle=None):
     """
 
     # Slice full curve through p1-p2-p3
-    fullcurve = slicesurf(node, elem, np.vstack((p1, p2, p3)))
+    fullcurve, curveloop, _ = slicesurf(node, elem, np.vstack((p1, p2, p3)), full=True)
+    if kwargs.get("maxloop", 0):
+        fullcurve = fullcurve[maxloop(curveloop) - 1, :]
 
     # Optional simplification
     if minangle is not None and minangle > 0:

setup.py

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