ENH: default read scaling from slope dtype

The apply_read_scaling routine was previously ignoring the dtype of
slope and inter when choosing the floating point type to scale with.
This change makes sure that the floating point type chosen for scaling
has precision at least as high as the input slope dtype.

Author: matthew-brett

nibabel/tests/test_utils.py

@@ -312,9 +312,16 @@ def test_apply_scaling():
     # Upcasting does occur with this routine
     assert_equal(apply_read_scaling(i16_arr, big).dtype, np.float64)
     assert_equal(apply_read_scaling(i16_arr, big_delta, big).dtype, np.float64)
-    assert_equal(apply_read_scaling(np.int8(0), -1.0, 0.0).dtype, np.float32)
-    assert_equal(apply_read_scaling(np.int8(0), 1e38, 0.0).dtype, np.float64)
-    assert_equal(apply_read_scaling(np.int8(0), -1e38, 0.0).dtype, np.float64)
+    # If float32 passed, no overflow, float32 returned
+    assert_equal(apply_read_scaling(np.int8(0), f32(-1.0), f32(0.0)).dtype,
+                 np.float32)
+    # float64 passed, float64 returned
+    assert_equal(apply_read_scaling(np.int8(0), -1.0, 0.0).dtype, np.float64)
+    # float32 passed, overflow, float64 returned
+    assert_equal(apply_read_scaling(np.int8(0), f32(1e38), f32(0.0)).dtype,
+                 np.float64)
+    assert_equal(apply_read_scaling(np.int8(0), f32(-1e38), f32(0.0)).dtype,
+                 np.float64)
 
 
 def test_int_scinter():

nibabel/volumeutils.py

@@ -747,7 +747,9 @@ def apply_read_scaling(arr, slope = 1.0, inter = 0.0):
     # Force float / float upcasting by promoting to arrays
     arr, slope, inter = [np.atleast_1d(v) for v in arr, slope, inter]
     if arr.dtype.kind in 'iu':
-        ftype = int_scinter_ftype(arr.dtype, slope, inter)
+        # Find floating point type for which scaling does not overflow, starting
+        # at given type
+        ftype = int_scinter_ftype(arr.dtype, slope, inter, slope.dtype.type)
         slope = slope.astype(ftype)
         inter = inter.astype(ftype)
     if slope != 1.0: