From 3f14ce4b73d41f13f60b1c953d89904bcd9a414a Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Thu, 14 Mar 2024 12:07:19 -0400
Subject: [PATCH 01/14] volume_apply_transform: avoid loading fdata twice. Move
 warning

---
 scilpy/image/volume_operations.py      |  4 ++++
 scripts/scil_volume_apply_transform.py | 17 ++++++-----------
 2 files changed, 10 insertions(+), 11 deletions(-)

diff --git a/scilpy/image/volume_operations.py b/scilpy/image/volume_operations.py
index e8ede37ad..0ade41040 100644
--- a/scilpy/image/volume_operations.py
+++ b/scilpy/image/volume_operations.py
@@ -152,6 +152,10 @@ def apply_transform(transfo, reference, moving,
     elif moving_data.ndim == 4:
         if isinstance(moving_data[0, 0, 0], np.void):
             raise ValueError('Does not support TrackVis RGB')
+        else:
+            logging.warning('You are applying a transform to a 4D volume. If'
+                            'it is a DWI volume, make sure to rotate your '
+                            'bvecs with scil_gradients_apply_transform.py')
 
         affine_map = AffineMap(np.linalg.inv(transfo),
                                dim[0:3], grid2world,
diff --git a/scripts/scil_volume_apply_transform.py b/scripts/scil_volume_apply_transform.py
index 6a6085dde..7ec3c058b 100755
--- a/scripts/scil_volume_apply_transform.py
+++ b/scripts/scil_volume_apply_transform.py
@@ -53,14 +53,11 @@ def main():
     args = parser.parse_args()
     logging.getLogger().setLevel(logging.getLevelName(args.verbose))
 
+    # Verifications
     assert_inputs_exist(parser, [args.in_file, args.in_target_file,
                                  args.in_transfo])
     assert_outputs_exist(parser, args, args.out_name)
 
-    transfo = load_matrix_in_any_format(args.in_transfo)
-    if args.inverse:
-        transfo = np.linalg.inv(transfo)
-
     _, ref_extension = split_name_with_nii(args.in_target_file)
     _, in_extension = split_name_with_nii(args.in_file)
     if ref_extension not in ['.nii', '.nii.gz']:
@@ -69,16 +66,14 @@ def main():
     if in_extension not in ['.nii', '.nii.gz']:
         parser.error('{} is an unsupported format.'.format(args.in_file))
 
-    # Load images and validate input type.
+    # Loading
+    transfo = load_matrix_in_any_format(args.in_transfo)
+    if args.inverse:
+        transfo = np.linalg.inv(transfo)
     moving = nib.load(args.in_file)
-
-    if moving.get_fdata().ndim == 4:
-        logging.warning('You are applying a transform to a 4D dwi volume, '
-                        'make sure to rotate your bvecs with '
-                        'scil_gradients_apply_transform.py')
-
     reference = nib.load(args.in_target_file)
 
+    # Processing, saving
     warped_img = apply_transform(
         transfo, reference, moving, keep_dtype=args.keep_dtype)
 

From cd8f772edfca96213c832729778a102803dd43b6 Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Thu, 14 Mar 2024 12:24:25 -0400
Subject: [PATCH 02/14] Encapsulate volume_b0_synthesis

---
 scilpy/image/volume_b0_synthesis.py | 91 +++++++++++++++++++++++++++
 scripts/scil_volume_b0_synthesis.py | 98 ++++++++---------------------
 2 files changed, 117 insertions(+), 72 deletions(-)
 create mode 100644 scilpy/image/volume_b0_synthesis.py

diff --git a/scilpy/image/volume_b0_synthesis.py b/scilpy/image/volume_b0_synthesis.py
new file mode 100644
index 000000000..01b5dd4da
--- /dev/null
+++ b/scilpy/image/volume_b0_synthesis.py
@@ -0,0 +1,91 @@
+# -*- coding: utf-8 -*-
+import logging
+import os
+import warnings
+
+# Disable tensorflow warnings
+with warnings.catch_warnings():
+    os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3'
+    warnings.simplefilter("ignore")
+    from dipy.nn.synb0 import Synb0
+
+import numpy as np
+from dipy.align.imaffine import AffineMap
+from dipy.segment.tissue import TissueClassifierHMRF
+from scipy.ndimage import gaussian_filter
+
+from scilpy.image.volume_operations import register_image
+
+
+def compute_b0_synthesis(t1_data, t1_bet_data, b0_data, b0_bet_data,
+                         template_data, t1_affine, b0_affine, template_affine,
+                         verbose):
+    """
+    Parameters
+    ----------
+    t1_data: np.ndarary
+        The T1 (wholebrain, with the skull).
+    t1_bet_data: np.ndarray
+        The mask.
+    b0_data: np.ndarray
+        The b0 (wholebrain, with the skull).
+    b0_bet_data: np.ndarray
+        The mask.
+    template_data: np.ndarray
+        The template for typical usage.
+    t1_affine: np.ndarray
+        The T1's affine.
+    b0_affine: np.ndarray
+        The b0's affine.
+    template_affine: np.ndarray
+        The template's affine
+    verbose: bool
+        Whether to make dipy's Synb0 verbose.
+
+    Returns
+    -------
+    rev_b0: np.ndarray
+        The synthetized b0.
+    """
+    logging.info('The usage of synthetic b0 is not fully tested.'
+                 'Be careful when using it.')
+
+    # Crude estimation of the WM mean intensity and normalization
+    logging.info('Estimating WM mean intensity')
+    hmrf = TissueClassifierHMRF()
+    t1_bet_data = gaussian_filter(t1_bet_data, 2)
+    _, final_segmentation, _ = hmrf.classify(t1_bet_data, 3, 0.25,
+                                             tolerance=1e-4, max_iter=5)
+    avg_wm = np.mean(t1_data[final_segmentation == 3])
+
+    # Modifying t1
+    t1_data /= avg_wm
+    t1_data *= 110
+
+    # SyNB0 works only in a standard space, so we need to register the images
+    logging.info('Registering images')
+    # Use the BET image for registration
+    t1_bet_to_b0, t1_bet_to_b0_transform = register_image(
+        b0_bet_data, b0_affine, t1_bet_data, t1_affine, fine=True)
+    affine_map = AffineMap(t1_bet_to_b0_transform,
+                           b0_data.shape, b0_affine,
+                           t1_data.shape, t1_affine)
+    t1_skull_to_b0 = affine_map.transform(t1_data.astype(np.float64))
+
+    # Then register to MNI (using the BET again)
+    _, t1_bet_to_b0_to_mni_transform = register_image(
+        template_data, template_affine, t1_bet_to_b0, b0_affine, fine=True)
+    affine_map = AffineMap(t1_bet_to_b0_to_mni_transform,
+                           template_data.shape, template_affine,
+                           b0_data.shape, b0_affine)
+
+    # But for prediction, we want the skull
+    b0_skull_to_mni = affine_map.transform(b0_data.astype(np.float64))
+    t1_skull_to_mni = affine_map.transform(t1_skull_to_b0.astype(np.float64))
+
+    logging.info('Running SyN-B0')
+    SyNb0 = Synb0(verbose)
+    rev_b0 = SyNb0.predict(b0_skull_to_mni, t1_skull_to_mni)
+    rev_b0 = affine_map.transform_inverse(rev_b0.astype(np.float64))
+
+    return rev_b0
diff --git a/scripts/scil_volume_b0_synthesis.py b/scripts/scil_volume_b0_synthesis.py
index 79516b8a5..8ffb81fdb 100644
--- a/scripts/scil_volume_b0_synthesis.py
+++ b/scripts/scil_volume_b0_synthesis.py
@@ -8,7 +8,7 @@
 template, run SyNb0 and then transform the result back to the original space.
 
 SyNb0 is a deep learning model that predicts a synthetic a distortion-free
-b0 image from a distorted b0 and T1w
+b0 image from a distorted b0 and T1w.
 
 This script must be used carefully, as it is not meant to be used in an
 environment with the following dependencies already installed (not default
@@ -20,33 +20,22 @@
 
 import argparse
 import logging
-import os
-import warnings
 
-# Disable tensorflow warnings
-with warnings.catch_warnings():
-    os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3'
-    warnings.simplefilter("ignore")
-    from dipy.nn.synb0 import Synb0
-
-from dipy.align.imaffine import AffineMap
-from dipy.segment.tissue import TissueClassifierHMRF
 import nibabel as nib
-import numpy as np
-from scipy.ndimage import gaussian_filter
 
+from scilpy.image.volume_b0_synthesis import compute_b0_synthesis
+from scilpy.io.image import get_data_as_mask
 from scilpy.io.fetcher import get_synb0_template_path
-from scilpy.io.utils import (add_overwrite_arg,
-                             add_verbose_arg,
-                             assert_inputs_exist,
-                             assert_outputs_exist)
-from scilpy.image.volume_operations import register_image
+from scilpy.io.utils import (add_overwrite_arg, add_verbose_arg,
+                             assert_inputs_exist, assert_outputs_exist,
+                             assert_headers_compatible)
 
 EPILOG = """
 [1] Schilling, Kurt G., et al. "Synthesized b0 for diffusion distortion
   correction (Synb0-DisCo)." Magnetic resonance imaging 64 (2019): 62-70.
 """
 
+
 def _build_arg_parser():
     p = argparse.ArgumentParser(
         description=__doc__,
@@ -72,73 +61,38 @@ def _build_arg_parser():
 def main():
     parser = _build_arg_parser()
     args = parser.parse_args()
-    assert_inputs_exist(parser, [args.in_b0, args.in_t1])
+    assert_inputs_exist(parser, [args.in_b0, args.in_b0_mask,
+                                 args.in_t1, args.in_t1_mask])
     assert_outputs_exist(parser, args, args.out_b0)
+    assert_headers_compatible(parser, [args.in_b0, args.in_b0_mask])
+    assert_headers_compatible(parser, [args.in_t1, args.in_t1_mask])
 
     logging.getLogger().setLevel(logging.getLevelName(args.verbose))
     logging.info('The usage of synthetic b0 is not fully tested.'
                  'Be careful when using it.')
 
+    # Loading
     template_img = nib.load(get_synb0_template_path())
     template_data = template_img.get_fdata()
 
     b0_img = nib.load(args.in_b0)
-    b0_skull_data = b0_img.get_fdata()
-    b0_mask_img = nib.load(args.in_b0_mask)
-    b0_mask_data = b0_mask_img.get_fdata()
+    b0_data = b0_img.get_fdata()
+    b0_mask_data = get_data_as_mask(nib.load(args.in_b0_mask))
 
     t1_img = nib.load(args.in_t1)
-    t1_skull_data = t1_img.get_fdata()
-    t1_mask_img = nib.load(args.in_t1_mask)
-    t1_mask_data = t1_mask_img.get_fdata()
-
-    b0_bet_data = np.zeros(b0_skull_data.shape)
-    b0_bet_data[b0_mask_data > 0] = b0_skull_data[b0_mask_data > 0]
-    t1_bet_data = np.zeros(t1_skull_data.shape)
-    t1_bet_data[t1_mask_data > 0] = t1_skull_data[t1_mask_data > 0]
-
-    # Crude estimation of the WM mean intensity and normalization
-    logging.info('Estimating WM mean intensity')
-    hmrf = TissueClassifierHMRF()
-    t1_bet_data = gaussian_filter(t1_bet_data, 2)
-    _, final_segmentation, _ = hmrf.classify(t1_bet_data, 3, 0.25,
-                                             tolerance=1e-4, max_iter=5)
-    avg_wm = np.mean(t1_skull_data[final_segmentation == 3])
-    t1_skull_data /= avg_wm
-    t1_skull_data *= 110
-
-    # SyNB0 works only in a standard space, so we need to register the images
-    logging.info('Registering images')
-    # Use the BET image for registration
-    t1_bet_to_b0, t1_bet_to_b0_transform = register_image(b0_bet_data,
-                                                          b0_img.affine,
-                                                          t1_bet_data,
-                                                          t1_img.affine,
-                                                          fine=True)
-    affine_map = AffineMap(t1_bet_to_b0_transform,
-                           b0_skull_data.shape, b0_img.affine,
-                           t1_skull_data.shape, t1_img.affine)
-    t1_skull_to_b0 = affine_map.transform(t1_skull_data.astype(np.float64))
-
-    # Then register to MNI (using the BET again)
-    _, t1_bet_to_b0_to_mni_transform = register_image(template_data,
-                                                      template_img.affine,
-                                                      t1_bet_to_b0,
-                                                      b0_img.affine,
-                                                      fine=True)
-    affine_map = AffineMap(t1_bet_to_b0_to_mni_transform,
-                           template_data.shape, template_img.affine,
-                           b0_skull_data.shape, b0_img.affine)
-
-    # But for prediction, we want the skull
-    b0_skull_to_mni = affine_map.transform(b0_skull_data.astype(np.float64))
-    t1_skull_to_mni = affine_map.transform(t1_skull_to_b0.astype(np.float64))
-
-    logging.info('Running SyN-B0')
-    SyNb0 = Synb0(args.verbose)
-    rev_b0 = SyNb0.predict(b0_skull_to_mni, t1_skull_to_mni)
-    rev_b0 = affine_map.transform_inverse(rev_b0.astype(np.float64))
+    t1_data = t1_img.get_fdata()
+    t1_mask_data = get_data_as_mask(nib.load(args.in_t1_mask))
+
+    b0_bet_data = b0_data * b0_mask_data
+    t1_bet_data = t1_data * t1_mask_data
+
+    # Processing
+    verbose = True if args.verbose in ['INFO', 'DEBUG'] else False
+    rev_b0 = compute_b0_synthesis(t1_data, t1_bet_data, b0_data, b0_bet_data,
+                                  template_data, t1_img.affine, b0_img.affine,
+                                  template_img.affine, verbose)
 
+    # Saving
     dtype = b0_img.get_data_dtype()
     nib.save(nib.Nifti1Image(rev_b0.astype(dtype), b0_img.affine),
              args.out_b0)

From 089b52202783bbe8d1af572d870ccd00977f08ae Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Thu, 14 Mar 2024 13:35:11 -0400
Subject: [PATCH 03/14] pytest skipping --help test. Fixing. Local run: fix
 t1_data (>0 but negative data)

---
 scilpy/image/volume_b0_synthesis.py       |  5 +++
 scripts/scil_volume_b0_synthesis.py       |  6 +--
 scripts/tests/test_volume_b0_synthesis.py | 52 ++++++++++++-----------
 3 files changed, 35 insertions(+), 28 deletions(-)

diff --git a/scilpy/image/volume_b0_synthesis.py b/scilpy/image/volume_b0_synthesis.py
index 01b5dd4da..303ce8ef7 100644
--- a/scilpy/image/volume_b0_synthesis.py
+++ b/scilpy/image/volume_b0_synthesis.py
@@ -21,6 +21,11 @@ def compute_b0_synthesis(t1_data, t1_bet_data, b0_data, b0_bet_data,
                          template_data, t1_affine, b0_affine, template_affine,
                          verbose):
     """
+    Note. Tensorflow is required here, through dipy.Synb0. If not installed,
+    dipy will raise an error, like:
+    >> dipy.utils.tripwire.TripWireError: We need package tensorflow_addons for
+    these functions, but ``import tensorflow_addons`` raised an ImportError.
+
     Parameters
     ----------
     t1_data: np.ndarary
diff --git a/scripts/scil_volume_b0_synthesis.py b/scripts/scil_volume_b0_synthesis.py
index 8ffb81fdb..c12ba1a91 100644
--- a/scripts/scil_volume_b0_synthesis.py
+++ b/scripts/scil_volume_b0_synthesis.py
@@ -10,9 +10,9 @@
 SyNb0 is a deep learning model that predicts a synthetic a distortion-free
 b0 image from a distorted b0 and T1w.
 
-This script must be used carefully, as it is not meant to be used in an
-environment with the following dependencies already installed (not default
-in Scilpy):
+This script must be used carefully, as it is meant to be used in an
+environment with the following dependencies already installed (not installed by
+default in Scilpy):
 - tensorflow-addons
 - tensorrt
 - tensorflow
diff --git a/scripts/tests/test_volume_b0_synthesis.py b/scripts/tests/test_volume_b0_synthesis.py
index 003083740..a9c811dc7 100644
--- a/scripts/tests/test_volume_b0_synthesis.py
+++ b/scripts/tests/test_volume_b0_synthesis.py
@@ -7,12 +7,16 @@
 import nibabel as nib
 import numpy as np
 import pytest
+from _pytest.outcomes import Skipped
 
 from scilpy import SCILPY_HOME
 from scilpy.io.fetcher import fetch_data, get_testing_files_dict
 
-tensorflow = pytest.importorskip("tensorflow")
-
+try:
+    tensorflow = pytest.importorskip("tensorflow")
+    have_tensorflow = True
+except Skipped:
+    have_tensorflow = False
 
 # If they already exist, this only takes 5 seconds (check md5sum)
 fetch_data(get_testing_files_dict(), keys=['others.zip', 'processing.zip'])
@@ -24,27 +28,25 @@ def test_help_option(script_runner):
     assert ret.success
 
 
-@pytest.mark.skipif(tensorflow is None, reason="Tensorflow not installed")
 def test_synthesis(script_runner, monkeypatch):
-    monkeypatch.chdir(os.path.expanduser(tmp_dir.name))
-    in_t1 = os.path.join(SCILPY_HOME, 'others',
-                         't1.nii.gz')
-    in_b0 = os.path.join(SCILPY_HOME, 'processing',
-                         'b0_mean.nii.gz')
-
-    t1_img = nib.load(in_t1)
-    b0_img = nib.load(in_b0)
-    t1_data = t1_img.get_fdata()
-    b0_data = b0_img.get_fdata()
-    t1_data[t1_data > 0] = 1
-    b0_data[b0_data > 0] = 1
-    nib.save(nib.Nifti1Image(t1_data.astype(np.uint8), t1_img.affine),
-             't1_mask.nii.gz')
-    nib.save(nib.Nifti1Image(b0_data.astype(np.uint8), b0_img.affine),
-             'b0_mask.nii.gz')
-
-    ret = script_runner.run('scil_volume_b0_synthesis.py',
-                            in_t1, 't1_mask.nii.gz',
-                            in_b0, 'b0_mask.nii.gz',
-                            'b0_synthesized.nii.gz', '-v')
-    assert ret.success
+    if have_tensorflow:
+        monkeypatch.chdir(os.path.expanduser(tmp_dir.name))
+        in_t1 = os.path.join(SCILPY_HOME, 'others', 't1.nii.gz')
+        in_b0 = os.path.join(SCILPY_HOME, 'processing', 'b0_mean.nii.gz')
+
+        t1_img = nib.load(in_t1)
+        b0_img = nib.load(in_b0)
+        t1_data = t1_img.get_fdata()
+        b0_data = b0_img.get_fdata()
+        t1_data[t1_data != 0] = 1
+        b0_data[b0_data > 0] = 1
+        nib.save(nib.Nifti1Image(t1_data.astype(np.uint8), t1_img.affine),
+                 't1_mask.nii.gz')
+        nib.save(nib.Nifti1Image(b0_data.astype(np.uint8), b0_img.affine),
+                 'b0_mask.nii.gz')
+
+        ret = script_runner.run('scil_volume_b0_synthesis.py',
+                                in_t1, 't1_mask.nii.gz',
+                                in_b0, 'b0_mask.nii.gz',
+                                'b0_synthesized.nii.gz', '-v')
+        assert ret.success

From 0cc93de04d7294402349f4df6b8c148f46ed53a3 Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Thu, 14 Mar 2024 14:33:50 -0400
Subject: [PATCH 04/14] Add more coverage for count_non_zero_voxel

---
 scripts/scil_volume_count_non_zero_voxels.py  | 22 +++++++++++++------
 .../test_volume_count_non_zero_voxels.py      | 18 ++++++++++++---
 2 files changed, 30 insertions(+), 10 deletions(-)

diff --git a/scripts/scil_volume_count_non_zero_voxels.py b/scripts/scil_volume_count_non_zero_voxels.py
index c8b984df9..e58cb7184 100755
--- a/scripts/scil_volume_count_non_zero_voxels.py
+++ b/scripts/scil_volume_count_non_zero_voxels.py
@@ -20,7 +20,8 @@
 import numpy as np
 
 from scilpy.image.volume_operations import count_non_zero_voxels
-from scilpy.io.utils import assert_inputs_exist, add_verbose_arg
+from scilpy.io.utils import assert_inputs_exist, add_verbose_arg, \
+    assert_outputs_exist, add_overwrite_arg
 
 
 def _build_arg_parser():
@@ -31,13 +32,14 @@ def _build_arg_parser():
 
     p.add_argument(
         '--out', metavar='OUT_FILE', dest='out_filename',
-        help='name of the output file, which will be saved as a text file.')
+        help='Name of the output file, which will be saved as a text file.')
     p.add_argument(
         '--stats', action='store_true', dest='stats_format',
-        help='output the value using a stats format. Using this syntax will '
-             'append\na line to the output file, instead of creating a file '
-             'with only one line.\nThis is useful to create a file to be used '
-             'as the source of data for a graph.\nCan be combined with --id')
+        help='If set, output the value using a stats format. Using this '
+             'synthax will append\na line to the output file, instead of '
+             'creating a file with only one line.\nThis is useful to create '
+             'a file to be used as the source of data for a graph.\nCan be '
+             'combined with --id')
     p.add_argument(
         '--id', dest='value_id',
         help='Id of the current count. If used, the value of this argument '
@@ -45,6 +47,7 @@ def _build_arg_parser():
              'Mostly useful with --stats.')
 
     add_verbose_arg(p)
+    add_overwrite_arg(p)
 
     return p
 
@@ -55,7 +58,12 @@ def main():
     logging.getLogger().setLevel(logging.getLevelName(args.verbose))
 
     assert_inputs_exist(parser, args.in_image)
-    # out_filename can exist or not
+    if not args.stats_format:
+        assert_outputs_exist(parser, args, [], args.out_filename)
+    elif args.overwrite:
+        parser.error("Using -f together with --stats is unclear. Please "
+                     "either remove --stats to start anew, or remove -f to "
+                     "append to existing file.")
 
     # Load image file
     im = nib.load(args.in_image).get_fdata(dtype=np.float32)
diff --git a/scripts/tests/test_volume_count_non_zero_voxels.py b/scripts/tests/test_volume_count_non_zero_voxels.py
index 6b000e08d..cec9c23bc 100644
--- a/scripts/tests/test_volume_count_non_zero_voxels.py
+++ b/scripts/tests/test_volume_count_non_zero_voxels.py
@@ -17,9 +17,21 @@ def test_help_option(script_runner):
     assert ret.success
 
 
-def test_execution_others(script_runner, monkeypatch):
+def test_execution_simple_print(script_runner, monkeypatch):
     monkeypatch.chdir(os.path.expanduser(tmp_dir.name))
-    in_img = os.path.join(SCILPY_HOME, 'others',
-                          'rgb.nii.gz')
+    in_img = os.path.join(SCILPY_HOME, 'others', 'rgb.nii.gz')
     ret = script_runner.run('scil_volume_count_non_zero_voxels.py', in_img)
     assert ret.success
+
+
+def test_execution_save_in_file(script_runner):
+    os.chdir(os.path.expanduser(tmp_dir.name))
+    in_img = os.path.join(SCILPY_HOME, 'others', 'rgb.nii.gz')
+    ret = script_runner.run('scil_volume_count_non_zero_voxels.py', in_img,
+                            '--out', 'printed.txt')
+    assert ret.success
+
+    # Then re-use the same out file with --stats
+    ret = script_runner.run('scil_volume_count_non_zero_voxels.py', in_img,
+                            '--out', 'printed.txt', '--stats')
+    assert ret.success

From 9d4263352709ada545724ed7125685218e16ae09 Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Thu, 14 Mar 2024 14:39:19 -0400
Subject: [PATCH 05/14] volume_crop: add more coverage to test

---
 scripts/scil_volume_crop.py       | 2 +-
 scripts/tests/test_volume_crop.py | 8 +++++++-
 2 files changed, 8 insertions(+), 2 deletions(-)

diff --git a/scripts/scil_volume_crop.py b/scripts/scil_volume_crop.py
index 71230fc08..b4982d813 100755
--- a/scripts/scil_volume_crop.py
+++ b/scripts/scil_volume_crop.py
@@ -50,7 +50,7 @@ def _build_arg_parser():
                          'the bounding box to crop input file.')
     g1.add_argument('--output_bbox',
                     help='Path of the pickle file where to write the '
-                         'computed bounding box.')
+                         'computed bounding box. (.pickle extension)')
     return p
 
 
diff --git a/scripts/tests/test_volume_crop.py b/scripts/tests/test_volume_crop.py
index f4ad89a8d..13890880b 100644
--- a/scripts/tests/test_volume_crop.py
+++ b/scripts/tests/test_volume_crop.py
@@ -20,5 +20,11 @@ def test_help_option(script_runner):
 def test_execution_processing(script_runner, monkeypatch):
     monkeypatch.chdir(os.path.expanduser(tmp_dir.name))
     in_dwi = os.path.join(SCILPY_HOME, 'processing', 'dwi.nii.gz')
-    ret = script_runner.run('scil_volume_crop.py', in_dwi, 'dwi_crop.nii.gz')
+    ret = script_runner.run('scil_volume_crop.py', in_dwi, 'dwi_crop.nii.gz',
+                            '--output_bbox', 'bbox.pickle')
+    assert ret.success
+
+    # Then try to load back the same box
+    ret = script_runner.run('scil_crop_volume.py', in_dwi, 'dwi_crop.nii.gz',
+                            '--input_bbox', 'bbox.pickle')
     assert ret.success

From f992649f7f1b7f32966575319be4bd8e17bf1147 Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Thu, 14 Mar 2024 14:49:19 -0400
Subject: [PATCH 06/14] volume_outliers_ransac: cleaning

---
 scilpy/image/volume_operations.py             | 42 +++++++++++++++++++
 scripts/scil_volume_remove_outliers_ransac.py | 39 ++++++-----------
 2 files changed, 55 insertions(+), 26 deletions(-)

diff --git a/scilpy/image/volume_operations.py b/scilpy/image/volume_operations.py
index 0ade41040..8ec3a2826 100644
--- a/scilpy/image/volume_operations.py
+++ b/scilpy/image/volume_operations.py
@@ -16,6 +16,7 @@
 import numpy as np
 from numpy import ma
 from scipy.ndimage import binary_dilation, gaussian_filter
+from sklearn import linear_model
 
 from scilpy.image.reslice import reslice  # Don't use Dipy's reslice. Buggy.
 from scilpy.io.image import get_data_as_mask
@@ -372,6 +373,47 @@ def compute_snr(dwi, bval, bvec, b0_thr, mask,
     return val
 
 
+def remove_outliers_ransac(in_data, min_fit, fit_thr, max_iter):
+    """
+    Remove outliers from image using the RANSAC algorithm.
+
+    Parameters
+    ----------
+    in_data: np.ndarray
+        The input.
+    min_fit: int
+        The minimum number of data values required to fit the model.
+    fit_thr: float
+        Threshold value for determining when a data point fits a model.
+    max_iter: int
+        The maximum number of iterations allowed in the algorithm.
+
+    Returns
+    -------
+    out_data: np.ndarray
+        Data without outliers.
+    """
+    init_shape = in_data.shape
+    in_data_flat = in_data.flatten()
+    in_nzr_ind = np.nonzero(in_data_flat)
+    in_nzr_val = np.array(in_data_flat[in_nzr_ind])
+
+    X = in_nzr_ind[0][:, np.newaxis]
+    model_ransac = linear_model.RANSACRegressor(
+        base_estimator=linear_model.LinearRegression(), min_samples=min_fit,
+        residual_threshold=fit_thr, max_trials=max_iter)
+    model_ransac.fit(X, in_nzr_val)
+
+    outlier_mask = np.logical_not(model_ransac.inlier_mask_)
+    outliers = X[outlier_mask]
+    logging.info('# outliers: {}'.format(len(outliers)))
+
+    in_data_flat[outliers] = 0
+    out_data = np.reshape(in_data_flat, init_shape)
+
+    return out_data
+
+
 def smooth_to_fwhm(data, fwhm):
     """
     Smooth a volume to given FWHM.
diff --git a/scripts/scil_volume_remove_outliers_ransac.py b/scripts/scil_volume_remove_outliers_ransac.py
index 7c3377b35..14da34c21 100755
--- a/scripts/scil_volume_remove_outliers_ransac.py
+++ b/scripts/scil_volume_remove_outliers_ransac.py
@@ -15,8 +15,8 @@
 
 import nibabel as nib
 import numpy as np
-from sklearn import linear_model
 
+from scilpy.image.volume_operations import remove_outliers_ransac
 from scilpy.io.utils import (add_overwrite_arg,
                              add_verbose_arg,
                              assert_inputs_exist,
@@ -33,14 +33,14 @@ def _build_arg_parser():
                    help='Corrected Nifti image.')
 
     p.add_argument('--min_fit', type=int, default=50,
-                   help='The minimum number of data values required ' +
-                        'to fit the model. [%(default)s]')
+                   help='The minimum number of data values required to fit '
+                        'the model. [%(default)s]')
     p.add_argument('--max_iter', type=int, default=1000,
-                   help='The maximum number of iterations allowed ' +
-                        'in the algorithm. [%(default)s]')
+                   help='The maximum number of iterations allowed in the '
+                        'algorithm. [%(default)s]')
     p.add_argument('--fit_thr', type=float, default=1e-2,
-                   help='Threshold value for determining when a data ' +
-                        'point fits a model. [%(default)s]')
+                   help='Threshold value for determining when a data point '
+                        'fits a model. [%(default)s]')
 
     add_verbose_arg(p)
     add_overwrite_arg(p)
@@ -53,6 +53,7 @@ def main():
     args = parser.parse_args()
     logging.getLogger().setLevel(logging.getLevelName(args.verbose))
 
+    # Verifications
     assert_inputs_exist(parser, args.in_image)
     assert_outputs_exist(parser, args, args.out_image)
 
@@ -66,6 +67,7 @@ def main():
         parser.error('--fit_thr should be greater than 0. Current value: {}'
                      .format(args.fit_thr))
 
+    # Loading
     in_img = nib.load(args.in_image)
     in_data = in_img.get_fdata(dtype=np.float32)
 
@@ -73,26 +75,11 @@ def main():
         logging.warning('Be careful, your image doesn\'t seem to be an ad, '
                         'md or rd.')
 
-    in_data_flat = in_data.flatten()
-    in_nzr_ind = np.nonzero(in_data_flat)
-    in_nzr_val = np.array(in_data_flat[in_nzr_ind])
+    # Processing
+    out_data = remove_outliers_ransac(in_data, args.min_fit, args.fit_thr,
+                                      args.max_iter)
 
-    X = in_nzr_ind[0][:, np.newaxis]
-    model_ransac = linear_model.RANSACRegressor(
-        base_estimator=linear_model.LinearRegression(),
-        min_samples=args.min_fit,
-        residual_threshold=args.fit_thr,
-        max_trials=args.max_iter)
-    model_ransac.fit(X, in_nzr_val)
-
-    outlier_mask = np.logical_not(model_ransac.inlier_mask_)
-    outliers = X[outlier_mask]
-
-    logging.info('# outliers: {}'.format(len(outliers)))
-
-    in_data_flat[outliers] = 0
-
-    out_data = np.reshape(in_data_flat, in_img.shape)
+    # Saving
     nib.save(nib.Nifti1Image(out_data, in_img.affine, in_img.header),
              args.out_image)
 

From 5dd516482cfd54a8ed59f5fb8a1b2d7eab774beb Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Thu, 14 Mar 2024 14:51:57 -0400
Subject: [PATCH 07/14] volume_outliers_ransac: cleaning

---
 scilpy/image/tests/test_volume_operations.py | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/scilpy/image/tests/test_volume_operations.py b/scilpy/image/tests/test_volume_operations.py
index 02d52e3b1..2c7e7fb47 100644
--- a/scilpy/image/tests/test_volume_operations.py
+++ b/scilpy/image/tests/test_volume_operations.py
@@ -119,6 +119,16 @@ def test_compute_snr():
     assert np.allclose(snr[0]['snr'], target_val, atol=0.00005)
 
 
+def test_remove_outliers_ransac():
+    # Could test, but uses mainly sklearn. Not testing again.
+    pass
+
+
+def smooth_to_fwhm():
+    # toDo
+    pass
+
+
 def test_resample_volume():
     moving3d = np.pad(np.ones((4, 4, 4)), pad_width=1,
                       mode='constant', constant_values=0)

From 2892bd9ba779c2f047f97d3383a74769d9c97dd8 Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Thu, 14 Mar 2024 14:53:43 -0400
Subject: [PATCH 08/14] volume_resample: ok

---
 scripts/scil_volume_resample.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/scripts/scil_volume_resample.py b/scripts/scil_volume_resample.py
index 42ed7d88a..fbe241053 100755
--- a/scripts/scil_volume_resample.py
+++ b/scripts/scil_volume_resample.py
@@ -66,6 +66,9 @@ def main():
     # Checking args
     assert_inputs_exist(parser, args.in_image, args.ref)
     assert_outputs_exist(parser, args, args.out_image)
+    # toDo after PR937
+    #assert_headers_compatible(parser, args.in_image, args.ref)
+
     if args.enforce_dimensions and not args.ref:
         parser.error("Cannot enforce dimensions without a reference image")
 

From 1d47c97ea8fdd1f148a7084bce48c7378451f6ad Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Thu, 14 Mar 2024 15:08:07 -0400
Subject: [PATCH 09/14] volume_stats_in_labels: reformat

---
 scilpy/image/labels.py                 | 47 ++++++++++++++++++-
 scilpy/image/tests/test_labels.py      |  5 ++
 scripts/scil_volume_stats_in_labels.py | 65 +++++++++-----------------
 3 files changed, 73 insertions(+), 44 deletions(-)

diff --git a/scilpy/image/labels.py b/scilpy/image/labels.py
index 90d73c3c1..4a8702102 100644
--- a/scilpy/image/labels.py
+++ b/scilpy/image/labels.py
@@ -24,7 +24,7 @@ def get_data_as_labels(in_img):
     curr_type = in_img.get_data_dtype()
 
     if np.issubdtype(curr_type, np.signedinteger) or \
-       np.issubdtype(curr_type, np.unsignedinteger):
+            np.issubdtype(curr_type, np.unsignedinteger):
         return np.asanyarray(in_img.dataobj).astype(np.uint16)
     else:
         basename = os.path.basename(in_img.get_filename())
@@ -300,3 +300,48 @@ def dilate_labels(data, vox_size, distance, nbr_processes,
     data = data.reshape(img_shape)
 
     return data
+
+
+def get_stats_in_label(map_data, label_data, label_lut):
+    """
+    Get statistics about a map for each label in an atlas.
+
+    Parameters
+    ----------
+    map_data: np.ndarray
+        The map from which to get statistics.
+    label_data: np.ndarray
+        The loaded atlas.
+    label_lut: dict
+        The loaded label LUT (look-up table).
+
+    Returns
+    -------
+    out_dict: dict
+        A dict with one key per label name, and its values are the computed
+        statistics.
+    """
+    (label_indices, label_names) = zip(*label_lut.items())
+
+    out_dict = {}
+    for label, name in zip(label_indices, label_names):
+        label = int(label)
+        if label != 0:
+            curr_data = (map_data[np.where(label_data == label)])
+            nb_vx_roi = np.count_nonzero(label_data == label)
+            nb_seed_vx = np.count_nonzero(curr_data)
+
+            if nb_seed_vx != 0:
+                mean_seed = np.sum(curr_data) / nb_seed_vx
+                max_seed = np.max(curr_data)
+                std_seed = np.sqrt(np.mean(abs(curr_data[curr_data != 0] -
+                                               mean_seed) ** 2))
+
+                out_dict[name] = {'ROI-idx': label,
+                                  'ROI-name': str(name),
+                                  'nb-vx-roi': int(nb_vx_roi),
+                                  'nb-vx-seed': int(nb_seed_vx),
+                                  'max': int(max_seed),
+                                  'mean': float(mean_seed),
+                                  'std': float(std_seed)}
+    return out_dict
diff --git a/scilpy/image/tests/test_labels.py b/scilpy/image/tests/test_labels.py
index bb7d1129a..e30685580 100644
--- a/scilpy/image/tests/test_labels.py
+++ b/scilpy/image/tests/test_labels.py
@@ -157,3 +157,8 @@ def test_split_labels():
     assert len(out_labels) == 3
     assert_equal(np.unique(out_labels[0]), [0, 6])
     assert_equal(np.unique(out_labels[1]), [0])
+
+
+def test_stats_in_labels():
+    # toDO. Will need to create a fake LUT.
+    pass
diff --git a/scripts/scil_volume_stats_in_labels.py b/scripts/scil_volume_stats_in_labels.py
index c4b926afb..ebefde0fa 100755
--- a/scripts/scil_volume_stats_in_labels.py
+++ b/scripts/scil_volume_stats_in_labels.py
@@ -1,10 +1,11 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 """
-Computes the information from the seeding map for each cortical region
-(corresponding to an atlas) associated with a specific bundle.
-Here we want to estimate the seeding attribution to cortical area
-affected by the bundle
+Computes the information from the input map for each cortical region
+(corresponding to an atlas).
+
+Hint: For instance, this script could be useful if you have a seed map from a
+specific bundle, to know from which regions it originated.
 
 Formerly: scil_compute_seed_by_labels.py
 """
@@ -16,23 +17,22 @@
 import nibabel as nib
 import numpy as np
 
-from scilpy.image.labels import get_data_as_labels
+from scilpy.image.labels import get_data_as_labels, get_stats_in_label
 from scilpy.io.utils import (add_overwrite_arg,
                              add_verbose_arg,
                              assert_inputs_exist)
 
 
 def _build_arg_parser():
-    p = argparse.ArgumentParser(
-        description=__doc__,
-        formatter_class=argparse.RawTextHelpFormatter)
+    p = argparse.ArgumentParser(description=__doc__,
+                                formatter_class=argparse.RawTextHelpFormatter)
     p.add_argument('in_labels',
                    help='Path of the input label file.')
     p.add_argument('in_labels_lut',
                    help='Path of the LUT file corresponding to labels,'
                         'used to name the regions of interest.')
-    p.add_argument('in_seed_maps',
-                   help='Path of the input seed map file.')
+    p.add_argument('in_map',
+                   help='Path of the input map file. Expecting a 3D file.')
 
     add_verbose_arg(p)
     add_overwrite_arg(p)
@@ -45,42 +45,21 @@ def main():
     args = parser.parse_args()
     logging.getLogger().setLevel(logging.getLevelName(args.verbose))
 
-    required = args.in_labels, args.in_seed_maps, args.in_labels_lut
-    assert_inputs_exist(parser, required)
-
-    # Load atlas image
-    label_img = nib.load(args.in_labels)
-    label_img_data = get_data_as_labels(label_img)
+    # Verifications
+    assert_inputs_exist(parser,
+                        [args.in_labels, args.in_map, args.in_labels_lut])
+    # toDo after PR937
+    # assert_headers_compatible(parser, [args.in_labels, args.in_map])
 
-    # Load atlas lut
+    # Loading
+    label_data = get_data_as_labels(nib.load(args.in_labels))
     with open(args.in_labels_lut) as f:
         label_dict = json.load(f)
-    (label_indices, label_names) = zip(*label_dict.items())
-
-    # Load seed image
-    seed_img = nib.load(args.in_seed_maps)
-    seed_img_data = seed_img.get_fdata(dtype=np.float32)
-
-    for label, name in zip(label_indices, label_names):
-        label = int(label)
-        if label != 0:
-            curr_data = (seed_img_data[np.where(label_img_data == label)])
-            nb_vx_roi = np.count_nonzero(label_img_data == label)
-            nb_seed_vx = np.count_nonzero(curr_data)
-
-            if nb_seed_vx != 0:
-                mean_seed = np.sum(curr_data)/nb_seed_vx
-                max_seed = np.max(curr_data)
-                std_seed = np.sqrt(np.mean(abs(curr_data[curr_data != 0] -
-                                               mean_seed)**2))
-
-                print(json.dumps({'ROI-idx': label,
-                                  'ROI-name': str(name),
-                                  'nb-vx-roi': int(nb_vx_roi),
-                                  'nb-vx-seed': int(nb_seed_vx),
-                                  'max': int(max_seed),
-                                  'mean': float(mean_seed),
-                                  'std': float(std_seed)}))
+    map_data = nib.load(args.in_map).get_fdata(dtype=np.float32)
+
+    # Process
+    out_dict = get_stats_in_label(map_data, label_data, label_dict)
+    print(json.dumps(out_dict))
 
 
 if __name__ == "__main__":

From aa7dd257bec3cf640e85b87fd78fde1a2331d50d Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Thu, 14 Mar 2024 15:13:50 -0400
Subject: [PATCH 10/14] volume_stats_in_labels: add execution test

---
 scripts/tests/test_volume_stats_in_labels.py | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/scripts/tests/test_volume_stats_in_labels.py b/scripts/tests/test_volume_stats_in_labels.py
index e0e1c4bd8..6190f3d46 100644
--- a/scripts/tests/test_volume_stats_in_labels.py
+++ b/scripts/tests/test_volume_stats_in_labels.py
@@ -1,7 +1,19 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
+import os
+
+from scilpy import SCILPY_HOME
 
 
 def test_help_option(script_runner):
     ret = script_runner.run('scil_volume_stats_in_labels.py', '--help')
     assert ret.success
+
+
+def test_execution(script_runner):
+    in_map = os.path.join(SCILPY_HOME, 'plot', 'fa.nii.gz')
+    in_atlas = os.path.join(SCILPY_HOME, 'plot', 'atlas_brainnetome.nii.gz')
+    atlas_lut = os.path.join(SCILPY_HOME, 'plot', 'atlas_brainnetome.json')
+    ret = script_runner.run('scil_volume_stats_in_labels.py',
+                            in_atlas, atlas_lut, in_map)
+    assert ret.success

From bdd226a63cc0c19bf5d01bc11f507231967dddd9 Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Thu, 14 Mar 2024 16:23:17 -0400
Subject: [PATCH 11/14] volume_stats_in_roi: simplify

---
 scilpy/utils/metrics_tools.py          |  27 -----
 scripts/scil_volume_stats_in_ROI.py    | 147 +++++++++++--------------
 scripts/scil_volume_stats_in_labels.py |   2 +
 3 files changed, 68 insertions(+), 108 deletions(-)

diff --git a/scilpy/utils/metrics_tools.py b/scilpy/utils/metrics_tools.py
index 66860549d..d4e5be865 100644
--- a/scilpy/utils/metrics_tools.py
+++ b/scilpy/utils/metrics_tools.py
@@ -377,30 +377,3 @@ def plot_metrics_stats(means, stds, title=None, xlabel=None,
 
     plt.close(fig)
     return fig
-
-
-def get_roi_metrics_mean_std(density_map, metrics_files):
-    """
-    Returns the mean and standard deviation of each metric, using the
-    provided density map. This can be a binary mask,
-    or contain weighted values between 0 and 1.
-
-    Parameters
-    ------------
-    density_map : ndarray
-        3D numpy array containing a density map.
-    metrics_files : sequence
-        list of nibabel objects representing the metrics files.
-
-    Returns
-    ---------
-    stats : list
-        list of tuples where the first element of the tuple is the mean
-        of a metric, and the second element is the standard deviation.
-
-    """
-
-    return map(lambda metric_file:
-               weighted_mean_std(density_map,
-                                 metric_file.get_fdata(dtype=np.float64)),
-               metrics_files)
diff --git a/scripts/scil_volume_stats_in_ROI.py b/scripts/scil_volume_stats_in_ROI.py
index 9de905a4c..28a66a9b0 100755
--- a/scripts/scil_volume_stats_in_ROI.py
+++ b/scripts/scil_volume_stats_in_ROI.py
@@ -3,14 +3,12 @@
 
 """
 Compute the statistics (mean, std) of scalar maps, which can represent
-diffusion metrics, in a ROI.
+diffusion metrics, in a ROI. Prints the results.
 
 The mask can either be a binary mask, or a weighting mask. If the mask is
 a weighting mask it should either contain floats between 0 and 1 or should be
-normalized with --normalize_weights.
-
-IMPORTANT: if the mask contains weights (and not 0 and 1 exclusively), the
-standard deviation will also be weighted.
+normalized with --normalize_weights. IMPORTANT: if the mask contains weights
+(and not 0 and 1 exclusively), the standard deviation will also be weighted.
 """
 
 import argparse
@@ -28,7 +26,7 @@
                              assert_inputs_exist,
                              assert_headers_compatible)
 from scilpy.utils.filenames import split_name_with_nii
-from scilpy.utils.metrics_tools import get_roi_metrics_mean_std
+from scilpy.utils.metrics_tools import weighted_mean_std
 
 
 def _build_arg_parser():
@@ -39,23 +37,23 @@ def _build_arg_parser():
                    help='Mask volume filename.\nCan be a binary mask or a '
                         'weighted mask.')
 
-    p_metric = p.add_argument_group('Metrics input options')
-    g_metric = p_metric.add_mutually_exclusive_group(required=True)
-    g_metric.add_argument('--metrics_dir',
-                          help='Metrics files directory. Name of the '
-                               'directory containing the metrics files.')
-    g_metric.add_argument('--metrics', dest='metrics_file_list', nargs='+',
-                          help='Metrics nifti filename. List of the names of '
-                               'the metrics file, in nifti format.')
+    g = p.add_argument_group('Metrics input options')
+    gg = g.add_mutually_exclusive_group(required=True)
+    gg.add_argument('--metrics_dir', metavar='dir',
+                    help='Name of the directory containing metrics files: '
+                         'we will \nload all nifti files.')
+    gg.add_argument('--metrics', dest='metrics_file_list', nargs='+',
+                    metavar='file',
+                    help='Metrics nifti filename. List of the names of the '
+                         'metrics file, \nin nifti format.')
 
     p.add_argument('--bin', action='store_true',
-                   help='If set, will consider every value of the mask '
-                        'higher than 0 to be part of the mask, and set to 1 '
+                   help='If set, we will consider every value of the mask '
+                        'higher than 0 \nto be part of the mask '
                         '(equivalent weighting for every voxel).')
-
     p.add_argument('--normalize_weights', action='store_true',
                    help='If set, the weights will be normalized to the [0,1] '
-                        'range.')
+                        'range. \n(We divide by the max,  ')
 
     add_json_args(p)
     add_verbose_arg(p)
@@ -69,81 +67,68 @@ def main():
     args = parser.parse_args()
     logging.getLogger().setLevel(logging.getLevelName(args.verbose))
 
-    if args.metrics_dir and os.path.exists(args.metrics_dir):
-        list_metrics_files = glob.glob(os.path.join(args.metrics_dir,
-                                                    '*nii.gz'))
-        assert_inputs_exist(parser, [args.in_mask] + list_metrics_files)
-        assert_headers_compatible(parser, [args.in_mask] + list_metrics_files)
-    elif args.metrics_file_list:
-        assert_inputs_exist(parser, [args.in_mask] + args.metrics_file_list)
-        assert_headers_compatible(parser,
-                                  [args.in_mask] + args.metrics_file_list)
-
-    # Load mask and validate content depending on flags
-    mask_img = nib.load(args.in_mask)
+    # Verifications
 
-    if len(mask_img.shape) > 3:
-        logging.error('Mask should be a 3D image.')
-
-    # Can be a weighted image
-    mask_data = mask_img.get_fdata(dtype=np.float32)
+    # Get input list. Either all files in dir or given list.
+    if args.metrics_dir:
+        if not os.path.exists(args.metrics_dir):
+            parser.error("Metrics directory does not exist: {}"
+                         .format(args.metrics_dir))
+        assert_inputs_exist(parser, args.in_mask)
+
+        tmp_file_list = glob.glob(os.path.join(args.metrics_dir, '*nii.gz'))
+        args.metrics_file_list = [os.path.join(args.metrics_dir, f)
+                                  for f in tmp_file_list]
+    else:
+        assert_inputs_exist(parser, [args.in_mask] + args.metrics_file_list)
+    assert_headers_compatible(parser,
+                              [args.in_mask] + args.metrics_file_list)
 
+    # Loading
+    mask_data = nib.load(args.in_mask).get_fdata(dtype=np.float32)
+    if len(mask_data.shape) > 3:
+        parser.error('Mask should be a 3D image.')
     if np.min(mask_data) < 0:
-        logging.error('Mask should not contain negative values.')
+        parser.error('Mask should not contain negative values.')
+    roi_name = split_name_with_nii(os.path.basename(args.in_mask))[0]
 
     # Discussion about the way the normalization is done.
     # https://github.com/scilus/scilpy/pull/202#discussion_r411355609
+    # Summary:
+    # 1) We don't want to normalize with data = (data-min) / (max-min) because
+    # it zeroes out the minimal values of the array. This is not a large error
+    # source, but not preferable.
+    # 2) data = data / max(data) or data = data / sum(data): in practice, when
+    # we use them in numpy using their weights argument, leads to the same
+    # result.
     if args.normalize_weights:
         mask_data /= np.max(mask_data)
-
-    if np.min(mask_data) < 0.0 or np.max(mask_data) > 1.0:
+    elif args.bin:
+        mask_data[np.where(mask_data > 0.0)] = 1.0
+    elif np.min(mask_data) < 0.0 or np.max(mask_data) > 1.0:
         parser.error('Mask data should only contain values between 0 and 1. '
                      'Try --normalize_weights.')
 
-    if args.bin:
-        mask_data[np.where(mask_data > 0.0)] = 1.0
-
-    # Load all metrics files.
-    metrics_files = []
-    if args.metrics_dir:
-        for f in sorted(os.listdir(args.metrics_dir)):
-            metric_img = nib.load(os.path.join(args.metrics_dir, f))
-            if len(metric_img.shape)==3:
-                # Check if NaNs in metrics
-                if np.any(np.isnan(metric_img.get_fdata(dtype=np.float64))):
-                    logging.warning('Metric \"{}\" contains some NaN.'.format(metric_img.get_filename()) +
-                                    ' Ignoring voxels with NaN.')
-                metrics_files.append(metric_img)
-            else:
-                parser.error('Metric {} is not compatible ({}D image).'.format(os.path.join(args.metrics_dir, f),
-                                                                               len(metric_img.shape)))
-    elif args.metrics_file_list:
-        assert_headers_compatible(parser, [args.in_mask] +
-                                  args.metrics_file_list)
-        for f in args.metrics_file_list:
-            metric_img = nib.load(f)
-            if len(metric_img.shape)==3:
-                # Check if NaNs in metrics
-                if np.any(np.isnan(metric_img.get_fdata(dtype=np.float64))):
-                    logging.warning('Metric \"{}\" contains some NaN.'.format(metric_img.get_filename()) +
-                                    ' Ignoring voxels with NaN.')
-                metrics_files.append(metric_img)
-            else:
-                parser.error('Metric {} is not compatible ({}D image).'.format(f,
-                                                                               len(metric_img.shape)))
-    # Compute the mean values and standard deviations
-    stats = get_roi_metrics_mean_std(mask_data, metrics_files)
-
-    roi_name = split_name_with_nii(os.path.basename(args.in_mask))[0]
+    # Load and process all metrics files.
     json_stats = {roi_name: {}}
-    for metric_file, (mean, std) in zip(metrics_files, stats):
-        metric_name = split_name_with_nii(
-            os.path.basename(metric_file.get_filename()))[0]
-        json_stats[roi_name][metric_name] = {
-            'mean': mean.item(),
-            'std': std.item()
-        }
-
+    for f in args.metrics_file_list:
+        metric_img = nib.load(f)
+        metric_name = split_name_with_nii(os.path.basename(f))[0]
+        if len(metric_img.shape) == 3:
+            data = metric_img.get_fdata(dtype=np.float64)
+            if np.any(np.isnan(data)):
+                logging.warning("Metric '{}' contains some NaN. Ignoring "
+                                "voxels with NaN."
+                                .format(os.path.basename(f)))
+            mean, std = weighted_mean_std(mask_data, data)
+            json_stats[roi_name][metric_name] = {'mean': mean,
+                                                 'std': std}
+        else:
+            parser.error(
+                'Metric {} is not a 3D image ({}D shape).'
+                .format(f, len(metric_img.shape)))
+
+    # Print results
     print(json.dumps(json_stats, indent=args.indent, sort_keys=args.sort_keys))
 
 
diff --git a/scripts/scil_volume_stats_in_labels.py b/scripts/scil_volume_stats_in_labels.py
index ebefde0fa..8fad012f7 100755
--- a/scripts/scil_volume_stats_in_labels.py
+++ b/scripts/scil_volume_stats_in_labels.py
@@ -56,6 +56,8 @@ def main():
     with open(args.in_labels_lut) as f:
         label_dict = json.load(f)
     map_data = nib.load(args.in_map).get_fdata(dtype=np.float32)
+    if len(map_data.shape) > 3:
+        parser.error('Mask should be a 3D image.')
 
     # Process
     out_dict = get_stats_in_label(map_data, label_data, label_dict)

From d050abcd223d6acb06322a2c93634218d8f275ce Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Thu, 14 Mar 2024 16:26:23 -0400
Subject: [PATCH 12/14] Fix failed crop test

---
 scripts/tests/test_volume_crop.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/scripts/tests/test_volume_crop.py b/scripts/tests/test_volume_crop.py
index 13890880b..145cd528b 100644
--- a/scripts/tests/test_volume_crop.py
+++ b/scripts/tests/test_volume_crop.py
@@ -25,6 +25,6 @@ def test_execution_processing(script_runner, monkeypatch):
     assert ret.success
 
     # Then try to load back the same box
-    ret = script_runner.run('scil_crop_volume.py', in_dwi, 'dwi_crop.nii.gz',
+    ret = script_runner.run('scil_crop_volume.py', in_dwi, 'dwi_crop2.nii.gz',
                             '--input_bbox', 'bbox.pickle')
     assert ret.success

From af04fc49dd67dfcec28e3bc2430d0a3ec5193ce2 Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Mon, 18 Mar 2024 10:34:23 -0400
Subject: [PATCH 13/14] Use new asser_headers_compatible in volumes scripts

---
 scripts/scil_volume_resample.py        | 6 +++---
 scripts/scil_volume_stats_in_labels.py | 8 +++-----
 2 files changed, 6 insertions(+), 8 deletions(-)

diff --git a/scripts/scil_volume_resample.py b/scripts/scil_volume_resample.py
index fbe241053..e65e5b853 100755
--- a/scripts/scil_volume_resample.py
+++ b/scripts/scil_volume_resample.py
@@ -14,7 +14,8 @@
 import nibabel as nib
 
 from scilpy.io.utils import (add_verbose_arg, add_overwrite_arg,
-                             assert_inputs_exist, assert_outputs_exist)
+                             assert_inputs_exist, assert_outputs_exist,
+                             assert_headers_compatible)
 from scilpy.image.volume_operations import resample_volume
 
 
@@ -66,8 +67,7 @@ def main():
     # Checking args
     assert_inputs_exist(parser, args.in_image, args.ref)
     assert_outputs_exist(parser, args, args.out_image)
-    # toDo after PR937
-    #assert_headers_compatible(parser, args.in_image, args.ref)
+    assert_headers_compatible(parser, args.in_image, args.ref)
 
     if args.enforce_dimensions and not args.ref:
         parser.error("Cannot enforce dimensions without a reference image")
diff --git a/scripts/scil_volume_stats_in_labels.py b/scripts/scil_volume_stats_in_labels.py
index 8fad012f7..07fc290b7 100755
--- a/scripts/scil_volume_stats_in_labels.py
+++ b/scripts/scil_volume_stats_in_labels.py
@@ -18,9 +18,8 @@
 import numpy as np
 
 from scilpy.image.labels import get_data_as_labels, get_stats_in_label
-from scilpy.io.utils import (add_overwrite_arg,
-                             add_verbose_arg,
-                             assert_inputs_exist)
+from scilpy.io.utils import (add_overwrite_arg, add_verbose_arg,
+                             assert_inputs_exist, assert_headers_compatible)
 
 
 def _build_arg_parser():
@@ -48,8 +47,7 @@ def main():
     # Verifications
     assert_inputs_exist(parser,
                         [args.in_labels, args.in_map, args.in_labels_lut])
-    # toDo after PR937
-    # assert_headers_compatible(parser, [args.in_labels, args.in_map])
+    assert_headers_compatible(parser, [args.in_labels, args.in_map])
 
     # Loading
     label_data = get_data_as_labels(nib.load(args.in_labels))

From ef3859a997bd204a5ec1a997e17bd0e22854d433 Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Thu, 28 Mar 2024 13:36:06 -0400
Subject: [PATCH 14/14] Answer Phil's comment

---
 scripts/scil_volume_stats_in_ROI.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/scripts/scil_volume_stats_in_ROI.py b/scripts/scil_volume_stats_in_ROI.py
index 28a66a9b0..831373fdc 100755
--- a/scripts/scil_volume_stats_in_ROI.py
+++ b/scripts/scil_volume_stats_in_ROI.py
@@ -48,12 +48,12 @@ def _build_arg_parser():
                          'metrics file, \nin nifti format.')
 
     p.add_argument('--bin', action='store_true',
-                   help='If set, we will consider every value of the mask '
-                        'higher than 0 \nto be part of the mask '
-                        '(equivalent weighting for every voxel).')
+                   help='If set, will consider every value of the mask higher'
+                        'than 0 to be \npart of the mask (equivalent '
+                        'weighting for every voxel).')
     p.add_argument('--normalize_weights', action='store_true',
                    help='If set, the weights will be normalized to the [0,1] '
-                        'range. \n(We divide by the max,  ')
+                        'range.')
 
     add_json_args(p)
     add_verbose_arg(p)