BF: Supports reading TRK (version 1)

Author: MarcCote

nibabel/streamlines/tests/test_trk.py

@@ -156,6 +156,20 @@ def test_load_file_with_wrong_information(self):
         new_trk_file = trk_file[:254] + noise + trk_file[254+4:]
         assert_raises(HeaderError, TrkFile.load, BytesIO(new_trk_file))
 
+    def test_load_trk_version_1(self):
+        trk_file = open(DATA['simple_trk_fname'], 'rb').read()
+
+        # Simulate a TRK (version 1).
+        version = np.array(1, dtype=np.int32).tostring()
+        new_trk_file = trk_file[:992] + version + trk_file[992+4:]
+        with clear_and_catch_warnings(record=True, modules=[trk_module]) as w:
+            trk = TrkFile.load(BytesIO(new_trk_file))
+            assert_equal(len(w), 1)
+            assert_true(issubclass(w[0].category, HeaderWarning))
+            assert_true("identity" in str(w[0].message))
+            assert_array_equal(trk.affine, np.eye(4))
+            assert_array_equal(trk.header['version'], 1)
+
     def test_load_complex_file_in_big_endian(self):
         trk_file = open(DATA['complex_trk_big_endian_fname'], 'rb').read()
         # We use hdr_size as an indicator of little vs big endian.

nibabel/streamlines/trk.py

@@ -26,36 +26,11 @@
 MAX_NB_NAMED_SCALARS_PER_POINT = 10
 MAX_NB_NAMED_PROPERTIES_PER_STREAMLINE = 10
 
-# See http://docs.scipy.org/doc/numpy/reference/arrays.dtypes.html
-header_1_dtd = [(Field.MAGIC_NUMBER, 'S6'),
-                (Field.DIMENSIONS, 'h', 3),
-                (Field.VOXEL_SIZES, 'f4', 3),
-                (Field.ORIGIN, 'f4', 3),
-                (Field.NB_SCALARS_PER_POINT, 'h'),
-                ('scalar_name', 'S20', MAX_NB_NAMED_SCALARS_PER_POINT),
-                (Field.NB_PROPERTIES_PER_STREAMLINE, 'h'),
-                ('property_name', 'S20',
-                 MAX_NB_NAMED_PROPERTIES_PER_STREAMLINE),
-                ('reserved', 'S508'),
-                (Field.VOXEL_ORDER, 'S4'),
-                ('pad2', 'S4'),
-                ('image_orientation_patient', 'f4', 6),
-                ('pad1', 'S2'),
-                ('invert_x', 'S1'),
-                ('invert_y', 'S1'),
-                ('invert_z', 'S1'),
-                ('swap_xy', 'S1'),
-                ('swap_yz', 'S1'),
-                ('swap_zx', 'S1'),
-                (Field.NB_STREAMLINES, 'i4'),
-                ('version', 'i4'),
-                ('hdr_size', 'i4'),
-                ]
-
-# Version 2 adds a 4x4 matrix giving the affine transformtation going
+# Version 2 adds a 4x4 matrix giving the affine transformation going
 # from voxel coordinates in the referenced 3D voxel matrix, to xyz
 # coordinates (axes L->R, P->A, I->S). If (0 based) value [3, 3] from
 # this matrix is 0, this means the matrix is not recorded.
+# See http://docs.scipy.org/doc/numpy/reference/arrays.dtypes.html
 header_2_dtd = [(Field.MAGIC_NUMBER, 'S6'),
                 (Field.DIMENSIONS, 'h', 3),
                 (Field.VOXEL_SIZES, 'f4', 3),
@@ -83,7 +58,6 @@
                 ]
 
 # Full header numpy dtypes
-header_1_dtype = np.dtype(header_1_dtd)
 header_2_dtype = np.dtype(header_2_dtd)
 
 
@@ -584,8 +558,8 @@ def _read_header(fileobj):
                                                  TrkFile.HEADER_SIZE))
 
             if header_rec['version'] == 1:
-                header_rec = np.fromstring(string=header_str,
-                                           dtype=header_1_dtype)
+                # There is no 4x4 matrix for voxel to RAS transformation.
+                header_rec[Field.VOXEL_TO_RASMM] = np.zeros((4, 4))
             elif header_rec['version'] == 2:
                 pass  # Nothing more to do.
             else:
@@ -724,12 +698,17 @@ def __str__(self):
                 hdr_field = getattr(Field, attr)
                 if hdr_field in vars:
                     vars[attr] = vars[hdr_field]
-        vars['scalar_names'] = '\n  '.join([asstr(s)
-                                            for s in vars['scalar_name']
-                                            if len(s) > 0])
-        vars['property_names'] = "\n  ".join([asstr(s)
-                                              for s in vars['property_name']
-                                              if len(s) > 0])
+
+        nb_scalars = self.header[Field.NB_SCALARS_PER_POINT]
+        scalar_names = [asstr(s)
+                        for s in vars['scalar_name'][:nb_scalars]
+                        if len(s) > 0]
+        vars['scalar_names'] = '\n  '.join(scalar_names)
+        nb_properties = self.header[Field.NB_PROPERTIES_PER_STREAMLINE]
+        property_names = [asstr(s)
+                          for s in vars['property_name'][:nb_properties]
+                          if len(s) > 0]
+        vars['property_names'] = "\n  ".join(property_names)
         # Make all byte strings into strings
         # Fixes recursion error on Python 3.3
         vars = dict((k, asstr(v) if hasattr(v, 'decode') else v)
@@ -739,11 +718,11 @@ def __str__(self):
 v.{version}
 dim: {DIMENSIONS}
 voxel_sizes: {VOXEL_SIZES}
-orgin: {ORIGIN}
+origin: {ORIGIN}
 nb_scalars: {NB_SCALARS_PER_POINT}
-scalar_name:\n  {scalar_names}
+scalar_names:\n  {scalar_names}
 nb_properties: {NB_PROPERTIES_PER_STREAMLINE}
-property_name:\n  {property_names}
+property_names:\n  {property_names}
 vox_to_world:\n{VOXEL_TO_RASMM}
 voxel_order: {VOXEL_ORDER}
 image_orientation_patient: {image_orientation_patient}