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 <t2_path>` : T2 image path (optional)
  * `--seg_log` :  Path to file in which run logs will be saved. If not set logs will be stored in `<subject_dir>/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 `<subject_dir>/qc_snapshots`. The created QC snapshots are created to simplify the visual quality control process.
 ###  FastSurfer Technical parameters (see FastSurfer documentation)
  * `--device`
@@ -116,7 +116,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(
@@ -170,7 +170,7 @@ def main(
         hypo_statsfile: str = HYPVINN_STATS_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",
@@ -211,7 +211,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.
@@ -229,6 +229,10 @@ def main(
     """
     if threads is not None and threads > 1:
         set_num_threads(threads)
+    # here on the main thread, not in Inference, because the T2 registration below is submitted to
+    # a worker before Inference exists and neuroreg does no thread handling of its own. Without
+    # this it would run on torch's own default, ignoring --threads.
+    torch.set_num_threads(get_num_threads())
 
     from concurrent.futures import Future
 

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 the T2 image to the T1 reference 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.
@@ -61,50 +61,25 @@ 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.
+    ValueError
+        If `registration_type` names a mode this function does not implement.
     """
-    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.")
+        # on the CPU, coreg's default. Whether a GPU would be faster here has not been measured.
+        # The thread count is the process-wide torch one, set in run_prediction before this is
+        # submitted, since neuroreg does no thread handling of its own.
+        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
 

doc/overview/LONG.md

@@ -58,7 +58,7 @@ $FASTSURFER_HOME/long_fastsurfer.sh \
 
 Note, with a few exceptions, you can add additional flags that can be understood by `run_fastsurfer.sh`, which will be passed through, e.g. the `--3T` when working with 3T images. The exceptions: `--sid`, `--t1` and `--t2` are set per time point from `--tpids`, `--t1s` and `--t2s`; `--seg_only` and `--surf_only` are replaced by stages (see [Running Individual Stages](#running-individual-stages) below); and `--remove_suffix`, `--keepgeom` and `--native_image` are not supported in longitudinal processing.
 
-If you have T2-weighted images, pass them with `--t2s <t2_path_1> <t2_path_2> ...`, one per time point in the same order as `--tpids`. They are used by the hypothalamus module (HypVINN) in the segmentation of each time point, where each T2 is registered directly to that time point's T1 in template space. The template itself does not need them, as the hypothalamus module does not run on it. Either every time point gets a T2 or none does: a T2 changes what the hypothalamus module computes, so a series where only some time points had one would not be comparable over time. `--reg_mode` selects how each T2 reaches template space: with `coreg` (the default) or `robust`, it is registered to its time point's T1 there; with `none`, it must already be co-registered with the T1 you passed in `--t1s`, and it is then mapped into template space with the same transform as that T1.
+If you have T2-weighted images, pass them with `--t2s <t2_path_1> <t2_path_2> ...`, one per time point in the same order as `--tpids`. They are used by the hypothalamus module (HypVINN) in the segmentation of each time point, where each T2 is registered directly to that time point's T1 in template space. The template itself does not need them, as the hypothalamus module does not run on it. Either every time point gets a T2 or none does: a T2 changes what the hypothalamus module computes, so a series where only some time points had one would not be comparable over time. `--reg_mode` selects how each T2 reaches template space: with `coreg` (the default), it is registered to its time point's T1 there; with `none`, it must already be co-registered with the T1 you passed in `--t1s`, and it is then mapped into template space with the same transform as that T1.
 
 The above command will, of course, be slightly different when using your preferred installation in Singularity or Docker. For example, for Singularity:
 
@@ -161,7 +161,7 @@ Behind the Scenes
      --t1s <t1_path_1> <t1_path_2> ... \
      --tpids <tpid_1> <tpid_2>
    ```
-   This segments each time point once, only to obtain a brain mask, and then registers (aligns) all time point images into the unbiased mid-space using `mri_robust_template`. It will also create the template image, the voxel-wise median across time. For single time point cases, it will align the input into a standard upright position. Finally, it maps each time point's input into template space, where the time point stages below pick it up.
+   This segments each time point once, only to obtain a brain mask, and then registers (aligns) all time point images into the unbiased mid-space using neuroreg's robust multi-image registration. It will also create the template image, the voxel-wise median across time. For single time point cases, it will align the input into a standard upright position, obtained from a rigid registration of its segmentation to MNI centroids. Finally, it maps each time point's input into template space, where the time point stages below pick it up.
 2. **Template Seg** (stage `template_seg`): Next, the template image will be segmented via a call to `run_fastsurfer.sh --sid <template_id> --base --seg_only ...` where the `--base` flag indicates that the input image will be taken from the already existing template directory. The cerebellum and hypothalamus modules do not run on the template.
 3. **Template Surf** (stage `template_surf`): This is followed by the surface processing of the template `run_fastsurfer.sh --sid <template_id> --base --surf_only ...`.
 4. **Long Seg** (stage `long_seg`): Next, the segmentation of each time point is performed `run_fastsurfer.sh --sid <tpid_n> --long <template_id> --seg_only ...`. It only needs step 1, so it can run at the same time as steps 2 and 3.

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.12.0', # needed for registration / etiv
     'nibabel>=5.4.0', # needed to fix a bug in nibabel
     'numpy>=1.25',
     'packaging',