[brain2mesh] enable decoupling in pass 2 if pass 1 fails

Author: fangq

iso2mesh/brain2mesh.py

@@ -354,6 +354,7 @@ def normalizer(x):
             continue
         if (loop == 1) and ("label_elem" not in locals()):
             if "bone_n" in locals() and "skin_n" in locals():
+                print("decoupling bone and skin surfaces")
                 bone_n, bone_f = surfboolean(
                     bone_n[:, :3],
                     bone_f[:, :3],
@@ -362,14 +363,17 @@ def normalizer(x):
                     skin_f[:, :3],
                 )
             if "bone_n" in locals() and "csf_n" in locals():
+                print("decoupling csf and bone surfaces")
                 csf_n, csf_f = surfboolean(
                     csf_n[:, :3], csf_f[:, :3], "decouple", bone_n[:, :3], bone_f[:, :3]
                 )
             if "pial_n" in locals() and "csf_n" in locals():
+                print("decoupling pial and csf surfaces")
                 pial_n, pial_f = surfboolean(
                     pial_n[:, :3], pial_f[:, :3], "decouple", csf_n[:, :3], csf_f[:, :3]
                 )
             if "pial_n" in locals() and "wm_n" in locals():
+                print("decoupling wm and pial surfaces")
                 wm_n, wm_f = surfboolean(
                     wm_n[:, :3], wm_f[:, :3], "decouple", pial_n[:, :3], pial_f[:, :3]
                 )
@@ -452,18 +456,27 @@ def normalizer(x):
             return brain_n, brain_el, brain_f
 
         # Generates a coarse tetrahedral mesh of the combined tissues
+        final_surf_n, final_surf_f = meshcheckrepair(final_surf_n, final_surf_f, "dup")
         try:
             final_n, final_e, _ = s2m(
-                final_surf_n, final_surf_f, 1.0, maxvol, "tetgen1.5", None, None, "-A"
+                final_surf_n,
+                final_surf_f,
+                1.0,
+                maxvol,
+                "tetgen1.5",
+                None,
+                None,
+                "-YY -A",
             )
-        except:
+        except RuntimeError as e:
             print(
-                "volumetric mesh generation failed, returning the intermediate surface model only"
+                f"volumetric mesh generation failed with error: {e}, returning the intermediate surface model only"
             )
-            brain_n = final_surf_n
-            brain_f = final_surf_f
-            brain_el = np.array([])
-            return brain_n, brain_el, brain_f
+            continue
+            # brain_n = final_surf_n
+            # brain_f = final_surf_f
+            # brain_el = np.array([])
+            # return brain_n, brain_el, brain_f
 
         # Removes the elements that are part of the box, but not the brain/head
         if dotruncate == 1 or isinstance(dotruncate, str):
@@ -1761,7 +1774,7 @@ def brain1020(
 
 
 def label2tpm(
-    vol: np.ndarray, names: Optional[Union[List[str], Dict[int, str]]] = None
+    vol: np.ndarray, names: Optional[Union[List[str], Dict[int, str]]] = None, **kwargs
 ) -> Dict[str, np.ndarray]:
     """
     Converting a multi-label volume to binary tissue probabilistic maps (TPMs)
@@ -1824,6 +1837,17 @@ def label2tpm(
         # MATLAB: seg.(nm) = uint8(vol == val(i))
         seg[nm] = (vol == val[i]).astype(np.uint8)
 
+    if kwargs.get("sigma", 0) > 0:
+        from scipy.ndimage import gaussian_filter
+
+        segsum = np.zeros_like(seg[next(iter(seg))], dtype=np.float32)
+        for key in seg:
+            seg[key] = gaussian_filter(seg[key].astype(np.float32), sigma=1)
+            segsum += seg[key]
+
+        for key in seg:
+            np.divide(seg[key], segsum, out=seg[key], where=segsum != 0)
+
     return seg
 
 

iso2mesh/core.py

@@ -983,7 +983,7 @@ def getintersecttri(tmppath):
     eid = []
     if status == 0:
         ids = re.findall(r" #([0-9]+) ", str_output)
-        eid = [int(id[0]) for id in ids]
+        eid = [int(id) for id in ids]
 
     eid = np.unique(eid)
     return eid

iso2mesh/modify.py

@@ -491,8 +491,14 @@ def removeisolatednode(node, elem, face=None):
         Face list of the resulting mesh.
     """
 
+    zeroindex = None
+
     oid = np.arange(node.shape[0])  # Old node indices
-    elem = elem - 1
+    if not isinstance(elem, list):
+        elem = elem - 1
+        zeroindex = np.where(elem < 0)
+    else:
+        elem = [e - 1 for e in elem]
 
     if not isinstance(elem, list):
         idx = np.setdiff1d(oid, elem.ravel(order="F"))  # Indices of isolated nodes
@@ -511,12 +517,16 @@ def removeisolatednode(node, elem, face=None):
 
     if not isinstance(elem, list):
         el = oid[elem]  # Update element list with new indices
+        if isinstance(zeroindex, tuple) and len(zeroindex[0]) > 0:
+            el[zeroindex] = elem[zeroindex]
     else:
         el = [oid[e] for e in elem]
 
     if face is not None:
+        zeroindex = np.where(face == 0)
         if not isinstance(face, list):
             fa = oid[face - 1]  # Update face list with new indices
+            fa[zeroindex] = face[zeroindex] - 1
         else:
             fa = [oid[f - 1] for f in face]
         fa = fa + 1
@@ -866,7 +876,7 @@ def mergesurf(node, elem, *args):
     return newnode, newelem
 
 
-def surfboolean(node, elem, *varargin):
+def surfboolean(node, elem, *varargin, **kwargs):
     """
     Perform boolean operations on triangular meshes and resolve intersecting elements.
 
@@ -939,9 +949,12 @@ def surfboolean(node, elem, *varargin):
             if "node1" not in locals():
                 node1 = node
                 elem1 = elem
-                newnode[:, 3] = 1
-                newelem[:, 3] = 1
-            opstr = " --decouple-inin 1 --shells 2"
+                newnode = np.hstack((newnode, np.ones((newnode.shape[0], 1))))
+                newelem = np.hstack((newelem, np.ones((newelem.shape[0], 1))))
+            if kwargs.get("dir", "in"):
+                opstr = " --decouple-inin 1 --shells 2"
+            else:
+                opstr = " --decouple-outout 1 --shells 2"
             saveoff(node1[:, :3], elem1[:, :3], mwpath("pre_decouple1.off"))
             if no.shape[1] != 3:
                 opstr = f"-q --shells {no}"
@@ -1004,7 +1017,10 @@ def surfboolean(node, elem, *varargin):
             [newnode, np.hstack([nnode, np.ones((nnode.shape[0], 1)) * 4])]
         )
     else:
-        newnode, newelem = readoff(mwpath("post_surfbool.off"))
+        if op == "decouple":
+            newnode, newelem = readoff(mwpath("pre_decouple1_fixed.off"))
+        else:
+            newnode, newelem = readoff(mwpath("post_surfbool.off"))
 
     return newnode, newelem