From 054a56ff97b25b3e7fb0455d955a8f1753028f04 Mon Sep 17 00:00:00 2001 From: Clemens Pollak Date: Wed, 7 Aug 2024 15:54:35 +0200 Subject: [PATCH 1/4] Fallback for unoriented triangles in sample_parc --- recon_surf/sample_parc.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/recon_surf/sample_parc.py b/recon_surf/sample_parc.py index 17461821..3d784ad2 100644 --- a/recon_surf/sample_parc.py +++ b/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() From 0f635ce82481af92bf59ed5494bd8bf6812e6f87 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20K=C3=BCgler?= Date: Wed, 7 Aug 2024 17:37:53 +0200 Subject: [PATCH 2/4] Create a new rewrite_oriented_surface.py script to fix incorrectly oriented surfaces Integrate the script into recon-surf.sh --- doc/api/recon_surf.rst | 1 + recon_surf/recon-surf.sh | 8 ++ recon_surf/rewrite_oriented_surface.py | 102 +++++++++++++++++++++++++ 3 files changed, 111 insertions(+) create mode 100644 recon_surf/rewrite_oriented_surface.py diff --git a/doc/api/recon_surf.rst b/doc/api/recon_surf.rst index 5ece8a11..54403fb9 100644 --- a/doc/api/recon_surf.rst +++ b/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 diff --git a/recon_surf/recon-surf.sh b/recon_surf/recon-surf.sh index 1ede909b..4ef72ba1 100755 --- a/recon_surf/recon-surf.sh +++ b/recon_surf/recon-surf.sh @@ -703,6 +703,14 @@ for hemi in lh rh; do cmd="recon-all -subject $subject -hemi $hemi -smooth2 -inflate2 -curvHK -no-isrunning $hiresflag $fsthreads" RunIt "$cmd" $LF $CMDF + # rename the freesurfer preaparc surface + cmd="mv $sdir/$hemi.white.preaparc $sdir/$hemi.white.preaparc.nofix" + RunIt "$cmd" $LF $CMDF + + # fix the surfaces if they are corrupt + cmd="$python ${binpath}rewrite_oriented_surface.py -i $sdir/$hemi.white.preaparc.nofix -o $sdir/$hemi.white.preaparc" + RunIt "$cmd" $LF $CMDF + # ============================= MAP-DKT ========================================================== diff --git a/recon_surf/rewrite_oriented_surface.py b/recon_surf/rewrite_oriented_surface.py new file mode 100644 index 00000000..381805df --- /dev/null +++ b/recon_surf/rewrite_oriented_surface.py @@ -0,0 +1,102 @@ +# 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. + +# 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( + "--input", "-i", + type=Path, + dest="input_surf", + help="path to input surface", + required=True, + ) + parser.add_argument( + "--output", "-o", + type=Path, + dest="output_surf", + help="path to output surface", + required=True, + ) + parser.add_argument( + "--version", + action="version", + version=f"{__version__} 2024/08/08 12:20:10 kueglerd", + ) + return parser + + +def resafe_surface(insurf: Path | str, outsurf: Path | str) -> None: + """ + Take path to insurf and rewrite it to outsurf thereby fixing improperly oriented triangles. + + Parameters + ---------- + insurf : Path, str + Path and name of input surface. + outsurf : Path, str + Path and name of output surface. + """ + 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(insurf)) + + # make sure the triangles are oriented (normals pointing to the same direction + if not triamesh.is_oriented(): + print("Surface was not oriented, flipping corrupted normals.") + triamesh.orient_() + else: + print("Surface was oriented.") + + triamesh.write_fssurf(str(outsurf)) + + +if __name__ == "__main__": + # Command Line options are error checking done here + parser = make_parser() + args = parser.parse_args() + surf_in = args.input_surf + surf_out = args.output_surf + + print(f"Reading in surface: {surf_in} ...") + resafe_surface(surf_in, surf_out) + print(f"Outputting surface as: {surf_out}") + sys.exit(0) From 170435ede8ae216950059a8f8e4b196555c68031 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20K=C3=BCgler?= Date: Thu, 8 Aug 2024 18:57:31 +0200 Subject: [PATCH 3/4] Move the surface fixing up and apply it to rh/lh.orig.premesh --- recon_surf/recon-surf.sh | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/recon_surf/recon-surf.sh b/recon_surf/recon-surf.sh index 4ef72ba1..2118a078 100755 --- a/recon_surf/recon-surf.sh +++ b/recon_surf/recon-surf.sh @@ -687,7 +687,18 @@ 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 + + # rename the freesurfer preaparc surface + cmd="mv $sdir/$hemi.orig.premesh $sdir/$hemi.orig.premesh.noorient" + RunIt "$cmd" $LF $CMDF + + # fix the surfaces if they are corrupt + cmd="$python ${binpath}rewrite_oriented_surface.py -i $sdir/$hemi.orig.premesh.noorient -o $sdir/$hemi.orig.premesh" + 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. @@ -703,14 +714,6 @@ for hemi in lh rh; do cmd="recon-all -subject $subject -hemi $hemi -smooth2 -inflate2 -curvHK -no-isrunning $hiresflag $fsthreads" RunIt "$cmd" $LF $CMDF - # rename the freesurfer preaparc surface - cmd="mv $sdir/$hemi.white.preaparc $sdir/$hemi.white.preaparc.nofix" - RunIt "$cmd" $LF $CMDF - - # fix the surfaces if they are corrupt - cmd="$python ${binpath}rewrite_oriented_surface.py -i $sdir/$hemi.white.preaparc.nofix -o $sdir/$hemi.white.preaparc" - RunIt "$cmd" $LF $CMDF - # ============================= MAP-DKT ========================================================== From e818b443aa6ed5c5801054dc55e1800e8fa07039 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20K=C3=BCgler?= Date: Mon, 12 Aug 2024 11:51:39 +0200 Subject: [PATCH 4/4] Add version check and recover fsinfo, if not using the fixed lapy version. Update the names of the command line arguments (--file and --backup) To reduce unnecessary file duplication, only save noorient, if the surface was actually not oriented. --- recon_surf/recon-surf.sh | 8 ++-- recon_surf/rewrite_oriented_surface.py | 64 +++++++++++++++++--------- 2 files changed, 46 insertions(+), 26 deletions(-) diff --git a/recon_surf/recon-surf.sh b/recon_surf/recon-surf.sh index 2118a078..33a24200 100755 --- a/recon_surf/recon-surf.sh +++ b/recon_surf/recon-surf.sh @@ -690,12 +690,10 @@ for hemi in lh rh; do cmd="recon-all -subject $subject -hemi $hemi -fix -no-isrunning $hiresflag $fsthreads" RunIt "$cmd" $LF $CMDF - # rename the freesurfer preaparc surface - cmd="mv $sdir/$hemi.orig.premesh $sdir/$hemi.orig.premesh.noorient" - RunIt "$cmd" $LF $CMDF - # fix the surfaces if they are corrupt - cmd="$python ${binpath}rewrite_oriented_surface.py -i $sdir/$hemi.orig.premesh.noorient -o $sdir/$hemi.orig.premesh" + 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" diff --git a/recon_surf/rewrite_oriented_surface.py b/recon_surf/rewrite_oriented_surface.py index 381805df..fc9c4a37 100644 --- a/recon_surf/rewrite_oriented_surface.py +++ b/recon_surf/rewrite_oriented_surface.py @@ -11,7 +11,7 @@ # 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 @@ -36,18 +36,19 @@ def make_parser() -> argparse.ArgumentParser: "correctly oriented) under a given name", ) parser.add_argument( - "--input", "-i", + "--file", "-f", type=Path, - dest="input_surf", - help="path to input surface", + dest="file", + help="path to surface to check and fix", required=True, ) parser.add_argument( - "--output", "-o", + "--backup", type=Path, - dest="output_surf", - help="path to output surface", - required=True, + dest="backup", + help="if the surface is corrupted, create a backup of the original surface. " + "Default: no backup.", + default=None, ) parser.add_argument( "--version", @@ -57,16 +58,27 @@ def make_parser() -> argparse.ArgumentParser: return parser -def resafe_surface(insurf: Path | str, outsurf: Path | str) -> None: +def resafe_surface( + surface_file: Path | str, + surface_backup: Path | str | None = None, +) -> bool: """ - Take path to insurf and rewrite it to outsurf thereby fixing improperly oriented triangles. + 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 ---------- - insurf : Path, str + surface_file : Path, str Path and name of input surface. - outsurf : Path, str + surface_backup : Path, str, optional Path and name of output surface. + + Returns + ------- + bool + Whether the surface was rewritten. """ import getpass try: @@ -77,26 +89,36 @@ def resafe_surface(insurf: Path | str, outsurf: Path | str) -> None: from os import environ environ.setdefault("USERNAME", "UNKNOWN") - triamesh = lapy.TriaMesh.read_fssurf(str(insurf)) + 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_() - else: - print("Surface was oriented.") - triamesh.write_fssurf(str(outsurf)) + 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() - surf_in = args.input_surf - surf_out = args.output_surf - print(f"Reading in surface: {surf_in} ...") - resafe_surface(surf_in, surf_out) - print(f"Outputting surface as: {surf_out}") + print(f"Reading in surface: {args.file} ...") + if resafe_surface(args.file, args.backup): + print(f"Outputting surface as: {args.file}") sys.exit(0)