Repository: openmrslab/suspect

Issue:
LCModel IO - "VOLUME = NaN"
For 20% Siemens TWIX data I was unable to run LCModel analysis with an error 1 FATAL ZEROVX 4. In CONTROL file there was "VOLUME = NaN" instead of a some numerical value.
B.

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

Coverage decreased (-0.3%) to 74.189% when pulling **9f70dd6e6813049ae2429317e14726e47bc740af on 145_transform_issue** into **97d7a167cb6f139d5311cb194759285277ef5493 on master**.


Patch:
diff --git a/docs/changelog.rst b/docs/changelog.rst
index f4113ab..9e4d883 100644
--- a/docs/changelog.rst
+++ b/docs/changelog.rst
@@ -4,6 +4,7 @@
 Changelog
 =========
 
+* :bug:`145` fixed issue with calculating transform for sagittally planned twix voxels
 * :feature:`146` substantial overhaul of the singlet fitting module
 * :release:`0.4.2 <20/08/20>`
 * :bug:`143` changed .IMA loading to handle cases where data shape is calculated incorrectly
diff --git a/suspect/io/ge.py b/suspect/io/ge.py
index 7b597f9..cbc3ae8 100644
--- a/suspect/io/ge.py
+++ b/suspect/io/ge.py
@@ -186,9 +186,9 @@ def prepare_pfile_svs(pfile):
     mrs_data = MRSData(raw_data,
                        **header_params)
 
-    print(ref_frames)
-    print(metadata)
-    print(mrs_data.shape)
+    #print(ref_frames)
+    #print(metadata)
+    #print(mrs_data.shape)
 
     wref = mrs_data[:, :ref_frames].reshape(-1, metadata["channels"], metadata["acquiredXRes"]).squeeze()
     data = mrs_data[:, ref_frames:].squeeze()
diff --git a/suspect/io/twix.py b/suspect/io/twix.py
index a4bf8ef..e8f8879 100644
--- a/suspect/io/twix.py
+++ b/suspect/io/twix.py
@@ -72,6 +72,24 @@ class TwixBuilder(object):
         return mrs_data
 
 
+def calculate_orientation(normal):
+    """
+
+    Parameters
+    ----------
+    normal : the normal vector to the image plane
+
+    Returns
+    -------
+
+    """
+    if normal[2] >= normal[1] - 1e-6 and normal[2] >= normal[0] - 1e-6:
+        return "TRA"
+    elif normal[1] >= normal[0] - 1e-6:
+        return "COR"
+    return "SAG"
+
+
 def read_double(name, header_string):
     substring = re.search(r"<ParamDouble.\"{}\">  {{ <Precision> \d+(  -?[0-9\.]+)?  }}".format(name), header_string)
     if not substring:
@@ -154,8 +172,11 @@ def parse_twix_header(header_string):
     # 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])
+    if calculate_orientation(normal_vector) == "SAG":
+        x_vector = numpy.array([0, 0, 1])
+    else:
+        x_vector = numpy.array([-1, 0, 0])
     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)
@@ -552,3 +573,40 @@ def anonymize_twix(filename, anonymized_filename):
                 anonymize_twix_vd(fin, fout)
             else:
                 anonymize_twix_vb(fin, fout)
+
+
+def get_header(filename):
+    with open(filename, 'rb') as fin:
+
+        # we can tell the type of file from the first two uints in the header
+        first_uint, second_uint = struct.unpack("II", fin.read(8))
+
+        # reset the file pointer before giving to specific function
+        fin.seek(0)
+
+
+        if first_uint == 0 and second_uint <= 64:
+            print("twix file is VD/VE")
+            headers = []
+            for i in range(second_uint):
+                fin.seek(8 + 152 * i)
+
+                meas_id, file_id, offset, length, patient_name, protocol_name = struct.unpack("IIQQ64s64s",
+                                                                                              fin.read(152))
+                fin.seek(offset)
+                # read the header and anonymize it
+                header_size = struct.unpack("I", fin.read(4))[0]
+                header = fin.read(header_size - 4)
+                header_string = header[:-24].decode('latin-1')
+                headers.append(header_string)
+            return headers
+
+        else:
+            print("twix file is VB")
+            # read the rest of the header minus the four bytes we already read
+            fin.seek(4)
+            header = fin.read(first_uint - 4)
+            # for some reason the last 24 bytes of the header contain some stuff that
+            # is not a string, I don't know what it is
+            header_string = header[:-24].decode('latin-1')
+            return header_string


Test patch:
diff --git a/tests/test_mrs/test_twix.py b/tests/test_mrs/test_twix.py
index 0b6b9ef..7f218ee 100644
--- a/tests/test_mrs/test_twix.py
+++ b/tests/test_mrs/test_twix.py
@@ -37,3 +37,24 @@ def test_veriofile():
 #    assert data.metadata["patient_name"] == "x" * 13
 #    assert data.metadata["patient_id"] == "x" * 6
 #    assert data.metadata["patient_birthdate"] == "19700101"
+
+def test_calculate_orientation():
+    assert suspect.io.twix.calculate_orientation([1, 1, 1]) == "TRA"
+    assert suspect.io.twix.calculate_orientation([0.5, 0.5, 1]) == "TRA"
+    assert suspect.io.twix.calculate_orientation([1, 1, 0.5]) == "COR"
+    assert suspect.io.twix.calculate_orientation([0.5, 1, 0.5]) == "COR"
+    assert suspect.io.twix.calculate_orientation([1, 0.5, 1]) == "TRA"
+    assert suspect.io.twix.calculate_orientation([1, 0.5, 0.5]) == "SAG"
+    assert suspect.io.twix.calculate_orientation([0.5, 1, 1]) == "TRA"
+    assert suspect.io.twix.calculate_orientation([1, 0.5, 0.5]) == "SAG"
+    assert suspect.io.twix.calculate_orientation([0.5, 0, 1]) == "TRA"
+    assert suspect.io.twix.calculate_orientation([0, 1, 0.5]) == "COR"
+    assert suspect.io.twix.calculate_orientation([0, 0.5, 1]) == "TRA"
+    assert suspect.io.twix.calculate_orientation([0.5, 1, 0]) == "COR"
+    assert suspect.io.twix.calculate_orientation([1, 0, 0]) == "SAG"
+    assert suspect.io.twix.calculate_orientation([0, 1, 0.5]) == "COR"
+    assert suspect.io.twix.calculate_orientation([0.5, 0.5, 1]) == "TRA"
+    assert suspect.io.twix.calculate_orientation([1, 0.5, 0]) == "SAG"
+    assert suspect.io.twix.calculate_orientation([0, 1, 0]) == "COR"
+    assert suspect.io.twix.calculate_orientation([1, 0, 0.5]) == "SAG"
+    assert suspect.io.twix.calculate_orientation([0, 0, 1]) == "TRA"