Repository: openmrslab/suspect

Issue:
Get spatial information out of twix files
At the moment, only rda files load a transform to enable them to be co-registered to structural MRI images. It should in theory be possible to read this information out of the twix file as well.

Discussion:
coveralls: 
[![Coverage Status](https://coveralls.io/builds/10809685/badge)](https://coveralls.io/builds/10809685)

Coverage increased (+0.5%) to 75.646% when pulling **ff8d94bd9539a608b8588295309b45ac9a9ffb5d on 45_twix_transform** into **964f2460e46378c29e78d280f999128f34e829df on master**.


Patch:
diff --git a/suspect/_transforms.py b/suspect/_transforms.py
index 78bc68a..caddbfe 100644
--- a/suspect/_transforms.py
+++ b/suspect/_transforms.py
@@ -38,3 +38,33 @@ def transformation_matrix(x_vector, y_vector, translation, spacing):
         for j in range(4):
             matrix[i, j] *= spacing[j]
     return matrix
+
+
+def rotation_matrix(angle, axis):
+    """
+    Creates a 3x3 matrix which rotates `angle` radians around `axis`
+    
+    Parameters
+    ----------
+    angle : float
+        The angle in radians to rotate around the axis
+    axis : array
+        The unit vector around which to rotate
+        
+    Returns
+    -------
+    matrix : array
+    """
+    c = numpy.cos(angle)
+    s = numpy.sin(angle)
+    matrix = numpy.zeros((3, 3))
+    matrix[0, 0] = c + axis[0] ** 2 * (1 - c)
+    matrix[0, 1] = axis[0] * axis[1] * (1 - c) - axis[2] * s
+    matrix[0, 2] = axis[0] * axis[2] * (1 - c) + axis[1] * s
+    matrix[1, 0] = axis[1] * axis[0] * (1 - c) + axis[2] * s
+    matrix[1, 1] = c + axis[1] ** 2 * (1 - c)
+    matrix[1, 2] = axis[1] * axis[2] * (1 - c) - axis[0] * s
+    matrix[2, 0] = axis[2] * axis[0] * (1 - c) - axis[1] * s
+    matrix[2, 1] = axis[2] * axis[1] * (1 - c) + axis[0] * s
+    matrix[2, 2] = c + axis[2] ** 2 * (1 - c)
+    return matrix
diff --git a/suspect/_version.py b/suspect/_version.py
index d3ec452..acafb4d 100644
--- a/suspect/_version.py
+++ b/suspect/_version.py
@@ -1,1 +1,1 @@
-__version__ = "0.2.0"
+__version__ = "0.3.0a"
diff --git a/suspect/io/twix.py b/suspect/io/twix.py
index 6f5c335..9daa82b 100644
--- a/suspect/io/twix.py
+++ b/suspect/io/twix.py
@@ -1,7 +1,8 @@
-from suspect import MRSData
+from suspect import MRSData, transformation_matrix, rotation_matrix
 
 import struct
 import numpy
+#import quaternion
 import re
 
 # This file largely relies on information from Siemens regarding the structure
@@ -58,13 +59,25 @@ class TwixBuilder(object):
             "patient_id": self.header_params["patient_id"],
             "patient_birthdate": self.header_params["patient_birthdate"]
         }
-        mrs_data = MRSData(data, self.header_params["dt"], self.header_params["f0"], metadata=metadata)
+        mrs_data = MRSData(data,
+                           self.header_params["dt"],
+                           self.header_params["f0"],
+                           metadata=metadata,
+                           transform=self.header_params["transform"])
 
         return mrs_data
 
 
+def read_double(name, header_string):
+    substring = re.search(r"<ParamDouble.\"{}\">  {{ <Precision> \d+(  -?[0-9\.]+)?  }}".format(name), header_string)
+    if not substring:
+        raise KeyError(r'ParamDouble."{}" not found in header string'.format(name))
+    number_string = substring.group(1)
+    return float(number_string) if number_string else 0
+
+
 def parse_twix_header(header_string):
-    # print(header_string)
+    #print(header_string)
     # get the name of the protocol being acquired
     protocol_name_string = re.search(r"<ParamString.\"tProtocolName\">  { \".+\"  }\n", header_string).group()
     protocol_name = protocol_name_string.split("\"")[3]
@@ -100,12 +113,49 @@ def parse_twix_header(header_string):
             break
     else:
         raise KeyError("Unable to identify Dwell Time from header")
+
+    # get voxel size
+    ro_fov = read_double("VoI_RoFOV", header_string)
+    pe_fov = read_double("VoI_PeFOV", header_string)
+    slice_thickness = read_double("VoI_SliceThickness", header_string)
+
+    # get position information
+    pos_sag = read_double("VoI_Position_Sag", header_string)
+    pos_cor = read_double("VoI_Position_Cor", header_string)
+    pos_tra = read_double("VoI_Position_Tra", header_string)
+
+    # get orientation information
+    in_plane_rot = read_double("VoI_InPlaneRotAngle", header_string)
+    normal_sag = read_double("VoI_Normal_Sag", header_string)
+    normal_cor = read_double("VoI_Normal_Cor", header_string)
+    normal_tra = read_double("VoI_Normal_Tra", header_string)
+
+    # the orientation is stored in a somewhat strange way - a normal vector and
+    # a rotation angle. to get the row vector, we first use Gram-Schmidt to
+    # make [-1, 0, 0] (the default row vector) orthogonal to the normal, and
+    # then rotate that vector by the rotation angle (which we do here with a
+    # quaternion (not any more, quaternion library has issues with Travis)
+    x_vector = numpy.array([-1, 0, 0])
+    normal_vector = numpy.array([normal_sag, normal_cor, normal_tra])
+    orthogonal_x = x_vector - numpy.dot(x_vector, normal_vector) * normal_vector
+    orthonormal_x = orthogonal_x / numpy.linalg.norm(orthogonal_x)
+    #rotation_quaternion = quaternion.from_rotation_vector(in_plane_rot * normal_vector)
+    #row_vector2 = quaternion.rotate_vectors(rotation_quaternion, orthonormal_x)
+    rot_matrix = rotation_matrix(in_plane_rot, normal_vector)
+    row_vector = numpy.dot(rot_matrix, orthonormal_x)
+    column_vector = numpy.cross(row_vector, normal_vector)
+    transform = transformation_matrix(row_vector,
+                                      column_vector,
+                                      [pos_sag, pos_cor, pos_tra],
+                                      [ro_fov, pe_fov, slice_thickness])
+
     return {"protocol_name": protocol_name,
             "patient_name": patient_name,
             "patient_id": patient_id,
             "patient_birthdate": patient_birthday,
             "dt": dwell_time,
-            "f0": frequency
+            "f0": frequency,
+            "transform": transform
             }
 
 


Test patch:
diff --git a/tests/test_mrs/test_twix.py b/tests/test_mrs/test_twix.py
index 4c95f5f..2b3b374 100644
--- a/tests/test_mrs/test_twix.py
+++ b/tests/test_mrs/test_twix.py
@@ -15,6 +15,12 @@ def test_veriofile():
     assert data.np == 2048
     assert data.dt == 2.5e-4
     numpy.testing.assert_almost_equal(data.f0, 123.261716)
+    numpy.testing.assert_allclose(data.transform, numpy.array(
+        [[-20, 0, 0, 4.917676],
+         [0, 20, 0, 57.525424],
+         [0, 0, -20, 43.220339],
+         [0, 0, 0, 1]]
+    ))
 
 
 #def test_skyra():