Merge pull request #138 from moloney/csa_slice_norm

BF: Only use the cross product (possibly negated) for the slice normal.

The information for the slice normal can come from the crossproduct of the within-slice axes, or the private Siemens CSA header.  These changes check that these are more or less the same apart from a sign flip, and then flips the cross-product vector if there is a sign flip in the CSA header slice normal.

Author: matthew-brett

nibabel/nicom/dicomwrappers.py

@@ -407,14 +407,28 @@ def __init__(self, dcm_data=None, csa_header=None):
 
     @one_time
     def slice_normal(self):
-        slice_normal = csar.get_slice_normal(self.csa_header)
-        if not slice_normal is None:
-            return np.array(slice_normal)
-        iop = self.image_orient_patient
-        if iop is None:
+        #The std_slice_normal comes from the cross product of the directions 
+        #in the ImageOrientationPatient
+        std_slice_normal = super(SiemensWrapper, self).slice_normal
+        csa_slice_normal = csar.get_slice_normal(self.csa_header)
+        if std_slice_normal is None and csa_slice_normal is None:
             return None
-        return np.cross(*iop.T[:])
-
+        elif std_slice_normal is None:
+            return np.array(csa_slice_normal)
+        elif csa_slice_normal is None:
+            return std_slice_normal
+        else:
+            #Make sure the two normals are very close to parallel unit vectors
+            dot_prod = np.dot(csa_slice_normal, std_slice_normal)
+            assert np.allclose(np.fabs(dot_prod), 1.0, atol=1e-5)
+            #Use the slice normal computed with the cross product as it will 
+            #always be the most orthogonal, but take the sign from the CSA
+            #slice normal
+            if dot_prod < 0:
+                return -std_slice_normal
+            else:
+                return std_slice_normal
+    
     @one_time
     def series_signature(self):
         ''' Add ICE dims from CSA header to signature '''

nibabel/nicom/tests/data/csa_slice_norm.dcm

@@ -27,7 +27,8 @@ def test_read_dwi():
 
 @dicom_test
 def test_read_dwis():
-    data, aff, bs, gs = didr.read_mosaic_dwi_dir(IO_DATA_PATH, '*.dcm.gz')
+    data, aff, bs, gs = didr.read_mosaic_dwi_dir(IO_DATA_PATH, 
+                                                 'siemens_dwi_*.dcm.gz')
     assert_equal(data.ndim, 4)
     assert_array_almost_equal(aff, EXPECTED_AFFINE)
     assert_array_almost_equal(bs, (0, EXPECTED_PARAMS[0]))

nibabel/nicom/tests/test_dicomreaders.py

@@ -30,6 +30,7 @@
 else:
     DATA = None
 DATA_FILE_B0 = pjoin(IO_DATA_PATH, 'siemens_dwi_0.dcm.gz')
+DATA_FILE_SLC_NORM = pjoin(IO_DATA_PATH, 'csa_slice_norm.dcm')
 
 # This affine from our converted image was shown to match our image
 # spatially with an image from SPM DICOM conversion. We checked the
@@ -48,7 +49,6 @@
                                   0.99997450,
                                   -0.005023611)]
 
-
 @dicom_test
 def test_wrappers():
     # test direct wrapper calls
@@ -143,3 +143,31 @@ def test_slice_indicator():
     assert_equal(z, dw_1000.slice_indicator)
     dw_empty = didw.Wrapper()
     assert_true(dw_empty.slice_indicator is None)
+
+
+@dicom_test
+def test_orthogonal():
+    #Test that the slice normal is sufficiently orthogonal
+    dw = didw.wrapper_from_file(DATA_FILE_SLC_NORM)
+    R = dw.rotation_matrix
+    assert  np.allclose(np.eye(3),
+                        np.dot(R, R.T),
+                        atol=1e-6)
+                    
+@dicom_test
+def test_use_csa_sign():
+    #Test that we get the same slice normal, even after swapping the iop 
+    #directions
+    dw = didw.wrapper_from_file(DATA_FILE_SLC_NORM)
+    iop = dw.image_orient_patient
+    dw.image_orient_patient = np.c_[iop[:,1], iop[:,0]]
+    dw2 = didw.wrapper_from_file(DATA_FILE_SLC_NORM)
+    assert np.allclose(dw.slice_normal, dw2.slice_normal)
+
+@dicom_test
+def test_assert_parallel():
+    #Test that we get an AssertionError if the cross product and the CSA 
+    #slice normal are not parallel
+    dw = didw.wrapper_from_file(DATA_FILE_SLC_NORM)
+    dw.image_orient_patient = np.c_[[1., 0., 0.], [0., 1., 0.]]
+    assert_raises(AssertionError, dw.__getattribute__, 'slice_normal')