Transfert functions for reorient

.github/workflows/test-ml.yml

@@ -53,7 +53,8 @@ jobs:
             libfreetype6-dev \
             libdrm-dev \
             libgl1-mesa-dev \
-            libosmesa6-dev
+            libosmesa6-dev \
+            python3-dev \
 
       - name: stdlib checkout
         if: ${{ !contains(steps.python-selector.outputs.python-version, '3.12') }}

.github/workflows/test.yml

@@ -52,7 +52,8 @@ jobs:
             libfreetype6-dev \
             libdrm-dev \
             libgl1-mesa-dev \
-            libosmesa6-dev
+            libosmesa6-dev \
+            python3-dev \
 
       - name: stdlib checkout
         if: ${{ !contains(steps.python-selector.outputs.python-version, '3.12') }}

.github/workflows/test_tutorials.yml

@@ -50,8 +50,9 @@ jobs:
             libfreetype6-dev \
             libdrm-dev \
             libgl1-mesa-dev \
-           libosmesa6-dev \
-           wget
+            libosmesa6-dev \
+            python3-dev \
+            wget
 
       - name: stdlib checkout
         if: ${{ !contains(steps.python-selector.outputs.python-version, '3.12') }}

pyproject.toml

@@ -52,6 +52,7 @@ dependencies = [
 "numpy==1.26.*",
 "openpyxl==3.1.*",
 "packaging==24.*",
+"pathspec==0.12.*",
 "pybids==0.18.*",
 "PyMCubes==0.1.*",
 "pyparsing==3.2.*",

src/scilpy/cli/scil_bundle_explore_bundleseg.py

@@ -478,11 +478,16 @@ def main():
     offset = 0
     count = 0
     for bundle in mapping.keys():
-        filename = glob.glob(f'{os.path.join(args.in_folder, bundle)}.t?k')[0]
-
-        if not os.path.exists(filename):
-            logging.warning(f'File {filename} not found.')
+        files = glob.glob(f'{os.path.join(args.in_folder, bundle)}.t?k')
+        if len(files) == 0:
+            logging.warning("Could not find any file fitting pattern {}"
+                            .format(os.path.join(args.in_folder, bundle)))
             continue
+        elif len(files) > 1:
+            logging.warning("Found two files for bundle {}. Selecting the "
+                            "first one. Verify your files!".format(bundle))
+
+        filename = files[0]
         count += 1
 
         tmp_sft = load_tractogram(filename, ref_img)

src/scilpy/cli/scil_tractogram_compute_density_map.py

@@ -77,13 +77,15 @@ def main():
     transformation, dimensions, _, _ = sft.space_attributes
 
     # Processing
+    logging.info("Computing density map...")
     if args.endpoints_only:
         streamline_count = get_endpoints_density_map(sft)
     else:
         streamline_count = compute_tract_counts_map(sft.streamlines,
                                                     dimensions)
 
     # Saving
+    logging.info("Saving density map {}".format(args.out_img))
     dtype_to_use = np.int32
     if args.binary is not None:
         if args.binary == 1:

src/scilpy/cli/scil_tractogram_segment_with_bundleseg.py

@@ -55,7 +55,8 @@
                              add_reference_arg, add_verbose_arg,
                              assert_inputs_exist,
                              assert_output_dirs_exist_and_empty,
-                             load_matrix_in_any_format, ranged_type)
+                             load_matrix_in_any_format, ranged_type,
+                             assert_inputs_dirs_exist)
 from scilpy.segment.voting_scheme import VotingScheme
 from scilpy.version import version_string
 
@@ -136,6 +137,7 @@ def main():
     logging.getLogger().setLevel(logging.getLevelName('INFO'))
 
     # Verifications
+    assert_inputs_dirs_exist(parser, args.in_directory)
     in_models_directories = [
         os.path.join(args.in_directory, x)
         for x in os.listdir(args.in_directory)

src/scilpy/cli/scil_tractogram_segment_with_recobundles.py

@@ -112,8 +112,8 @@ def main():
     if args.tractogram_clustering_thr and args.in_pickle:
         parser.error("Option --tractogram_clustering_thr should not be "
                      "used with --in_pickle.")
-    else:
-        # Setting default value. (Will be ignored in args.in_pickle)
+    elif args.tractogram_clustering_thr is not None:
+        # Setting default value. (Will be ignored if args.in_pickle)
         args.tractogram_clustering_thr = 8.0
 
     # Loading

src/scilpy/image/volume_math.py

@@ -55,6 +55,7 @@ def get_array_ops():
         ('subtraction', subtraction),
         ('multiplication', multiplication),
         ('division', division),
+        ('maximum', maximum),
         ('mean', mean),
         ('std', std),
         ('correlation', neighborhood_correlation),
@@ -442,6 +443,27 @@ def convert(input_list, ref_img):
     return input_list[0].get_fdata(dtype=np.float64)
 
 
+def maximum(input_list, ref_img):
+    """
+    maximum: IMGs
+        Compute the voxel-wise maximum across images.
+    """
+    _validate_length(input_list, 2, at_least=True)
+    _validate_imgs_type(*input_list, all_imgs=False)
+    _validate_same_shape(*input_list, ref_img, all_imgs=False)
+
+    output_data = np.zeros(ref_img.header.get_data_shape(), dtype=np.float64)
+    for img in input_list:
+        if isinstance(img, nib.Nifti1Image):
+            data = img.get_fdata(dtype=np.float64)
+            output_data = np.maximum(output_data, data)
+            img.uncache()
+        else:
+            output_data = np.maximum(output_data, img)
+
+    return output_data
+
+
 def addition(input_list, ref_img):
     """
     addition: IMGs

src/scilpy/io/image.py

@@ -8,7 +8,6 @@
 
 from scilpy.utils import is_float
 
-
 def load_img(arg):
     """
     Function to create the variable for scil_volume_math main function.

src/scilpy/io/stateful_image.py

@@ -0,0 +1,182 @@
+# -*- coding: utf-8 -*-
+
+import nibabel as nib
+import numpy as np
+from dipy.io.utils import get_reference_info
+
+class StatefulImage(nib.Nifti1Image):
+    """
+    A class that extends nib.Nifti1Image to manage image orientation state.
+
+    This class ensures that image data loaded into memory is always in a
+    consistent orientation (RAS by default), while preserving the original
+    on-disk orientation information. When saving, the image is automatically
+    reverted to its original orientation, ensuring non-destructive operations.
+    """
+
+    def __init__(self, dataobj, affine, header=None, extra=None,
+                 file_map=None, original_affine=None,
+                 original_dimensions=None, original_voxel_sizes=None,
+                 original_axcodes=None):
+        """
+        Initialize a StatefulImage object.
+
+        Extends the Nifti1Image constructor to store original orientation info.
+        """
+        super().__init__(dataobj, affine, header, extra, file_map)
+
+        # Store original image information
+        self._original_affine = original_affine
+        self._original_dimensions = original_dimensions
+        self._original_voxel_sizes = original_voxel_sizes
+        self._original_axcodes = original_axcodes
+
+    @classmethod
+    def load(cls, filename, to_orientation='RAS'):
+        """
+        Load a NIfTI image, store its original orientation, and reorient it.
+
+        Parameters
+        ----------
+        filename : str
+            Path to the NIfTI file.
+        to_orientation : str or tuple, optional
+            The target orientation for the in-memory data. Default is 'RAS'.
+
+        Returns
+        -------
+        StatefulImage
+            An instance of StatefulImage with data in the target orientation.
+        """
+        img = nib.load(filename)
+
+        original_affine = img.affine.copy()
+        original_axcodes = nib.orientations.aff2axcodes(img.affine)
+        original_dims = img.header.get_data_shape()
+        original_voxel_sizes = img.header.get_zooms()
+
+        if to_orientation:
+            start_ornt = nib.orientations.io_orientation(img.affine)
+            target_ornt = nib.orientations.axcodes2ornt(to_orientation)
+            transform = nib.orientations.ornt_transform(start_ornt,
+                                                        target_ornt)
+            reoriented_img = img.as_reoriented(transform)
+        else:
+            reoriented_img = img
+
+        return cls(reoriented_img.dataobj, reoriented_img.affine,
+                   reoriented_img.header, original_affine=original_affine,
+                   original_dimensions=original_dims,
+                   original_voxel_sizes=original_voxel_sizes,
+                   original_axcodes=original_axcodes)
+
+    def save(self, filename):
+        """
+        Save the image to a file, reverting to its original orientation.
+
+        Parameters
+        ----------
+        filename : str
+            Path to save the NIfTI file.
+        """
+        self.reorient_to_original()
+        nib.save(self, filename)
+
+    def reorient_to_original(self):
+        """
+        Get a Nifti1Image object reoriented to its original orientation.
+
+        Returns
+        -------
+        nib.Nifti1Image
+            A new Nifti1Image instance in the original orientation.
+        """
+        self.reorient(self._original_axcodes)
+
+
+    def reorient(self, target_axcodes):
+        """
+        Reorient the in-memory image to a target orientation.
+
+        Parameters
+        ----------
+        target_axcodes : str or tuple
+            The target orientation axis codes (e.g., 'LPS', ('R', 'A', 'S')).
+        """
+        if target_axcodes is None:
+            raise ValueError("Target axis codes cannot be None.")
+
+        current_axcodes = nib.orientations.aff2axcodes(self.affine)
+        if current_axcodes == tuple(target_axcodes):
+            return
+
+        # Check unique are only valid axis codes
+        valid_codes = {'L', 'R', 'A', 'P', 'S', 'I'}
+        for code in target_axcodes:
+            if code not in valid_codes:
+                raise ValueError(f"Invalid axis code '{code}' in target.")
+
+        # Check L/R, A/P, S/I pairs are not both present
+        pairs = [('L', 'R'), ('A', 'P'), ('S', 'I')]
+        for pair in pairs:
+            if pair[0] in target_axcodes and pair[1] in target_axcodes:
+                raise ValueError(f"Conflicting axis codes '{pair[0]}' and "
+                                 f"'{pair[1]}' in target.")
+
+        # Check no repeated axis codes (LL, RR, etc.)
+        if len(set(target_axcodes)) != 3:
+            raise ValueError("Target axis codes must be unique.")
+
+        start_ornt = nib.orientations.axcodes2ornt(current_axcodes)
+        target_ornt = nib.orientations.axcodes2ornt(target_axcodes)
+        transform = nib.orientations.ornt_transform(start_ornt, target_ornt)
+
+        reoriented_img = self.as_reoriented(transform)
+        self.__init__(reoriented_img.dataobj, reoriented_img.affine,
+                      reoriented_img.header,
+                      original_affine=self._original_affine,
+                      original_dimensions=self._original_dimensions,
+                      original_voxel_sizes=self._original_voxel_sizes,
+                      original_axcodes=self._original_axcodes)
+
+    def to_ras(self):
+        """Convenience method to reorient in-memory data to RAS."""
+        self.reorient(('R', 'A', 'S'))
+
+    def to_lps(self):
+        """Convenience method to reorient in-memory data to LPS."""
+        self.reorient(('L', 'P', 'S'))
+
+    def to_reference(self, obj):
+        """Reorient the in-memory image to match the orientation of a reference
+           object."""
+        if isinstance(obj, StatefulImage):
+            raise TypeError('Reference object must not be a StatefulImage.')
+
+        _, _, _, voxel_order = get_reference_info(obj)
+        self.reorient(voxel_order)
+
+    @property
+    def axcodes(self):
+        """Get the axis codes for the current image orientation."""
+        return nib.orientations.aff2axcodes(self.affine)
+
+    @property
+    def original_axcodes(self):
+        """Get the axis codes for the original image orientation."""
+        return self._original_axcodes
+
+    def __str__(self):
+        """Return a string representation of the image, including orientation."""
+        base_str = super().__str__()
+        current_axcodes = self.axcodes
+        reoriented = current_axcodes != self._original_axcodes
+
+        orientation_info = (
+            f"Reorientation Information:\n"
+            f"  Original axis codes:    {self._original_axcodes}\n"
+            f"  Current axis codes:     {current_axcodes}\n"
+            f"  Reoriented from original: {reoriented}"
+        )
+
+        return f"{base_str}\n{orientation_info}"