Merge branch 'mnt/imcp-dont-modify-data' into 'main'

MNT: Work-around `nibabel` auto-rescaling logic so that data/slope/intercept in `imcp`/`immv`'d images are not modified

See merge request fsl/fslpy!462

CHANGELOG.rst

@@ -2,6 +2,20 @@ This document contains the ``fslpy`` release history in reverse chronological
 order.
 
 
+3.21.0 (Tuesday 23rd July 2024)
+-------------------------------
+
+
+Changed
+^^^^^^^
+
+
+* Behaviour of the ``imcp`` / ``immv`` commands has been adjusted so that
+  when an image file is copied/moved, and converted to a different format
+  (e.g. ``.nii`` to ``.nii.gz``), the image data and header slope/intercept
+  fields are not modified (!462).
+
+
 3.20.0 (Wednesday 10th July 2024)
 ---------------------------------
 

fsl/tests/__init__.py

@@ -10,6 +10,7 @@
 import              os
 import              sys
 import              glob
+import              hashlib
 import              shutil
 import              fnmatch
 import              logging
@@ -439,3 +440,10 @@ def make_random_mask(filename, shape, xform, premask=None, minones=1):
     img.save(filename)
 
     return img
+
+
+def sha256(filename):
+    hashobj = hashlib.sha256()
+    with open(filename, 'rb') as f:
+        hashobj.update(f.read())
+    return hashobj.hexdigest()

fsl/tests/test_dicom.py

@@ -151,7 +151,7 @@ def test_scanDir():
         datafile = op.join(datadir, 'example_dicom.tbz2')
 
         with tarfile.open(datafile) as f:
-            f.extractall()
+            f.extractall(filter='data')
 
         series = fsldcm.scanDir('.')
         assert len(series) == 2

fsl/tests/test_immv_imcp.py

@@ -7,22 +7,26 @@
 
 from __future__ import print_function
 
+import              gzip
+import itertools as it
+import os.path   as op
+import              os
+import              shutil
+import              tempfile
 
+from unittest import mock
 
-import os.path    as op
-import               os
-import               shutil
-import               tempfile
-
-import numpy as np
+import numpy   as np
 import nibabel as nib
 
-import fsl.utils.imcp        as imcp
-import fsl.data.image        as fslimage
+import fsl.utils.imcp    as imcp
+import fsl.utils.tempdir as tempdir
+import fsl.data.image    as fslimage
 
-from fsl.tests import make_random_image
-from fsl.tests import make_dummy_file
-from fsl.tests import looks_like_image
+from fsl.tests import (make_random_image,
+                       make_dummy_file,
+                       looks_like_image,
+                       sha256)
 
 
 real_print = print
@@ -315,3 +319,60 @@ def test_imcp_shouldPass(move=False):
 
 def test_immv_shouldPass():
     test_imcp_shouldPass(move=True)
+
+
+def test_imcp_data_unmodified():
+    """Test that the data in an imcp'd image file is not modified. """
+
+    dtypes = [
+       np.int16,
+       np.int32,
+       np.float32,
+       np.float64]
+
+    slints = [(None, None), (1, 0), (3, 1.5)]
+
+    for dtype, (slope, inter) in it.product(dtypes, slints):
+        with tempdir.tempdir():
+            data = np.random.randint(1, 100, (10, 10, 10)).astype(dtype)
+            hdr  = nib.Nifti1Header()
+            hdr.set_data_dtype(dtype)
+            hdr.set_data_shape((10, 10, 10))
+            hdr.set_slope_inter(slope, inter)
+            hdr.set_sform(np.eye(4))
+
+            # write header/data separately, as otherwise
+            # nibabel will automatically rescale the data
+            with open('image.nii', 'wb') as f:
+                hdr.write_to(f)
+                f.write(data.tobytes())
+
+            # Input/output formats the same,
+            # should induce a straight file copy
+            imcp.imcp('image.nii', 'copied.nii', useDefaultExt=False)
+
+            # uncompresed->compressed will cause imcp
+            # to load in the image, rather than doing a
+            # file copy
+            with mock.patch.dict(os.environ, FSLOUTPUTTYPE='NIFTI_GZ'):
+                imcp.imcp('image.nii', 'converted.nii.gz', useDefaultExt=True)
+
+            # copied files should be identical
+            assert sha256('image.nii') == sha256('copied.nii')
+
+            # Converted files should have the same
+            # data, slope, and intercept. Read result
+            # header/data separately to avoid nibabel
+            # auto-rescaling.
+            with gzip.open('converted.nii.gz', 'rb') as f:
+                gothdr    = nib.Nifti1Header.from_fileobj(f)
+                databytes = f.read()
+
+            gotdata = np.frombuffer(databytes, dtype=dtype).reshape((10, 10, 10))
+
+            # Data should be binary identical
+            assert np.all(gotdata == data)
+
+            if slope is None: slope = 1
+            if inter is None: inter = 0
+            assert np.all(np.isclose(gothdr.get_slope_inter(), (slope, inter)))

fsl/utils/imcp.py

@@ -50,7 +50,8 @@ def imcp(src,
                         :func:`~fsl.data.image.defaultOutputType`. If the
                         source file does not have the same type as the default
                         extension, it will be converted. If ``False``, the
-                        source file type is not changed.
+                        source file type is used, and the destination file type
+                        (if specified) is ignored.
 
     :arg move:          If ``True``, the files are moved, instead of being
                         copied. See :func:`immv`.
@@ -143,7 +144,7 @@ def imcp(src,
 
         img = nib.load(src)
 
-        # Force conversion to specific data type if
+        # Force conversion to specific file format if
         # necessary.  The file format (pair, gzipped
         # or not) is taken care of automatically by
         # nibabel
@@ -151,7 +152,26 @@ def imcp(src,
         elif 'NIFTI2'  in destType.name: cls = nib.Nifti2Image
         elif 'NIFTI'   in destType.name: cls = nib.Nifti1Image
 
-        img = cls(np.asanyarray(img.dataobj), None, header=img.header)
+        # The default behaviour of nibabel when saving
+        # is to rescale the image data to the full range
+        # of the data type, and then set the scl_slope/
+        # inter header fields accordingly. This is highly
+        # disruptive in many circumstances. Fortunately:
+        #  - The nibabel ArrayProxy class provides a
+        #    get_unscaled method, which allows us to
+        #    bypass the rescaling at load time.
+        #  - Explicitly setting the slope and intercept
+        #    on the header allows us to bypass rescaling
+        #    at save time.
+        #
+        # https://github.com/nipy/nibabel/issues/1001
+        # https://neurostars.org/t/preserve-datatype-and-precision-with-nibabel/27641/2
+        slope = img.dataobj.slope
+        inter = img.dataobj.inter
+        data  = np.asanyarray(img.dataobj.get_unscaled(),
+                              dtype=img.get_data_dtype())
+        img   = cls(data, None, header=img.header)
+        img.header.set_slope_inter(slope, inter)
         nib.save(img, dest)
 
         # Make sure the image reference is cleared, and

fsl/version.py

@@ -47,7 +47,7 @@ import            re
 import            string
 
 
-__version__ = '3.21.0.dev0'
+__version__ = '3.22.0.dev0'
 """Current version number, as a string. """