BF+TEST - I was using the intercept when scale was 0, because I was fixing scale=0 to scale=1

Author: matthew-brett

nibabel/analyze.py

@@ -550,6 +550,11 @@ def raw_data_from_fileobj(self, fileobj):
     def data_from_fileobj(self, fileobj):
         ''' Read scaled data array from `fileobj`
 
+        Use this routine to get the scaled image data from an image file
+        `fileobj`, given a header `self`.  "Scaled" means, with any header
+        scaling factors applied to the raw data in the file.  Use
+        `raw_data_from_fileobj` to get the raw data.
+
         Parameters
         ----------
         fileobj : file-like
@@ -559,12 +564,19 @@ def data_from_fileobj(self, fileobj):
         -------
         arr : ndarray
            scaled data array
+
+        Notes
+        -----
+        We use the header to get any scale or intercept values to apply to the
+        data.  Raw Analyze files don't have scale factors or intercepts, but
+        this routine also works with formats based on Analyze, that do have
+        scaling, such as SPM analyze formats and NIfTI.
         '''
         # read unscaled data
         data = self.raw_data_from_fileobj(fileobj)
-        # get scalings from header
+        # 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 not 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.
@@ -1086,33 +1098,30 @@ def get_slope_inter(self):
 
         These are not implemented for basic Analyze
         '''
-        return 1.0, 0.0
+        return None, None
 
-    def set_slope_inter(self, slope, inter=0.0):
+    def set_slope_inter(self, slope, inter=None):
         ''' Set slope and / or intercept into header
 
         Set slope and intercept for image data, such that, if the image
         data is ``arr``, then the scaled image data will be ``(arr *
         slope) + inter``
 
-        Note that trying to set not-default values raises error for
-        Analyze header - which cannot contain slope or intercept terms.
+        In this case, for Analyze images, we can't store the slope or the
+        intercept, so this method only checks that `slope` is None or 1.0, and
+        that `inter` is None or 0.
 
         Parameters
         ----------
         slope : None or float
-           If None, implies `slope` of 1.0, `inter` of 0.0 (i.e. no
-           scaling of the image data).  If `slope` is None, we ignore
-           the passed value of `inter`
-        inter : float, optional
-           intercept
+            If float, value must be 1.0 or we raise a ``HeaderTypeError``
+        inter : None or float, optional
+            If float, value must be 0.0 or we raise a ``HeaderTypeError``
         '''
-        if slope is None:
-            slope = 1.0
-            inter = 0.0
-        if slope != 1.0 or inter:
-            raise HeaderTypeError('Cannot set slope or intercept '
-                                  'for Analyze headers')
+        if (slope is None or slope == 1.0) and (inter is None or inter == 0):
+            return
+        raise HeaderTypeError('Cannot set slope != 1 or intercept != 0 '
+                              'for Analyze headers')
 
     def scaling_from_data(self, data):
         ''' Calculate slope, intercept, min, max from data given header

nibabel/nifti1.py

@@ -17,7 +17,6 @@
 from .quaternions import fillpositive, quat2mat, mat2quat
 from . import analyze # module import
 from .spm99analyze import SpmAnalyzeHeader
-from .fileholders import copy_file_map
 
 # nifti1 flat header definition for Analyze-like first 348 bytes
 # first number in comments indicates offset in file header in bytes
@@ -752,14 +751,13 @@ def get_slope_inter(self):
         Returns
         -------
         slope : None or float
-           scaling (slope).  None if there is no valid scaling from
-           these fields
+           scaling (slope).  None if there is no valid scaling from these fields
         inter : None or float
-           offset (intercept).  Also None if there is no valid scaling, offset
+           offset (intercept). None if there is no valid scaling or if offset is
+           not finite.
 
         Examples
         --------
-        >>> fields = {'scl_slope':1,'scl_inter':0}
         >>> hdr = Nifti1Header()
         >>> hdr.get_slope_inter()
         (1.0, 0.0)
@@ -775,14 +773,14 @@ def get_slope_inter(self):
         (1.0, 1.0)
         >>> hdr['scl_inter'] = np.inf
         >>> hdr.get_slope_inter()
-        (1.0, 0.0)
+        (1.0, None)
         '''
         scale = float(self['scl_slope'])
         dc_offset = float(self['scl_inter'])
-        if not scale or not np.isfinite(scale):
+        if scale == 0 or not np.isfinite(scale):
             return None, None
         if not np.isfinite(dc_offset):
-            dc_offset = 0.0
+            dc_offset = None
         return scale, dc_offset
 
     def set_slope_inter(self, slope, inter=0.0):
@@ -795,14 +793,15 @@ def set_slope_inter(self, slope, inter=0.0):
         Parameters
         ----------
         slope : None or float
-           If None, implies `slope` of 1.0, `inter` of 0.0 (i.e. no
-           scaling of the image data).  If `slope` is not, we ignore the
-           passed value of `inter`
-        inter : float, optional
-           intercept
+           If None, implies `slope`  of 0. When the slope is set to 0 or a
+           not-finite value, ``get_slope_inter`` returns (None, None), i.e.
+           `inter` is ignored unless there is a valid value for `slope`.
+        inter : None or float, optional
+           intercept.  None implies inter value of 0.
         '''
         if slope is None:
-            slope = 1.0
+            slope = 0.0
+        if inter is None:
             inter = 0.0
         self._header_data['scl_slope'] = slope
         self._header_data['scl_inter'] = inter
@@ -1224,37 +1223,54 @@ def _get_checks(klass):
                 klass._chk_datatype,
                 klass._chk_bitpix,
                 klass._chk_pixdims,
-                klass._chk_scale_slope,
                 klass._chk_scale_inter,
                 klass._chk_qfac,
                 klass._chk_magic_offset,
                 klass._chk_qform_code,
                 klass._chk_sform_code)
 
     @staticmethod
-    def _chk_scale_slope(hdr, fix=False):
+    def _chk_scale_inter(hdr, fix=False):
         rep = Report(HeaderDataError)
         scale = hdr['scl_slope']
-        if scale and np.isfinite(scale):
+        offset = hdr['scl_inter']
+        usable_scale = np.isfinite(scale) and scale !=0
+        # Nonzero finite scale, and valid offset
+        if usable_scale and np.isfinite(offset) or (offset, scale) == (0, 0):
             return hdr, rep
-        rep.problem_level = 30
-        rep.problem_msg = '"scl_slope" is %s; should !=0 and be finite' % scale
-        if fix:
-            hdr['scl_slope'] = 1
-            rep.fix_msg = 'setting "scl_slope" to 1'
-        return hdr, rep
-
-    @staticmethod
-    def _chk_scale_inter(hdr, fix=False):
-        rep = Report(HeaderDataError)
-        scale = hdr['scl_inter']
-        if np.isfinite(scale):
+        # If scale is usable but the intercept is not finite, that's a serious
+        # problem
+        if usable_scale and not np.isfinite(offset):
+            rep.problem_level = 40
+            rep.problem_msg = ('"scl_slope" is %s; but "scl_inter" is %s; '
+                               '"scl_inter" should be finite'
+                               % (scale, offset))
+            if fix:
+                hdr['scl_inter'] = 0
+                rep.fix_msg = 'setting "scl_inter" to 0'
             return hdr, rep
-        rep.problem_level = 30
-        rep.problem_msg = '"scl_inter" is %s; should be finite' % scale
-        if fix:
-            hdr['scl_inter'] = 0
-            rep.fix_msg = 'setting "scl_inter" to 0'
+        level = 0
+        msgs = []
+        fix_msgs = []
+        # Non-finite scale is obviously an error.  We still need to check the
+        # intercept though
+        if not np.isfinite(scale):
+            level = 30
+            msgs.append('"scl_slope" is %s; should be finite' % scale)
+            if fix:
+                hdr['scl_slope'] = 0
+                fix_msgs.append('setting "scl_slope" to 0 (no scaling)')
+        # We've established scale is not usable, so inter will be ignored.  That
+        # means we can go a bit easy on bad intercepts
+        if offset != 0:
+            if level == 0: level = 20
+            msgs.append('Unused "scl_inter" is %s; should be 0' % offset)
+            if fix:
+                hdr['scl_inter'] = 0
+                fix_msgs.append('setting "scl_inter" to 0')
+        rep.problem_level = level
+        rep.problem_msg = '; '.join(msgs)
+        rep.fix_msg = '; '.join(fix_msgs)
         return hdr, rep
 
     @staticmethod

nibabel/spm2analyze.py

@@ -115,7 +115,8 @@ def get_slope_inter(self):
     def _chk_scale(klass, hdr, fix=True):
         rep = Report(HeaderDataError)
         scale, offset = hdr.get_slope_inter()
-        if not scale is None:
+        # scl_slope of 0 is valid and implies no scaling OR intercept
+        if not scale is None or hdr['scl_slope'] == 0:
             return hdr, rep
         rep.problem_level = 30
         rep.problem_msg = ('no valid scaling in scalefactor (=%s) '

nibabel/spm99analyze.py

@@ -50,12 +50,17 @@ def _empty_headerdata(self, endianness=None):
         return hdr_data
 
     def get_slope_inter(self):
-        ''' Get scalefactor and intercept '''
+        ''' Get scalefactor and intercept 
+
+        If scalefactor is 0.0 return None to indicate no scalefactor.  Intercept
+        is always None because SPM99 analyze cannot store intercepts.
+        '''
         slope = self._header_data['scl_slope']
-        inter = 0.0
-        return slope, inter
+        if slope == 0.0:
+            return None, None
+        return slope, None
 
-    def set_slope_inter(self, slope, inter=0.0):
+    def set_slope_inter(self, slope, inter=None):
         ''' Set slope and / or intercept into header
 
         Set slope and intercept for image data, such that, if the image
@@ -71,16 +76,17 @@ def set_slope_inter(self, slope, inter=0.0):
            If None, implies `slope` of 1.0, `inter` of 0.0 (i.e. no
            scaling of the image data).  If `slope` is not, we ignore the
            passed value of `inter`
-        inter : float, optional
-           intercept
+        inter : None or float, optional
+           intercept (dc offset).  If float, must be 0, because SPM99 cannot
+           store intercepts.
         '''
         if slope is None:
-            slope = 1.0
-            inter = 0.0
+            slope = 0.0
         self._header_data['scl_slope'] = slope
-        if inter:
-            raise HeaderTypeError('Cannot set non-zero intercept '
-                                  'for SPM headers')
+        if inter is None or inter == 0:
+            return
+        raise HeaderTypeError('Cannot set non-zero intercept '
+                              'for SPM headers')
 
     @classmethod
     def _get_checks(klass):
@@ -91,10 +97,10 @@ def _get_checks(klass):
     def _chk_scale(hdr, fix=False):
         rep = Report(HeaderDataError)
         scale = hdr['scl_slope']
-        if scale and np.isfinite(scale):
+        if np.isfinite(scale):
             return hdr, rep
         rep.problem_level = 30
-        rep.problem_msg = ('scale slope is %s; should !=0 and be finite'
+        rep.problem_msg = ('scale slope is %s; should be finite'
                            % scale)
         if fix:
             hdr['scl_slope'] = 1

nibabel/tests/test_analyze.py

@@ -67,20 +67,31 @@
 
 def _log_chk(hdr, level):
     # utility function to check header checking / logging
+    # If level == 0, this header should always be OK
     str_io = StringIO()
     logger = logging.getLogger('test.logger')
     handler = logging.StreamHandler(str_io)
     logger.addHandler(handler)
     str_io.truncate(0)
+    hdrc = hdr.copy()
+    if level == 0: # Should never log or raise error
+        logger.setLevel(0)
+        hdrc.check_fix(logger=logger, error_level=0)
+        assert_true(str_io.getvalue() == '')
+        logger.removeHandler(handler)
+        return hdrc, '', ()
+    # Non zero level, test above and below threshold
+    # Logging level above threshold, no log
     logger.setLevel(level+1)
     e_lev = level+1
-    hdrc = hdr.copy()
     hdrc.check_fix(logger=logger, error_level=e_lev)
-    assert(str_io.getvalue() == '')
+    assert_true(str_io.getvalue() == '')
+    # Logging level below threshold, log appears
+    logger.setLevel(level+1)
     logger.setLevel(level-1)
     hdrc = hdr.copy()
     hdrc.check_fix(logger=logger, error_level=e_lev)
-    assert(str_io.getvalue() != '')
+    assert_true(str_io.getvalue() != '')
     message = str_io.getvalue().strip()
     logger.removeHandler(handler)
     hdrc2 = hdr.copy()
@@ -91,7 +102,6 @@ def _log_chk(hdr, level):
     return hdrc, message, raiser
 
 
-
 class TestAnalyzeHeader(tb._TestBinaryHeader):
     header_class = AnalyzeHeader
     example_file = header_file
@@ -338,15 +348,19 @@ def test_scaling():
     yield assert_false(np.allclose(data_p5, rdata))
 
 
-@parametric
 def test_slope_inter():
     hdr = AnalyzeHeader()
-    yield assert_equal(hdr.get_slope_inter(), (1.0, 0.0))
-    hdr.set_slope_inter(None)
-    yield assert_equal(hdr.get_slope_inter(), (1.0, 0.0))
-    hdr.set_slope_inter(1.0)
-    yield assert_equal(hdr.get_slope_inter(), (1.0, 0.0))
-    yield assert_raises(HeaderTypeError, hdr.set_slope_inter, 1.1)
+    assert_equal(hdr.get_slope_inter(), (None, None))
+    for slinter in ((None,),
+                    (None, None),
+                    (1.0,),
+                    (1.0, None),
+                    (None, 0),
+                    (1.0, 0)):
+        hdr.set_slope_inter(*slinter)
+        assert_equal(hdr.get_slope_inter(), (None, None))
+    assert_raises(HeaderTypeError, hdr.set_slope_inter, 1.1)
+    assert_raises(HeaderTypeError, hdr.set_slope_inter, 1.0, 0.1)
 
 
 class TestAnalyzeImage(tsi.TestSpatialImage):