Use neuroreg for registration, mapping, conversion in several places (HypVINN, longitudinal primarily)

HypVINN/README.md

@@ -7,7 +7,7 @@ Hypothalamic subfields segmentation pipeline
 
 ### Requirements
 * Same as FastSurfer.
-* If the T1w and T2w images are available and not co-registered, FreeSurfer should be sourced to run the registration code, and the mri_coreg and mri_vol2vol binaries should also be available.
+* If the T1w and T2w images are available and not co-registered, HypVINN can register the T2w image to the T1w reference internally.
 
 ### Model weights
 * EUDAT (FZ Jülich) data repository: https://b2share.fz-juelich.de/records/2af6da63d5c1414b832c1f606bbd068a
@@ -30,7 +30,7 @@ Note: These weights (version 1.1) are retrained compared to paper ([version 1.0]
  *  `--t2 </dir/T2**.nii.gz>` : T2 image path
  * `--seg_log` :  Path to file in which run logs will be saved. If not set logs will be stored in `/sd/sid/scripts/hypvinn_seg.log` 
 ### Image processing options
- * `--reg_mode` :  Ignored, if no T2 image is passed. Specifies the registration method used to register T1 and T2 images. Options are 'coreg' (default) for mri_coreg, 'robust' for mri_robust_register, and 'none' to skip registration (this requires T1 and T2 are externally co-registered).
+ * `--reg_mode` :  Ignored, if no T2 image is passed. Specifies the registration method used to register T1 and T2 images. Options are 'coreg' (default) for `neuroreg.coreg` and 'none' to skip registration (this requires T1 and T2 are externally co-registered).
  * `--qc_snap`: Activate the creation of QC snapshots of the predicted HypVINN segmentation in `/sd/sid/qc_snapshots`. The created QC snapshots are created to simplify the visual quality control process.
 ###  FastSurfer Technical parameters (see FastSurfer documentation)
  * `--device`
@@ -104,7 +104,7 @@ Note: These weights (version 1.1) are retrained compared to paper ([version 1.0]
             |--hypothalamus.HypVINN.nii.gz(Hypothalamus Segmentation)
             |-- hypothalamus_mask.HypVINN.nii.gz (Hypothalamus Segmentation Mask)
             |-- transforms
-                |-- t2tot1.lta (FreeSurfer registration file, only available if registration is performed)
+                |-- t2tot1.lta (T2-to-T1 registration transform, only available if registration is performed)
         |-- qc_snapshots : QC outputs (optional)
             |-- hypothalamus.HypVINN_qc_screenshoot.png (Coronal quality control image)
         |-- stats : Statistics outputs                                                 

HypVINN/run_prediction.py

@@ -109,11 +109,11 @@ def option_parse() -> argparse.ArgumentParser:
         "--reg_mode",
         type=str,
         default="coreg",
-        choices=["none", "coreg", "robust"],
-        help="Freesurfer Registration type to run. coreg: mri_coreg, "
-             "robust : mri_robust_register, none: entirely deactivates "
-             "registration of T2 to T1, if both images are passed, "
-             "images need to be register properly externally.",
+        choices=["none", "coreg"],
+        help="Registration type to run for T2-to-T1 alignment when both images "
+             "are passed. coreg: neuroreg.coreg, none: entirely deactivates "
+             "registration of T2 to T1, so images need to be registered "
+             "properly externally.",
     )
 
     parser.add_argument(
@@ -156,7 +156,7 @@ def main(
         hypo_maskfile: str = HYPVINN_MASK_NAME,
         qc_snapshots: bool = False,
         threads: int | None = None,
-        reg_mode: Literal["coreg", "robust", "none"] = "coreg",
+        reg_mode: Literal["coreg", "none"] = "coreg",
         batch_size: int = 1,
         async_io: bool = False,
         device: str = "auto",
@@ -195,7 +195,7 @@ def main(
         Whether to create QC snapshots.
     threads : int, optional
         If not None, updates the FastSurfer global setting in `FastSurfer.utils.parallel`.
-    reg_mode : "coreg", "robust", "none", default="coreg"
+    reg_mode : "coreg", "none", default="coreg"
         The registration mode to use.
     batch_size : int, default=1
         The batch size to use.

HypVINN/utils/__init__.py

@@ -14,4 +14,4 @@ class ViewOperationDefinition(TypedDict):
 ViewOperations = dict[Plane, ViewOperationDefinition | None]
 ModalityMode = Literal["t1", "t2", "t1t2"]
 ModalityDict = dict[Literal["t1", "t2"], ndarray]
-RegistrationMode = Literal["robust", "coreg", "none"]
+RegistrationMode = Literal["coreg", "none"]

HypVINN/utils/preproc.py

@@ -36,7 +36,7 @@ def t1_to_t2_registration(
         threads: int = -1,
 ) -> Path:
     """
-    Register T1 to T2 images using either mri_coreg or mri_robust_register.
+    Register the T2 image to the T1 reference using neuroreg.
 
     Parameters
     ----------
@@ -49,7 +49,7 @@ def t1_to_t2_registration(
     lta_path : Path
         The path to the lta transform.
     registration_type : RegistrationMode, default="coreg"
-        The type of registration to be used. It can be either "coreg" or "robust".
+        The type of registration to be used. It can be either "coreg" or "none".
     threads : int, default=-1
         The number of threads to be used. If it is less than or equal to 0, the number
         of threads will be automatically determined.
@@ -62,49 +62,21 @@ def t1_to_t2_registration(
     Raises
     ------
     RuntimeError
-        If mri_coreg, mri_vol2vol, or mri_robust_register fails to run or if they cannot
-        be found.
+        If the requested neuroreg registration backend fails.
     """
-    import shutil
-
-    from FastSurferCNN.utils.parallel import get_num_threads
-    from FastSurferCNN.utils.run_tools import Popen
-
-    if threads <= 0:
-        threads = get_num_threads()
-
-    def from_freesurfer_home(fs_binary: str) -> str:
-        if not os.environ.get("FREESURFER_HOME", ""):
-            raise RuntimeError(
-                f"Could not find {fs_binary}, source FreeSurfer or set the FREESURFER_HOME environment variable"
-            )
-        return os.environ["FREESURFER_HOME"] + "/bin/" + fs_binary
-
-    def run_fs_binary(fs_binary: str, args: list[str]) -> int:
-        fs_binary = shutil.which(fs_binary) or from_freesurfer_home(fs_binary)
-        args = [fs_binary] + list(map(str, args))
-        LOGGER.info("Running " + " ".join(args))
-        retval = Popen(args).finish()
-        if retval.retcode != 0:
-            LOGGER.error(f"{fs_binary} failed with error code {retval.retcode}.")
-            raise RuntimeError(f"{fs_binary} failed")
-
-        LOGGER.info(f"{fs_binary} finished in {retval.runtime}!")
+    from neuroreg import coreg
 
     if registration_type == "coreg":
-        run_fs_binary(
-            "mri_coreg",
-            ["--mov", t2_path, "--targ", t1_path, "--reg", lta_path, "--threads", str(threads)],
-        )
-        run_fs_binary(
-            "mri_vol2vol",
-            ["--mov", t2_path, "--targ", t1_path, "--reg", lta_path, "--o", output_path, "--cubic", "--keep-precision"],
+        LOGGER.info("Running neuroreg.coreg for T2-to-T1 registration.")
+        coreg(
+            str(t2_path),
+            str(t1_path),
+            lta_name=str(lta_path),
+            mapped_name=str(output_path),
+            keep_dtype=True,
         )
     else:
-        run_fs_binary(
-            "mri_robust_register",
-            ["--mov", t2_path, "--dst", t1_path, "--lta", lta_path, "--mapmov", output_path, "--cost", "NMI"],
-        )
+        raise ValueError(f"Unknown registration type: {registration_type}")
 
     return output_path
 
@@ -125,7 +97,8 @@ def hypvinn_preproc(
     mode : ModalityMode
         The mode for HypVINN. It should be "t1t2".
     reg_mode : RegistrationMode
-        The registration mode. If it is not "none", the function will register T1 to T2 images.
+        The registration mode. If it is not "none", the function will register the
+        T2 image to the T1 reference.
     t1_path : Path
         The path to the T1 image.
     t2_path : Path
@@ -162,7 +135,7 @@ def hypvinn_preproc(
                 f"T2 image will be interpolated to the resolution of the T1 image."
             )
 
-        LOGGER.info("Registering T1 to T2 ...")
+        LOGGER.info("Registering T2 to T1 ...")
         t1_to_t2_registration(
             t1_path=t1_path,
             t2_path=t2_path,

pyproject.toml

@@ -38,7 +38,7 @@ dependencies = [
     'matplotlib>=3.7.1',
     'meshpy>=2025.1.1', # needed for FastSurfer-CC
     'monai>=1.4.0', # needed for FastSurfer-CC
-    'neuroreg>=0.6.2', # needed for registration / etiv
+    'neuroreg>=0.10.1', # needed for registration / etiv
     'nibabel>=5.4.0', # needed to fix a bug in nibabel
     'numpy>=1.25',
     'packaging',

recon_surf/long_prepare_template.sh

@@ -24,7 +24,7 @@
 #  to setup a base template for processing with a slightly modified version of run_fastsurfer.
 # 
 #  1. It obtains a brainmask on all time points for a specific subject
-#  2. It uses mri_robust_template to co-register all time points into a mid-space
+#  2. It uses neuroreg's multireg to co-register all time points into a mid-space
 #     and sets up the base (subject template) directory with the aligned images the median image
 #     and the forward and backward rigid transformations
 # 
@@ -38,8 +38,7 @@
 #   - Maybe use intersection of tp masks as brainmask for base, as done in FreeSurfer.
 #   - Add flag for adding a new time point to an existing base/template.
 #
-#  FreeSurfer requirements: 
-#  mri_convert, mri_robust_template, mri_mask, mri_concatenate_lta, make_upright
+#  neuroreg requirements: multireg, segreg, vol2vol, lta, mri (mask, diff, info)
 #
 ###################################################################################################
 
@@ -68,6 +67,7 @@ robust_template_avg_arg=1  # median for template creation (if more than 1 time p
 # default arguments
 python="python3 -s" # avoid user-directory package inclusion
 sd="$SUBJECTS_DIR"
+device="" # torch device for neuroreg (empty/auto = let neuroreg pick GPU else CPU)
 
 # init variables that need to be passed
 tid=""
@@ -163,7 +163,9 @@ case $key in
   --t1s) while [[ $# -gt 0 ]] && [[ $1 != -* ]] ; do t1s+=("$1") ; shift ; done ;;
   --sd) export SUBJECTS_DIR="$1" ; shift  ;;
   # these flags are passed through to run_prediction.py
-  --vox_size|--device|--viewagg_device|--conform_to_1mm_threshold) run_pred_flags+=("$key" "$1") ; shift ;;
+  --vox_size|--viewagg_device|--conform_to_1mm_threshold) run_pred_flags+=("$key" "$1") ; shift ;;
+  # device is passed through to run_prediction.py and also forwarded to neuroreg
+  --device) run_pred_flags+=("$key" "$1") ; device="$1" ; shift ;;
   --threads|--threads_seg) run_pred_flags+=("--threads" "$1") ; shift ;;
   --batch) run_pred_flags+=("--batch_size" "$1") ; shift ;;
   # these known arguments get ignored
@@ -226,6 +228,8 @@ fi
 # check that SUBJECTS_DIR exists
 check_create_subjects_dir_properties "$SUBJECTS_DIR"
 
+auto_detect_fs_license "the longitudinal template preparation" || exit $?
+
 ################################## SETUP and LOGFILE ##############################
 
 
@@ -283,7 +287,7 @@ do
   # check if geometry differs across time
   if [[ "$s" != "${t1s[0]}" ]]
   then
-    cmda=(mri_diff --notallow-pix --notallow-geo "$s" "${t1s[0]}" --res-thresh "0.000001")
+    cmda=($python -m neuroreg.cli.mri diff --notallow-pix --notallow-geo "$s" "${t1s[0]}" --res-thresh "0.000001")
     difftext=$("${cmda[@]}")
     retcode=${PIPESTATUS[0]}
     if [[ "$retcode" != 0 ]] ; then geodiff+="Comparing $s and ${t1s[0]} (code $retcode):\n$difftext\n" ; fi
@@ -296,7 +300,7 @@ then
     echo "*******************************************************************************" 
     echo "WARNING: Image parameters differ across time, maybe due to acquisition changes?"
     echo "         Consistent changes in, e.g., resolution can potentially bias a "
-    echo "         longitudinal study! You can check image parameters by running mri_info"
+    echo "         longitudinal study! You can check image parameters by running neuroreg mri info"
     echo "         on each input image."
     echo "*******************************************************************************"
     echo "$geodiff"
@@ -331,7 +335,7 @@ for ((i=0;i<${#tpids[@]};++i)); do
   mkdir -p "$mdir"
   # Import (copy) raw inputs (convert to extension format)
   t1input=$mdir/cross_input${extension}
-  cmd="mri_convert ${t1s[i]} $t1input"
+  cmd="$python -m neuroreg.cli.vol2vol --in ${t1s[i]} --out $t1input"
   RunIt "$cmd" "$LF"
 
   # conform !!!!!!! should we conform to some common value, determined from all time points?? !!!!!!
@@ -360,7 +364,7 @@ for ((i=0;i<${#tpids[@]};++i)); do
   run_it "$LF" "${cmda[@]}"
 
   # mask is binary, we need to use on conformed image:
-  cmda=(mri_mask "$conformed_name" "$mask_name" "$mdir/cross_brainmask${extension}")
+  cmda=($python -m neuroreg.cli.mri mask "$conformed_name" "$mask_name" "$mdir/cross_brainmask${extension}")
   run_it "$LF" "${cmda[@]}"
 done
 
@@ -401,38 +405,60 @@ done
 
 if [ ${#tpids[@]} == 1 ]
 then
-  # if only a single time point, create fake 'base' by making the image upright
-  # this assures that also subjects with a single time point get processed as the other
-  # subjects in the longitudinal stream
+  # If only a single time point, we still create a 'base' so that single-tp subjects are
+  # processed like the others. Instead of FreeSurfer's make_upright (mirror-symmetry
+  # uprighting), we rigidly register this time point's segmentation to the reference centroid
+  # template - the same way the talairach initialization prealigns - which yields a
+  # reproducible, template-aligned (upright) base pose without scaling the anatomy. The
+  # segmentation always exists here (it was created above to obtain the brainmask). The base
+  # lives in the time point's own (conformed) geometry.
+
+  seg0="$SUBJECTS_DIR/$tid/long-inputs/${tpids[0]}/cross_aparc+aseg.orig${extension}"
+  reference_centroids="mni_icbm152_t1_tal_nlin_asym_09c"
+
+  # 1. rigid segmentation-to-template registration: tp -> base pose (RAS-to-RAS transform)
+  cmd="$python -m neuroreg.cli.segreg --seg $seg0 --centroids $reference_centroids"
+  cmd="$cmd --dof 6 --lta ${ltaXforms[0]}"
+  RunIt "$cmd" "$LF"
 
-  # 1. make the norm upright (base space)
-  cmd="make_upright ${normInVols[0]} \
-       ${SUBJECTS_DIR}/$tid/mri/base_brainmask${extension} ${ltaXforms[0]}"
+  # 2. create the base brainmask by mapping the norm into the base pose. The transform is
+  #    RAS-to-RAS, so the target geometry is taken from the input itself (--ref), keeping
+  #    the base in the time point geometry.
+  cmd="$python -m neuroreg.cli.vol2vol --in ${normInVols[0]} --transform ${ltaXforms[0]}"
+  cmd="$cmd --ref ${normInVols[0]} --interp cubic"
+  cmd="$cmd --out ${SUBJECTS_DIR}/$tid/mri/base_brainmask${extension}"
   RunIt "$cmd" "$LF"
 
-  # 2. create the upright orig volume
-  cmd="mri_convert -rt cubic \
-       -at ${ltaXforms[0]} ${subjInVols[0]} ${SUBJECTS_DIR}/$tid/mri/orig.mgz"
+  # 3. create the base orig volume the same way
+  cmd="$python -m neuroreg.cli.vol2vol --in ${subjInVols[0]} --transform ${ltaXforms[0]}"
+  cmd="$cmd --ref ${subjInVols[0]} --interp cubic"
+  cmd="$cmd --out ${SUBJECTS_DIR}/$tid/mri/orig.mgz"
   RunIt "$cmd" "$LF"
 
 else #more than 1 time point:
 
+  # forward the device to neuroreg if explicitly requested (empty/auto = neuroreg auto-selects)
+  device_opt=""
+  if [[ -n "$device" ]] && [[ "$device" != "auto" ]] ; then device_opt="--device $device" ; fi
 
-  # create the 'mean/median' norm volume:
-  cmd="mri_robust_template --mov ${normInVols[*]}"
+  # robust co-registration of all time points into an unbiased mid-space, creating the
+  # 'mean/median' norm (brainmask) volume and the forward transforms (tp -> base):
+  cmd="$python -m neuroreg.cli.multireg --mov ${normInVols[*]}"
   cmd="$cmd --lta ${ltaXforms[*]}"
   cmd="$cmd --template ${SUBJECTS_DIR}/$tid/mri/base_brainmask${extension}"
   cmd="$cmd --average ${robust_template_avg_arg}"
-  cmd="$cmd --sat 4.685"
+  cmd="$cmd --sat 4.685 $device_opt"
   RunIt "$cmd" "$LF"
 
-  # create the 'mean/median' input (orig) volume:
-  cmd="mri_robust_template --mov ${subjInVols[*]}"
+  # create the 'mean/median' input (orig) volume by reusing the transforms above (no
+  # registration, --noit); the per-tp LTAs carry the base geometry so orig.mgz lands on
+  # the same grid as base_brainmask:
+  cmd="$python -m neuroreg.cli.multireg --mov ${subjInVols[*]}"
   cmd="$cmd --average ${robust_template_avg_arg}"
   cmd="$cmd --ixforms ${ltaXforms[*]}"
   cmd="$cmd --noit"
   t1=${SUBJECTS_DIR}/$tid/mri/orig.mgz
-  cmd="$cmd --template $t1"
+  cmd="$cmd --template $t1 $device_opt"
   RunIt "$cmd" "$LF"
 
 fi # more than one time point
@@ -451,9 +477,11 @@ done
 for ((i=0;i<${#tpids[@]};++i))
 do
   mdir="$SUBJECTS_DIR/$tid/long-inputs/${tpids[i]}"
-  # map orig to base space
-  cmd="mri_convert -at ${ltaXforms[$i]} -rt $interpol $mdir/cross_input${extension} $mdir/long_conform${extension}"
+  # map orig to base space; the transforms are RAS-to-RAS, so the base geometry is supplied
+  # explicitly via --ref (rather than relying on a target geometry stored in the transform)
+  cmd="$python -m neuroreg.cli.vol2vol --in $mdir/cross_input${extension}"
+  cmd="$cmd --transform ${ltaXforms[$i]} --ref ${SUBJECTS_DIR}/$tid/mri/base_brainmask${extension}"
+  cmd="$cmd --interp $interpol --out $mdir/long_conform${extension}"
   RunIt "$cmd" "$LF"
 done
 
-

recon_surf/recon-surf.sh

@@ -468,7 +468,7 @@ fi
 
 # create orig.mgz and aparc.DKTatlas+aseg.orig.mgz (copy of T1 and segmentation)
 # also ensures .mgz format (in case inputs are nifti)
-cmd="mri_convert $t1 $mdir/orig.mgz"
+cmd="$python -m neuroreg.cli.vol2vol --in $t1 --out $mdir/orig.mgz"
 RunIt "$cmd" "$LF"
 
 asegdkt_segfile_manedit=$(add_file_suffix "$asegdkt_segfile" "manedit")
@@ -477,7 +477,7 @@ if [[ ! "$asegdkt_segfile_manedit" =~ (\.manedit){2,}\. ]] && [[ -f "$asegdkt_se
 then
     asegdkt_segfile="$asegdkt_segfile_manedit" # use the manedit file
 fi
-cmd="mri_convert $asegdkt_segfile $mdir/aparc.DKTatlas+aseg.orig.mgz"
+cmd="$python -m neuroreg.cli.vol2vol --in $asegdkt_segfile --out $mdir/aparc.DKTatlas+aseg.orig.mgz"
 RunIt "$cmd" "$LF"
 
 # link original T1 input to rawavg (needed by pctsurfcon)
@@ -597,7 +597,7 @@ fi
 
 # the difference between nu and orig_nu is the fact that nu has the talairach-registration header
 # create norm by masking nu (supports manedit-ed mask)
-cmda=(mri_mask "$mdir/nu.mgz" "$mask" "$mdir/norm.mgz")
+cmda=($python -m neuroreg.cli.mri mask "$mdir/nu.mgz" "$mask" "$mdir/norm.mgz")
 run_it "$LF" "${cmda[@]}"
 if [[ "$get_t1" == "true" ]]
 then
@@ -618,7 +618,7 @@ then
   cmda=(mri_normalize -g 1 -seed 1234 -mprage "$mdir/nu.mgz" "$mdir/T1.mgz" "${noconform_if_hires[@]}")
   run_it "$LF" "${cmda[@]}"
   # create brainmask by masking T1 (supports manedit-ed mask)
-  cmda=(mri_mask "$mdir/T1.mgz" "$mask" "$mdir/brainmask.mgz")
+  cmda=($python -m neuroreg.cli.mri mask "$mdir/T1.mgz" "$mask" "$mdir/brainmask.mgz")
   run_it "$LF" "${cmda[@]}"
 else
   # create brainmask by linkage to norm.mgz (masked nu.mgz)

recon_surf/talairach-reg.sh

@@ -162,7 +162,7 @@ else
 
     prealigned_name=$mdir/segreg_prealigned.mgz
     target_geom="$FREESURFER_HOME/average/mni305.cor.mgz" # only used for its conformed 1mm header; does not imply registration to this template
-    cmd=(mri_convert --apply_transform "$prealigned_lta" --reslice_like "$target_geom" "$norm_name" "$prealigned_name")
+    cmd=($python -m neuroreg.cli.vol2vol --in "$norm_name" --transform "$prealigned_lta" --ref "$target_geom" --interp linear --out "$prealigned_name")
     run_it "$LF" "${cmd[@]}"
 
     # talairach.xfm: compute talairach full head (25sec)