Merge pull request #554 from ClePol/dev

Fallback for unoriented triangles in sample_parc

doc/api/recon_surf.rst

@@ -15,6 +15,7 @@ recon_surf
     map_surf_label
     N4_bias_correct
     paint_cc_into_pred
+    rewrite_oriented_surface
     rewrite_mc_surface
     rotate_sphere
     sample_parc

recon_surf/recon-surf.sh

@@ -687,7 +687,16 @@ for hemi in lh rh; do
   echo "echo " | tee -a $CMDF
   echo "echo \"=================== Creating surfaces $hemi - fix ========================\"" | tee -a $CMDF
   echo "echo " | tee -a $CMDF
-  cmd="recon-all -subject $subject -hemi $hemi -fix -autodetgwstats -white-preaparc -cortex-label -no-isrunning $hiresflag $fsthreads"
+  cmd="recon-all -subject $subject -hemi $hemi -fix -no-isrunning $hiresflag $fsthreads"
+  RunIt "$cmd" $LF $CMDF
+
+  # fix the surfaces if they are corrupt
+  cmd="$python ${binpath}rewrite_oriented_surface.py --file $sdir/$hemi.orig.premesh --backup $sdir/$hemi.orig.premesh.noorient"
+  RunIt "$cmd" $LF $CMDF
+  cmd="$python ${binpath}rewrite_oriented_surface.py --file $sdir/$hemi.orig --backup $sdir/$hemi.orig.noorient"
+  RunIt "$cmd" $LF $CMDF
+
+  cmd="recon-all -subject $subject -hemi $hemi -autodetgwstats -white-preaparc -cortex-label -no-isrunning $hiresflag $fsthreads"
   RunIt "$cmd" $LF $CMDF
   ## copy nofix to orig and inflated for next step
   # -white (don't know how to call this from recon-all as it needs -whiteonly setting and by default it also creates the pial.

recon_surf/rewrite_oriented_surface.py

@@ -0,0 +1,124 @@
+# Copyright 2024 Image Analysis Lab, German Center for Neurodegenerative Diseases (DZNE), Bonn
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+import shutil
+# IMPORTS
+import sys
+import argparse
+from pathlib import Path
+
+import lapy
+
+__version__ = "1.0"
+
+
+def make_parser() -> argparse.ArgumentParser:
+    """
+    Create a command line interface and return command line options.
+
+    Returns
+    -------
+    options
+        Namespace object holding options.
+    """
+    parser = argparse.ArgumentParser(
+        description="Script to load and safe surface (that are guaranteed to be "
+                    "correctly oriented) under a given name",
+    )
+    parser.add_argument(
+        "--file", "-f",
+        type=Path,
+        dest="file",
+        help="path to surface to check and fix",
+        required=True,
+    )
+    parser.add_argument(
+        "--backup",
+        type=Path,
+        dest="backup",
+        help="if the surface is corrupted, create a backup of the original surface. "
+             "Default: no backup.",
+        default=None,
+    )
+    parser.add_argument(
+        "--version",
+        action="version",
+        version=f"{__version__} 2024/08/08 12:20:10 kueglerd",
+    )
+    return parser
+
+
+def resafe_surface(
+        surface_file: Path | str,
+        surface_backup: Path | str | None = None,
+) -> bool:
+    """
+    Take path to surface_file and rewrite it to fix improperly oriented triangles.
+
+    If the surface is not oriented and surface_backup is set, rename the old
+    surface_file to surface_backup. Else just overwrite with the corrected surface.
+
+    Parameters
+    ----------
+    surface_file : Path, str
+        Path and name of input surface.
+    surface_backup : Path, str, optional
+        Path and name of output surface.
+
+    Returns
+    -------
+    bool
+        Whether the surface was rewritten.
+    """
+    import getpass
+    try:
+        getpass.getuser()
+    except Exception:
+        # nibabel crashes in write_geometry, if getpass.getuser does not return
+        # make sure the process has a username
+        from os import environ
+        environ.setdefault("USERNAME", "UNKNOWN")
+
+    triamesh = lapy.TriaMesh.read_fssurf(str(surface_file))
+    fsinfo = triamesh.fsinfo
+
+    # make sure the triangles are oriented (normals pointing to the same direction
+    if not triamesh.is_oriented():
+        if surface_backup is not None:
+            print(f"Renaming {surface_file} to {surface_backup}")
+            shutil.move(surface_file, surface_backup)
+
+        print("Surface was not oriented, flipping corrupted normals.")
+        triamesh.orient_()
+
+        from packaging.version import Version
+        if Version(lapy.__version__) <= Version("1.0.1"):
+            print(f"lapy version {lapy.__version__}<=1.0.1 detected, fixing fsinfo.")
+            triamesh.fsinfo = fsinfo
+
+        triamesh.write_fssurf(str(surface_file))
+        return True
+    else:
+        print("Surface was already oriented.")
+        return False
+
+
+if __name__ == "__main__":
+    # Command Line options are error checking done here
+    parser = make_parser()
+    args = parser.parse_args()
+
+    print(f"Reading in surface: {args.file} ...")
+    if resafe_surface(args.file, args.backup):
+        print(f"Outputting surface as: {args.file}")
+    sys.exit(0)

recon_surf/sample_parc.py

@@ -295,6 +295,12 @@ def sample_img(surf, img, cortex=None, projmm=0.0, radius=None):
     data = np.asarray(img.dataobj)
     # Use LaPy TriaMesh for vertex normal computation
     T = TriaMesh(surf[0], surf[1])
+
+    # make sure the triangles are oriented (normals pointing to the same direction
+    if not T.is_oriented():
+        print("WARNING: Surface is not oriented, flipping corrupted normals.")
+        T.orient_()
+
     # compute sample coordinates projmm mm along the surface normal
     # in surface RAS coordiante system:
     x = T.v + projmm * T.vertex_normals()