NF: Enable data scaling within the target dtype

Author: effigies

nibabel/arrayproxy.py

@@ -336,36 +336,68 @@ def _get_fileobj(self):
                     self.file_like, keep_open=False) as opener:
                 yield opener
 
+    def _get_unscaled(self, slicer):
+        with self._get_fileobj() as fileobj:
+            return fileslice(fileobj,
+                             slicer,
+                             self._shape,
+                             self._dtype,
+                             self._offset,
+                             order=self.order,
+                             lock=self._lock)
+
+    def _get_scaled(self, dtype, slicer):
+        # Ensure scale factors have dtypes
+        scl_slope = np.asanyarray(self._slope)
+        scl_inter = np.asanyarray(self._inter)
+        if dtype is None:
+            dtype = max(scl_slope.dtype, scl_inter.dtype)
+        slope = scl_slope.astype(dtype)
+        inter = scl_inter.astype(dtype)
+        # Read array
+        raw_data = self.get_unscaled(slicer=slicer)
+        # Upcast as necessary for big slopes, intercepts
+        return apply_read_scaling(raw_data, slope, inter)
+
     def get_unscaled(self):
-        """ Read of data from file
+        """ Read data from file
 
         This is an optional part of the proxy API
         """
-        with self._get_fileobj() as fileobj, self._lock:
-            raw_data = array_from_file(self._shape,
-                                       self._dtype,
-                                       fileobj,
-                                       offset=self._offset,
-                                       order=self.order,
-                                       mmap=self._mmap)
-        return raw_data
+        return self._get_unscaled(slicer=())
+
+    def get_scaled(self, dtype=None):
+        """ Read data from file and apply scaling
+
+        The dtype of the returned array is the narrowest dtype that can
+        represent the data without overflow, and is at least as wide as
+        the dtype parameter.
+
+        If dtype is unspecified, it is the wider of the dtypes of the slope
+        or intercept. This will generally be determined by the parameter
+        size in the image header, and so should be consistent for a given
+        image format, but may vary across formats. Notably, these factors
+        are single-precision (32-bit) floats for NIfTI-1 and double-precision
+        (64-bit) floats for NIfTI-2.
+
+        Parameters
+        ----------
+        dtype : numpy dtype specifier
+            A numpy dtype specifier specifying the narrowest acceptable
+            dtype.
+
+        Returns
+        -------
+        array
+            Scaled of image data of data type `dtype`.
+        """
+        return self._get_scaled(dtype=dtype, slicer=())
 
     def __array__(self):
-        # Read array and scale
-        raw_data = self.get_unscaled()
-        return apply_read_scaling(raw_data, self._slope, self._inter)
+        return self._get_scaled(dtype=None, slicer=())
 
     def __getitem__(self, slicer):
-        with self._get_fileobj() as fileobj:
-            raw_data = fileslice(fileobj,
-                                 slicer,
-                                 self._shape,
-                                 self._dtype,
-                                 self._offset,
-                                 order=self.order,
-                                 lock=self._lock)
-        # Upcast as necessary for big slopes, intercepts
-        return apply_read_scaling(raw_data, self._slope, self._inter)
+        return self._get_scaled(dtype=None, slicer=slicer)
 
     def reshape(self, shape):
         """ Return an ArrayProxy with a new shape, without modifying data """

nibabel/dataobj_images.py

@@ -10,6 +10,7 @@
 
 import numpy as np
 
+from .arrayproxy import is_proxy
 from .filebasedimages import FileBasedImage
 from .keywordonly import kw_only_meth
 from .deprecated import deprecate_with_version
@@ -350,7 +351,14 @@ def get_fdata(self, caching='fill', dtype=np.float64):
         if self._fdata_cache is not None:
             if self._fdata_cache.dtype.type == dtype.type:
                 return self._fdata_cache
-        data = np.asanyarray(self._dataobj).astype(dtype, copy=False)
+        dataobj = self._dataobj
+        # Attempt to confine data array to dtype during scaling
+        # On overflow, may still upcast
+        if is_proxy(dataobj):
+            dataobj = dataobj.get_scaled(dtype=dtype)
+        # Always return requested data type
+        # For array proxies, will only copy on overflow
+        data = np.asanyarray(dataobj, dtype=dtype)
         if caching == 'fill':
             self._fdata_cache = data
         return data