BF - quaternion fill threshold was not adjusted for float32 data

Author: matthew-brett

nibabel/nifti1.py

@@ -18,6 +18,9 @@
 from . import analyze # module import
 from .spm99analyze import SpmAnalyzeHeader
 
+# Needed for quaternion calculation
+FLOAT32_EPS_3 = -np.finfo(np.float32).eps * 3
+
 # nifti1 flat header definition for Analyze-like first 348 bytes
 # first number in comments indicates offset in file header in bytes
 header_dtd = [
@@ -578,7 +581,8 @@ def get_qform_quaternion(self):
         '''
         hdr = self._header_data
         bcd = [hdr['quatern_b'], hdr['quatern_c'], hdr['quatern_d']]
-        return fillpositive(bcd)
+        # Adjust threshold to fact that source data was float32
+        return fillpositive(bcd, FLOAT32_EPS_3)
 
     def get_qform(self):
         ''' Return 4x4 affine matrix from qform parameters in header '''

nibabel/quaternions.py

@@ -87,9 +87,9 @@ def fillpositive(xyz, w2_thresh=None):
     # If necessary, guess precision of input
     if w2_thresh is None:
         try: # trap errors for non-array, integer array
-            w2_thresh = -np.finfo(xyz.dtype).eps
+            w2_thresh = -np.finfo(xyz.dtype).eps * 3
         except (AttributeError, ValueError):
-            w2_thresh = -FLOAT_EPS
+            w2_thresh = -FLOAT_EPS * 3
     # Use maximum precision
     xyz = np.asarray(xyz, dtype=MAX_FLOAT)
     # Calculate w

nibabel/tests/test_nifti1.py

@@ -205,13 +205,14 @@ def test_quaternion():
     hdr['quatern_b'] = 0
     hdr['quatern_c'] = 0
     hdr['quatern_d'] = 0
-    yield assert_true, np.allclose(hdr.get_qform_quaternion(),
-                       [1.0, 0, 0, 0])
+    assert_true(np.allclose(hdr.get_qform_quaternion(), [1.0, 0, 0, 0]))
     hdr['quatern_b'] = 1
     hdr['quatern_c'] = 0
     hdr['quatern_d'] = 0
-    yield assert_true, np.allclose(hdr.get_qform_quaternion(),
-                       [0, 1, 0, 0])
+    assert_true(np.allclose(hdr.get_qform_quaternion(), [0, 1, 0, 0]))
+    # Check threshold set correctly for float32
+    hdr['quatern_b'] = 1+np.finfo(np.float32).eps
+    assert_array_almost_equal(hdr.get_qform_quaternion(), [0, 1, 0, 0])
 
 
 def test_qform():