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.*",
@@ -274,6 +275,7 @@ scil_viz_volume_screenshot_mosaic  = "scilpy.cli.scil_viz_volume_screenshot_mosa
 scil_viz_volume_screenshot  = "scilpy.cli.scil_viz_volume_screenshot:main"
 scil_volume_apply_transform  = "scilpy.cli.scil_volume_apply_transform:main"
 scil_volume_b0_synthesis  = "scilpy.cli.scil_volume_b0_synthesis:main"
+scil_volume_modify_voxel_order = "scilpy.cli.scil_volume_modify_voxel_order:main"
 scil_volume_count_non_zero_voxels  = "scilpy.cli.scil_volume_count_non_zero_voxels:main"
 scil_volume_crop  = "scilpy.cli.scil_volume_crop:main"
 scil_volume_distance_map  = "scilpy.cli.scil_volume_distance_map:main"

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/cli/scil_volume_apply_transform.py

@@ -17,6 +17,7 @@
 from scilpy.io.utils import (add_overwrite_arg, assert_inputs_exist,
                              assert_outputs_exist, add_verbose_arg,
                              load_matrix_in_any_format)
+from scilpy.io.stateful_image import StatefulImage
 from scilpy.utils.filenames import split_name_with_nii
 from scilpy.version import version_string
 
@@ -72,15 +73,15 @@ def main():
     transfo = load_matrix_in_any_format(args.in_transfo)
     if args.inverse:
         transfo = np.linalg.inv(transfo)
-    moving = nib.load(args.in_file)
-    reference = nib.load(args.in_target_file)
+    moving = StatefulImage.load(args.in_file)
+    reference = StatefulImage.load(args.in_target_file)
 
     # Processing, saving
     warped_img = apply_transform(
         transfo, reference, moving, keep_dtype=args.keep_dtype,
         interp=args.interpolation)
 
-    nib.save(warped_img, args.out_name)
+    warped_img.save(args.out_name)
 
 
 if __name__ == "__main__":

src/scilpy/cli/scil_volume_crop.py

@@ -21,7 +21,6 @@
 import argparse
 import logging
 
-import nibabel as nib
 import numpy as np
 
 from scilpy.io.utils import (add_overwrite_arg,
@@ -32,6 +31,7 @@
 from scilpy.image.utils import compute_nifti_bounding_box
 from scilpy.image.volume_operations import crop_volume
 from scilpy.version import version_string
+from scilpy.io.stateful_image import StatefulImage
 
 
 def _build_arg_parser():
@@ -73,23 +73,23 @@ def main():
     assert_inputs_exist(parser, args.in_image, args.input_bbox)
     assert_outputs_exist(parser, args, args.out_image, args.output_bbox)
 
-    img = nib.load(args.in_image)
+    simg = StatefulImage.load(args.in_image)
     if args.input_bbox:
         wbbox = WorldBoundingBox.load(args.input_bbox,
                                       args.use_deprecated_pickle)
         if not args.ignore_voxel_size:
-            voxel_size = img.header.get_zooms()[0:3]
+            voxel_size = simg.header.get_zooms()[0:3]
             if not np.allclose(voxel_size, wbbox.voxel_size[0:3], atol=1e-03):
                 raise IOError("Bounding box and data voxel sizes are not "
                               "compatible. Use option --ignore_voxel_size "
                               "to ignore this test.")
     else:
-        wbbox = compute_nifti_bounding_box(img)
+        wbbox = compute_nifti_bounding_box(simg)
         if args.output_bbox:
             wbbox.dump(args.output_bbox, args.use_deprecated_pickle)
 
-    out_nifti_file = crop_volume(img, wbbox)
-    nib.save(out_nifti_file, args.out_image)
+    out_simg = crop_volume(simg, wbbox)
+    out_simg.save(args.out_image)
 
 
 if __name__ == "__main__":

src/scilpy/cli/scil_volume_flip.py

@@ -1,7 +1,17 @@
 #! /usr/bin/env python3
-
+# -*- coding: utf-8 -*-
 """
-Flip the volume according to the specified axis.
+Flip the volume according to the specified axis. In this script, axes are
+referred to as 'x', 'y' and 'z', but they simply correspond to the first,
+second and third dimensions of the data array.
+
+This script only flips the data array in memory and does not modify the
+image's strides or orientation information in the header. It simply
+flips the numpy data.
+
+
+In contrast, `scil_volume_modify_voxel_order` modifies the image header's
+voxel order, but does not modify the data array.
 """
 
 import argparse

src/scilpy/cli/scil_volume_modify_voxel_order.py

@@ -0,0 +1,83 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Change the voxel order (strides) of a NIfTI image.
+
+This script allows you to change the voxel order of a NIfTI image by modifying
+its header. The voxel order, also known as strides, defines the orientation of
+the image data in memory. This can be useful for compatibility with different
+software packages that expect a specific voxel order.
+In contrast, `scil_volume_flip` only flips the data array in memory,
+without changing the header's orientation information.
+
+The new voxel order can be specified in several ways:
+- As a string of 3 characters, e.g., 'RAS', 'LPS', 'ASR'.
+- As a comma-separated string of 3 characters, e.g., 'R,A,S'.
+- As a string of 3 numbers, e.g., '123', '231', '-12-3'.
+- As a comma-separated string of 3 numbers, e.g., '1,2,3', '-1,2,-3'.
+
+For numeric input, 1, 2, and 3 correspond to the R, A, and S axes of the
+image when loaded in RAS orientation. A negative sign flips the axis.
+For example., '-1,2,-3' would correspond to a voxel order of 'LAS'.
+
+For 4D images, the voxel order must be specified numerically.
+e.g., '1,2,3,4' or '1,2,3' (if the 4th dimension is time and does not
+need to be reordered). The 4th dimension must be 4 or -4.
+
+To change the header of a tractogram (.trk), we recommend converting it to a
+.tck file, then converting it back to .trk with the target NIfTI image as a
+reference.
+"""
+
+import argparse
+import logging
+import nibabel as nib
+
+from scilpy.io.utils import (add_overwrite_arg,
+                             add_verbose_arg,
+                             assert_inputs_exist,
+                             assert_outputs_exist)
+from scilpy.utils.orientation import parse_voxel_order
+from scilpy.io.stateful_image import StatefulImage
+from scilpy.version import version_string
+
+
+def _build_arg_parser():
+    p = argparse.ArgumentParser(description=__doc__,
+                                formatter_class=argparse.RawTextHelpFormatter,
+                                epilog=version_string)
+
+    p.add_argument('in_image',
+                   help='Path of the NIfTI file to modify.')
+    p.add_argument('out_image',
+                   help='Path of the modified NIfTI file to write.')
+    p.add_argument('--new_voxel_order', required=True,
+                   help='The new voxel order (e.g., "RAS", "1,2,3").')
+
+    add_verbose_arg(p)
+    add_overwrite_arg(p)
+
+    return p
+
+
+def main():
+    parser = _build_arg_parser()
+    args = parser.parse_args()
+    logging.getLogger().setLevel(logging.getLevelName(args.verbose))
+
+    assert_inputs_exist(parser, args.in_image)
+    assert_outputs_exist(parser, args, args.out_image)
+
+    img = nib.load(args.in_image)
+    simg = StatefulImage.load(args.in_image)
+
+    parsed_voxel_order = parse_voxel_order(args.new_voxel_order,
+                                           dimensions=len(img.shape))
+
+    simg.reorient(parsed_voxel_order)
+
+    nib.save(simg, args.out_image)
+
+
+if __name__ == "__main__":
+    main()

src/scilpy/cli/scil_volume_resample.py

@@ -24,6 +24,7 @@
                              assert_inputs_exist, assert_outputs_exist)
 from scilpy.image.volume_operations import resample_volume
 from scilpy.version import version_string
+from scilpy.io.stateful_image import StatefulImage
 
 
 def _build_arg_parser():
@@ -95,39 +96,39 @@ def main():
 
     logging.info('Loading raw data from %s', args.in_image)
 
-    img = nib.load(args.in_image)
+    simg = StatefulImage.load(args.in_image)
 
     ref_img = None
     if args.ref:
         ref_img = nib.load(args.ref)
 
         # Must not verify that headers are compatible. But can verify that, at
         # least, the first columns of their affines are compatible.
-        img_zoom_invert = [1 / zoom for zoom in img.header.get_zooms()]
+        img_zoom_invert = [1 / zoom for zoom in simg.header.get_zooms()]
         ref_zoom_invert = [1 / zoom for zoom in ref_img.header.get_zooms()]
 
-        img_affine = np.dot(img.affine[:3, :3], img_zoom_invert)
+        img_affine = np.dot(simg.affine[:3, :3], img_zoom_invert)
         ref_affine = np.dot(ref_img.affine[:3, :3], ref_zoom_invert)
 
         if not np.allclose(img_affine, ref_affine):
             parser.error("The --ref image should have the same affine as the "
                          "input image (but with a different sampling).")
 
     # Resampling volume
-    resampled_img = resample_volume(img, ref_img=ref_img,
-                                    volume_shape=args.volume_size,
-                                    iso_min=args.iso_min,
-                                    voxel_res=args.voxel_size,
-                                    interp=args.interp,
-                                    enforce_dimensions=args.enforce_dimensions)
+    resampled_simg = resample_volume(simg, ref_img=ref_img,
+                                     volume_shape=args.volume_size,
+                                     iso_min=args.iso_min,
+                                     voxel_res=args.voxel_size,
+                                     interp=args.interp,
+                                     enforce_dimensions=args.enforce_dimensions)
 
     # Saving results
-    zooms = list(resampled_img.header.get_zooms())
+    zooms = list(resampled_simg.header.get_zooms())
     if args.voxel_size:
         if len(args.voxel_size) == 1:
             args.voxel_size = args.voxel_size * 3
 
-        if not np.array_equal(zooms[:3], args.voxel_size):
+        if not np.allclose(zooms[:3], args.voxel_size, atol=1e-3):
             logging.warning('Voxel size is different from expected.'
                             ' Got: %s, expected: %s',
                             tuple(zooms), tuple(args.voxel_size))
@@ -137,10 +138,10 @@ def main():
                 zooms[0] = args.voxel_size[0]
                 zooms[1] = args.voxel_size[1]
                 zooms[2] = args.voxel_size[2]
-                resampled_img.header.set_zooms(tuple(zooms))
+                resampled_simg.header.set_zooms(tuple(zooms))
 
     logging.info('Saving resampled data to %s', args.out_image)
-    nib.save(resampled_img, args.out_image)
+    resampled_simg.save(args.out_image)
 
 
 if __name__ == '__main__':

src/scilpy/cli/scil_volume_reshape.py

@@ -27,6 +27,7 @@
                              assert_inputs_exist, assert_outputs_exist)
 from scilpy.image.volume_operations import reshape_volume
 from scilpy.version import version_string
+from scilpy.io.stateful_image import StatefulImage
 
 
 def _build_arg_parser():
@@ -80,7 +81,7 @@ def main():
 
     logging.info('Loading raw data from %s', args.in_image)
 
-    img = nib.load(args.in_image)
+    simg = StatefulImage.load(args.in_image)
 
     ref_img = None
     if args.ref:
@@ -93,14 +94,14 @@ def main():
             volume_shape = args.volume_size
 
     # Resampling volume
-    reshaped_img = reshape_volume(img, volume_shape,
-                                  mode=args.mode,
-                                  cval=args.constant_value,
-                                  dtype=args.data_type)
+    reshaped_simg = reshape_volume(simg, volume_shape,
+                                   mode=args.mode,
+                                   cval=args.constant_value,
+                                   dtype=args.data_type)
 
     # Saving results
     logging.info('Saving reshaped data to %s', args.out_image)
-    nib.save(reshaped_img, args.out_image)
+    reshaped_simg.save(args.out_image)
 
 
 if __name__ == '__main__':