Preserve data dtype in concat_images to fix save/load round-trip (gh-986)

nibabel/funcs.py

@@ -115,13 +115,10 @@ def concat_images(images, check_affines=True, axis=None):
     klass = img0.__class__
     shape0 = img0.shape
     n_dim = len(shape0)
-    if axis is None:
-        # collect images in output array for efficiency
-        out_shape = (n_imgs,) + shape0
-        out_data = np.empty(out_shape)
-    else:
-        # collect images in list for use with np.concatenate
-        out_data = [None] * n_imgs
+    # collect images in a list; stacking/concatenating later preserves the input
+    # data dtype instead of upcasting to float64, so an integer image round-trips
+    # through save/load without precision loss (gh-986)
+    out_data = [None] * n_imgs
     # Get part of shape we need to check inside loop
     idx_mask = np.ones((n_dim,), dtype=bool)
     if axis is not None:
@@ -141,7 +138,7 @@ def concat_images(images, check_affines=True, axis=None):
         out_data[i] = np.asanyarray(img.dataobj)
 
     if axis is None:
-        out_data = np.rollaxis(out_data, 0, out_data.ndim)
+        out_data = np.stack(out_data, axis=-1)
     else:
         out_data = np.concatenate(out_data, axis=axis)
 

nibabel/tests/test_funcs.py

@@ -14,7 +14,7 @@
 
 from ..analyze import AnalyzeImage
 from ..funcs import OrientationError, as_closest_canonical, concat_images
-from ..loadsave import save
+from ..loadsave import load, save
 from ..nifti1 import Nifti1Image
 from ..tmpdirs import InTemporaryDirectory
 
@@ -192,3 +192,27 @@ def test_closest_canonical():
     img.header.set_dim_info(None, None, 2)
     xyz_img = as_closest_canonical(img)
     assert xyz_img.header.get_dim_info() == (None, None, 1)
+
+
+def test_concat_integer_roundtrip():
+    # Regression test for gh-986: concatenating integer images with the default
+    # ``axis=None`` must preserve the input dtype rather than upcasting to float64,
+    # so the saved-and-reloaded data equals the in-memory data instead of being
+    # quantized back to the integer storage dtype on save.
+    arr0 = np.arange(24, dtype=np.uint16).reshape(2, 3, 4)
+    arr1 = arr0 + 1000
+    img0 = Nifti1Image(arr0, np.eye(4))
+    img1 = Nifti1Image(arr1, np.eye(4))
+
+    concat = concat_images([img0, img1])
+    assert concat.get_data_dtype() == np.dtype(np.uint16)
+    in_memory = np.asarray(concat.dataobj)
+    assert in_memory.dtype == np.uint16
+
+    with InTemporaryDirectory():
+        save(concat, 'concat.nii')
+        # Load with mmap=False so the file handle is released before the
+        # temporary directory is torn down; Windows cannot delete a file
+        # that is still memory-mapped, which otherwise fails cleanup here.
+        reloaded = np.asarray(load('concat.nii', mmap=False).dataobj)
+    assert_array_equal(in_memory, reloaded)