RF+TST - use apply scaling upcast for image load

Use safer apply_read_scaling routine to - er - apply the scaling from
the header to the read data.

Make tests to check for upcasting.  The tests revealed that the
upcasting was mostly happening already due to accidental float casting
when getting the slope and intercept for nifti1 and spm2analyze.

Expand tests of the underly int_scinter_ftype routine

Author: matthew-brett

nibabel/analyze.py

@@ -85,7 +85,8 @@
 import numpy as np
 
 from .volumeutils import (native_code, swapped_code, make_dt_codes, allopen,
-                          shape_zoom_affine, array_from_file, seek_tell)
+                          shape_zoom_affine, array_from_file, seek_tell,
+                          apply_read_scaling)
 from .arraywriters import make_array_writer, get_slope_inter, WriterError
 from .wrapstruct import WrapStruct
 from .spatialimages import (HeaderDataError, HeaderTypeError,
@@ -485,24 +486,10 @@ def data_from_fileobj(self, fileobj):
         data = self.raw_data_from_fileobj(fileobj)
         # get scalings from header.  Value of None means not present in header
         slope, inter = self.get_slope_inter()
-        if slope is None or (slope==1.0 and (inter is None or inter == 0)):
-            return data
-        # in-place multiplication and addition on integer types leads to
-        # integer output types, and disastrous integer rounding.
-        # We'd like to do inplace if we can, to save memory
-        is_flt = data.dtype.kind in 'fc'
-        if slope != 1.0:
-            if is_flt:
-                data *= slope
-            else:
-                data = data * slope
-                is_flt = True
-        if inter:
-            if is_flt:
-                data += inter
-            else:
-                data = data + inter
-        return data
+        slope = 1.0 if slope is None else slope
+        inter = 0.0 if inter is None else inter
+        # Upcast as necessary for big slopes, intercepts
+        return apply_read_scaling(data, slope, inter)
 
     def data_to_fileobj(self, data, fileobj):
         ''' Write `data` to `fileobj`, maybe modifying `self`

nibabel/nifti1.py

@@ -825,6 +825,8 @@ def get_slope_inter(self):
         >>> hdr.get_slope_inter()
         (1.0, None)
         '''
+        # Note that we are returning float (float64) scalefactors and
+        # intercepts, although they are stored as np.float32.
         scale = float(self['scl_slope'])
         dc_offset = float(self['scl_inter'])
         if scale == 0 or not np.isfinite(scale):

nibabel/tests/test_nifti1.py

@@ -62,6 +62,21 @@ def test_from_eg_file(self):
         assert_equal(hdr['magic'], asbytes('ni1'))
         assert_equal(hdr['sizeof_hdr'], 348)
 
+    def test_big_scaling(self):
+        # Test that upcasting works for huge scalefactors
+        # See tests for apply_read_scaling in test_utils
+        hdr = self.header_class()
+        hdr.set_data_shape((2,1,1))
+        hdr.set_data_dtype(np.int16)
+        sio = BytesIO()
+        dtt = np.float32
+        # This will generate a huge scalefactor
+        finf = np.finfo(dtt)
+        data = np.array([finf.min, finf.max], dtype=dtt)[:,None, None]
+        hdr.data_to_fileobj(data, sio)
+        data_back = hdr.data_from_fileobj(sio)
+        assert_true(np.allclose(data, data_back))
+
     def test_nifti_log_checks(self):
         # in addition to analyze header checks
         HC = self.header_class

nibabel/tests/test_spm99analyze.py

@@ -53,6 +53,20 @@ def test_scaling(self):
         data_back2 = hdr.data_from_fileobj(S3)
         assert_array_equal(data_back, data_back2, 4)
 
+    def test_big_scaling(self):
+        # Test that upcasting works for huge scalefactors
+        # See tests for apply_read_scaling in test_utils
+        hdr = self.header_class()
+        hdr.set_data_shape((1,1,1))
+        hdr.set_data_dtype(np.int16)
+        sio = BytesIO()
+        dtt = np.float32
+        # This will generate a huge scalefactor
+        data = np.array([np.finfo(dtt).max], dtype=dtt)[:,None, None]
+        hdr.data_to_fileobj(data, sio)
+        data_back = hdr.data_from_fileobj(sio)
+        assert_true(np.allclose(data, data_back))
+
     def test_origin_checks(self):
         HC = self.header_class
         # origin

nibabel/tests/test_utils.py

@@ -28,6 +28,8 @@
                            shape_zoom_affine,
                            rec2dict)
 
+from ..casting import flt2nmant, floor_log2
+
 from numpy.testing import (assert_array_almost_equal,
                            assert_array_equal)
 
@@ -233,16 +235,28 @@ def test_apply_scaling():
     assert_true(apply_read_scaling(arr) is arr)
     assert_true(apply_read_scaling(arr, np.float64(1.0)) is arr)
     assert_true(apply_read_scaling(arr, inter=np.float64(0)) is arr)
-    # Check upcast
+    f32, f64 = np.float32, np.float64
+    f32_arr = np.zeros((1,), dtype=f32)
+    i16_arr = np.zeros((1,), dtype=np.int16)
+    # Check float upcast (not the normal numpy scalar rule)
+    # This is the normal rule - no upcast from scalar
+    assert_equal((f32_arr * f64(1)).dtype, np.float32)
+    assert_equal((f32_arr + f64(1)).dtype, np.float32)
+    # The function does upcast though
     ret = apply_read_scaling(np.float32(0), np.float64(2))
     assert_equal(ret.dtype, np.float64)
     ret = apply_read_scaling(np.float32(0), inter=np.float64(2))
     assert_equal(ret.dtype, np.float64)
     # Check integer inf upcast
-    arr = np.zeros((1,), dtype=np.int32)
-    big = np.array([np.finfo(np.float32).max], dtype=np.float32)
-    assert_equal(apply_read_scaling(arr, big).dtype, np.float64)
-    assert_equal(apply_read_scaling(arr, 1.0, big).dtype, np.float64)
+    big = f32(np.finfo(f32).max)
+    # Normally this would not upcast
+    assert_equal((i16_arr * big).dtype, np.float32)
+    # An equivalent case is a little hard to find for the intercept
+    big_delta = np.float32(2**(floor_log2(big)-flt2nmant(np.float32)))
+    assert_equal((i16_arr * big_delta + big).dtype, np.float32)
+    # 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)

nibabel/volumeutils.py

@@ -915,6 +915,16 @@ def int_scinter_ftype(ifmt, slope=1.0, inter=0.0, default=np.float32):
     True
     >>> int_scinter_ftype(np.int8, 1e38, 0.0) == np.float64
     True
+
+    Notes
+    -----
+    It is difficult to make floats overflow with just addition because the
+    deltas are so large at the extremes of floating point.  For example::
+
+        >>> arr = np.array([np.finfo(np.float32).max], dtype=np.float32)
+        >>> res = arr + np.iinfo(np.int16).max
+        >>> arr == res
+        array([ True], dtype=bool)
     """
     floats = np.sctypes['float']
     def_ind = floats.index(default)