Merge pull request #1082 from effigies/enh/image_dtype_arg

NF: Pass dtype to Analyze-like images at initialization/serialization, warn on creation of NIfTI images with 64-bit ints (API change)

Author: effigies

nibabel/analyze.py

@@ -914,12 +914,15 @@ class AnalyzeImage(SpatialImage):
     ImageArrayProxy = ArrayProxy
 
     def __init__(self, dataobj, affine, header=None,
-                 extra=None, file_map=None):
+                 extra=None, file_map=None, dtype=None):
         super(AnalyzeImage, self).__init__(
             dataobj, affine, header, extra, file_map)
         # Reset consumable values
         self._header.set_data_offset(0)
         self._header.set_slope_inter(None, None)
+
+        if dtype is not None:
+            self.set_data_dtype(dtype)
     __init__.__doc__ = SpatialImage.__init__.__doc__
 
     def get_data_dtype(self):
@@ -989,23 +992,29 @@ def _get_fileholders(file_map):
         """
         return file_map['header'], file_map['image']
 
-    def to_file_map(self, file_map=None):
+    def to_file_map(self, file_map=None, dtype=None):
         """ Write image to `file_map` or contained ``self.file_map``
 
         Parameters
         ----------
         file_map : None or mapping, optional
            files mapping.  If None (default) use object's ``file_map``
            attribute instead
+        dtype : dtype-like, optional
+           The on-disk data type to coerce the data array.
         """
         if file_map is None:
             file_map = self.file_map
         data = np.asanyarray(self.dataobj)
         self.update_header()
         hdr = self._header
-        out_dtype = self.get_data_dtype()
         # Store consumable values for later restore
         offset = hdr.get_data_offset()
+        data_dtype = hdr.get_data_dtype()
+        # Override dtype conditionally
+        if dtype is not None:
+            hdr.set_data_dtype(dtype)
+        out_dtype = hdr.get_data_dtype()
         # Scalars of slope, offset to get immutable values
         slope = hdr['scl_slope'].item() if hdr.has_data_slope else np.nan
         inter = hdr['scl_inter'].item() if hdr.has_data_intercept else np.nan
@@ -1045,6 +1054,7 @@ def to_file_map(self, file_map=None):
         self.file_map = file_map
         # Restore any changed consumable values
         hdr.set_data_offset(offset)
+        hdr.set_data_dtype(data_dtype)
         if hdr.has_data_slope:
             hdr['scl_slope'] = slope
         if hdr.has_data_intercept:

nibabel/arrayproxy.py

@@ -412,3 +412,17 @@ def reshape_dataobj(obj, shape):
     """
     return (obj.reshape(shape) if hasattr(obj, 'reshape')
             else np.reshape(obj, shape))
+
+
+def get_obj_dtype(obj):
+    """ Get the effective dtype of an array-like object """
+    if is_proxy(obj):
+        # Read and potentially apply scaling to one value
+        idx = (0,) * len(obj.shape)
+        return obj[idx].dtype
+    elif hasattr(obj, "dtype"):
+        # Trust the dtype (probably an ndarray)
+        return obj.dtype
+    else:
+        # Coerce; this could be expensive but we don't know what we can do with it
+        return np.asanyarray(obj).dtype

nibabel/cifti2/cifti2.py

@@ -1460,7 +1460,7 @@ def from_image(klass, img):
             return img
         raise NotImplementedError
 
-    def to_file_map(self, file_map=None):
+    def to_file_map(self, file_map=None, dtype=None):
         """ Write image to `file_map` or contained ``self.file_map``
 
         Parameters
@@ -1493,7 +1493,7 @@ def to_file_map(self, file_map=None):
         # If qform not set, reset pixdim values so Nifti2 does not complain
         if header['qform_code'] == 0:
             header['pixdim'][:4] = 1
-        img = Nifti2Image(data, None, header)
+        img = Nifti2Image(data, None, header, dtype=dtype)
         img.to_file_map(file_map or self.file_map)
 
     def update_headers(self):

nibabel/filebasedimages.py

@@ -299,24 +299,26 @@ def filespec_to_file_map(klass, filespec):
             file_map[key] = FileHolder(filename=fname)
         return file_map
 
-    def to_filename(self, filename):
-        """ Write image to files implied by filename string
+    def to_filename(self, filename, **kwargs):
+        r""" Write image to files implied by filename string
 
         Parameters
         ----------
         filename : str or os.PathLike
            filename to which to save image.  We will parse `filename`
            with ``filespec_to_file_map`` to work out names for image,
            header etc.
+        \*\*kwargs : keyword arguments
+           Keyword arguments to format-specific save
 
         Returns
         -------
         None
         """
         self.file_map = self.filespec_to_file_map(filename)
-        self.to_file_map()
+        self.to_file_map(**kwargs)
 
-    def to_file_map(self, file_map=None):
+    def to_file_map(self, file_map=None, **kwargs):
         raise NotImplementedError
 
     @classmethod
@@ -552,13 +554,14 @@ def from_bytes(klass, bytestring):
         file_map = klass.make_file_map({'image': bio, 'header': bio})
         return klass.from_file_map(file_map)
 
-    def to_bytes(self):
-        """ Return a ``bytes`` object with the contents of the file that would
+    def to_bytes(self, **kwargs):
+        r""" Return a ``bytes`` object with the contents of the file that would
         be written if the image were saved.
 
         Parameters
         ----------
-        None
+        \*\*kwargs : keyword arguments
+            Keyword arguments that may be passed to ``img.to_file_map()``
 
         Returns
         -------
@@ -569,5 +572,5 @@ def to_bytes(self):
             raise NotImplementedError("to_bytes() is undefined for multi-file images")
         bio = io.BytesIO()
         file_map = self.make_file_map({'image': bio, 'header': bio})
-        self.to_file_map(file_map)
+        self.to_file_map(file_map, **kwargs)
         return bio.getvalue()

nibabel/funcs.py

@@ -35,7 +35,7 @@ def squeeze_image(img):
     --------
     >>> import nibabel as nf
     >>> shape = (10,20,30,1,1)
-    >>> data = np.arange(np.prod(shape)).reshape(shape)
+    >>> data = np.arange(np.prod(shape), dtype='int32').reshape(shape)
     >>> affine = np.eye(4)
     >>> img = nf.Nifti1Image(data, affine)
     >>> img.shape == (10, 20, 30, 1, 1)
@@ -47,7 +47,7 @@ def squeeze_image(img):
     If the data are 3D then last dimensions of 1 are ignored
 
     >>> shape = (10,1,1)
-    >>> data = np.arange(np.prod(shape)).reshape(shape)
+    >>> data = np.arange(np.prod(shape), dtype='int32').reshape(shape)
     >>> img = nf.ni1.Nifti1Image(data, affine)
     >>> img.shape == (10, 1, 1)
     True

nibabel/loadsave.py

@@ -130,15 +130,17 @@ def guessed_image_type(filename):
     raise ImageFileError(f'Cannot work out file type of "{filename}"')
 
 
-def save(img, filename):
-    """ Save an image to file adapting format to `filename`
+def save(img, filename, **kwargs):
+    r""" Save an image to file adapting format to `filename`
 
     Parameters
     ----------
     img : ``SpatialImage``
        image to save
     filename : str or os.PathLike
        filename (often implying filenames) to which to save `img`.
+    \*\*kwargs : keyword arguments
+        Keyword arguments to format-specific save
 
     Returns
     -------
@@ -148,7 +150,7 @@ def save(img, filename):
 
     # Save the type as expected
     try:
-        img.to_filename(filename)
+        img.to_filename(filename, **kwargs)
     except ImageFileError:
         pass
     else:
@@ -196,7 +198,7 @@ def save(img, filename):
     # Here, we either have a klass or a converted image.
     if converted is None:
         converted = klass.from_image(img)
-    converted.to_filename(filename)
+    converted.to_filename(filename, **kwargs)
 
 
 @deprecate_with_version('read_img_data deprecated. '

nibabel/nifti1.py

@@ -17,6 +17,7 @@
 import numpy.linalg as npl
 from numpy.compat.py3k import asstr
 
+from .arrayproxy import get_obj_dtype
 from .optpkg import optional_package
 from .filebasedimages import SerializableImage
 from .volumeutils import Recoder, make_dt_codes, endian_codes
@@ -881,6 +882,51 @@ def set_data_shape(self, shape):
             shape = (-1, 1, 1) + shape[3:]
         super(Nifti1Header, self).set_data_shape(shape)
 
+    def set_data_dtype(self, datatype):
+        """ Set numpy dtype for data from code or dtype or type
+
+        Using :py:class:`int` or ``"int"`` is disallowed, as these types
+        will be interpreted as ``np.int64``, which is almost never desired.
+        ``np.int64`` is permitted for those intent on making poor choices.
+
+        Examples
+        --------
+        >>> hdr = Nifti1Header()
+        >>> hdr.set_data_dtype(np.uint8)
+        >>> hdr.get_data_dtype()
+        dtype('uint8')
+        >>> hdr.set_data_dtype(np.dtype(np.uint8))
+        >>> hdr.get_data_dtype()
+        dtype('uint8')
+        >>> hdr.set_data_dtype('implausible') #doctest: +IGNORE_EXCEPTION_DETAIL
+        Traceback (most recent call last):
+           ...
+        HeaderDataError: data dtype "implausible" not recognized
+        >>> hdr.set_data_dtype('none') #doctest: +IGNORE_EXCEPTION_DETAIL
+        Traceback (most recent call last):
+           ...
+        HeaderDataError: data dtype "none" known but not supported
+        >>> hdr.set_data_dtype(np.void) #doctest: +IGNORE_EXCEPTION_DETAIL
+        Traceback (most recent call last):
+           ...
+        HeaderDataError: data dtype "<type 'numpy.void'>" known but not supported
+        >>> hdr.set_data_dtype('int') #doctest: +IGNORE_EXCEPTION_DETAIL
+        Traceback (most recent call last):
+           ...
+        ValueError: Invalid data type 'int'. Specify a sized integer, e.g., 'uint8' or numpy.int16.
+        >>> hdr.set_data_dtype(int) #doctest: +IGNORE_EXCEPTION_DETAIL
+        Traceback (most recent call last):
+           ...
+        ValueError: Invalid data type 'int'. Specify a sized integer, e.g., 'uint8' or numpy.int16.
+        >>> hdr.set_data_dtype('int64')
+        >>> hdr.get_data_dtype() == np.dtype('int64')
+        True
+        """
+        if not isinstance(datatype, np.dtype) and datatype in (int, "int"):
+            raise ValueError(f"Invalid data type {datatype!r}. Specify a sized integer, "
+                             "e.g., 'uint8' or numpy.int16.")
+        super().set_data_dtype(datatype)
+
     def get_qform_quaternion(self):
         """ Compute quaternion from b, c, d of quaternion
 
@@ -1754,12 +1800,27 @@ class Nifti1Pair(analyze.AnalyzeImage):
     rw = True
 
     def __init__(self, dataobj, affine, header=None,
-                 extra=None, file_map=None):
+                 extra=None, file_map=None, dtype=None):
+        # Special carve-out for 64 bit integers
+        # See GitHub issues
+        #  * https://github.com/nipy/nibabel/issues/1046
+        #  * https://github.com/nipy/nibabel/issues/1089
+        # This only applies to NIfTI because the parent Analyze formats did
+        # not support 64-bit integer data, so `set_data_dtype(int64)` would
+        # already fail.
+        danger_dts = (np.dtype("int64"), np.dtype("uint64"))
+        if header is None and dtype is None and get_obj_dtype(dataobj) in danger_dts:
+            msg = (f"Image data has type {dataobj.dtype}, which may cause "
+                   "incompatibilities with other tools. This will error in "
+                   "NiBabel 5.0. This warning can be silenced "
+                   f"by passing the dtype argument to {self.__class__.__name__}().")
+            warnings.warn(msg, FutureWarning, stacklevel=2)
         super(Nifti1Pair, self).__init__(dataobj,
                                          affine,
                                          header,
                                          extra,
-                                         file_map)
+                                         file_map,
+                                         dtype)
         # Force set of s/q form when header is None unless affine is also None
         if header is None and affine is not None:
             self._affine2header()
@@ -1865,7 +1926,7 @@ def set_qform(self, affine, code=None, strip_shears=True, **kwargs):
 
         Examples
         --------
-        >>> data = np.arange(24).reshape((2,3,4))
+        >>> data = np.arange(24, dtype='f4').reshape((2,3,4))
         >>> aff = np.diag([2, 3, 4, 1])
         >>> img = Nifti1Pair(data, aff)
         >>> img.get_qform()
@@ -1948,7 +2009,7 @@ def set_sform(self, affine, code=None, **kwargs):
 
         Examples
         --------
-        >>> data = np.arange(24).reshape((2,3,4))
+        >>> data = np.arange(24, dtype='f4').reshape((2,3,4))
         >>> aff = np.diag([2, 3, 4, 1])
         >>> img = Nifti1Pair(data, aff)
         >>> img.get_sform()

nibabel/spm99analyze.py

@@ -308,7 +308,7 @@ def from_file_map(klass, file_map, *, mmap=True, keep_file_open=None):
         ret._affine = np.dot(ret._affine, to_111)
         return ret
 
-    def to_file_map(self, file_map=None):
+    def to_file_map(self, file_map=None, dtype=None):
         """ Write image to `file_map` or contained ``self.file_map``
 
         Extends Analyze ``to_file_map`` method by writing ``mat`` file
@@ -321,7 +321,7 @@ def to_file_map(self, file_map=None):
         """
         if file_map is None:
             file_map = self.file_map
-        super(Spm99AnalyzeImage, self).to_file_map(file_map)
+        super(Spm99AnalyzeImage, self).to_file_map(file_map, dtype=dtype)
         mat = self._affine
         if mat is None:
             return

nibabel/tests/test_analyze.py

@@ -290,7 +290,12 @@ def assert_set_dtype(dt_spec, np_dtype):
         # Test aliases to Python types
         assert_set_dtype(float, np.float64)  # float64 always supported
         np_sys_int = np.dtype(int).type  # int could be 32 or 64 bit
-        if np_sys_int in self.supported_np_types:  # no int64 for Analyze
+        if issubclass(self.header_class, Nifti1Header):
+            # We don't allow int aliases in Nifti
+            with pytest.raises(ValueError):
+                hdr = self.header_class()
+                hdr.set_data_dtype(int)
+        elif np_sys_int in self.supported_np_types:  # no int64 for Analyze
             assert_set_dtype(int, np_sys_int)
         hdr = self.header_class()
         for inp in all_unsupported_types:
@@ -759,6 +764,20 @@ def test_affine_44(self):
         with pytest.raises(ValueError):
             IC(data, np.diag([2, 3, 4]))
 
+    def test_dtype_init_arg(self):
+        # data_dtype can be set by argument in absence of header
+        img_klass = self.image_class
+        arr = np.arange(24, dtype=np.int16).reshape((2, 3, 4))
+        aff = np.eye(4)
+        for dtype in self.supported_np_types:
+            img = img_klass(arr, aff, dtype=dtype)
+            assert img.get_data_dtype() == dtype
+        # It can also override the header dtype
+        hdr = img.header
+        for dtype in self.supported_np_types:
+            img = img_klass(arr, aff, hdr, dtype=dtype)
+            assert img.get_data_dtype() == dtype
+
     def test_offset_to_zero(self):
         # Check offset is always set to zero when creating images
         img_klass = self.image_class
@@ -873,6 +892,21 @@ def test_no_finite_values(self):
         img_back = self.image_class.from_file_map(fm)
         assert_array_equal(img_back.dataobj, 0)
 
+    def test_dtype_to_filename_arg(self):
+        # data_dtype can be set by argument in absence of header
+        img_klass = self.image_class
+        arr = np.arange(24, dtype=np.int16).reshape((2, 3, 4))
+        aff = np.eye(4)
+        img = img_klass(arr, aff)
+        fname = 'test' + img_klass.files_types[0][1]
+        with InTemporaryDirectory():
+            for dtype in self.supported_np_types:
+                img.to_filename(fname, dtype=dtype)
+                new_img = img_klass.from_filename(fname)
+                assert new_img.get_data_dtype() == dtype
+                # data_type is reset after write
+                assert img.get_data_dtype() == np.int16
+
 
 def test_unsupported():
     # analyze does not support uint32