Loader for XA DICOM

suspect/io/siemens.py

@@ -6,6 +6,8 @@
 
 from suspect import MRSData, transformation_matrix, rotation_matrix
 from ._common import complex_array_from_iter
+from .dicom import load_dicom
+from .twix import calculate_orientation
 
 CSA1 = 0
 CSA2 = 1
@@ -94,18 +96,34 @@ def read_csa_header(csa_header_bytes):
 
 
 def load_siemens_dicom(filename):
-    """Imports a file in the Siemens .IMA format.
+    """Imports a file in the Siemens .IMA format for non-XA version and .dcm for XA 
+    version.
 
     Parameters
     ----------
     filename : str
         The name of the file to import
 
+    """
+    # Start by reading in the DICOM file completely.
+    dataset = pydicom.dcmread(filename)
+    software_version = dataset[0x0018, 0x1020].value
+    if "XA" in software_version:
+        return _load_siemens_dicom_xa(dataset)
+    else:
+        return _load_siemens_dicom_nonxa(dataset)
+
+def _load_siemens_dicom_nonxa(dataset):
+    """Imports a file in the Siemens .IMA format for older/non-XA version.
+
+    Parameters
+    ----------
+    dataset : pydicom.dataset.FileDataset
+        Loaded DICOM dataset object
+
     """
     # the .IMA format is a DICOM standard, unfortunately most of the information is contained inside a private and very
     # complicated header with its own data storage format, we have to get that information out along with the data
-    # start by reading in the DICOM file completely
-    dataset = pydicom.dcmread(filename)
     # now loop through the available element in group (0029) to work out the element 
     # number of the csa header
     header_index = 0
@@ -197,6 +215,57 @@ def load_siemens_dicom(filename):
                    metadata=metadata)
 
 
+def _load_siemens_dicom_xa(dataset):
+    """Imports a file in the Siemens .IMA format for XA version.
+
+    Parameters
+    ----------
+    dataset : pydicom.dataset.FileDataset
+        Loaded DICOM dataset object
+
+    """
+    # Newer XA version uses combination of DICOM Magnetic Resonance Spectroscopy format (SOP 1.2.840.10008.5.1.4.1.1.4.2)
+    # and a new Siemens-specific DICOM tags for some metadata (not CSA header anymore)
+    mrsdata = load_dicom(dataset.filename)
+    dt = int(dataset[0x5200, 0x9230].value[0][0x0021,0x10fe][0][0x0021, 0x1042].value) * 1e-9
+    assert dt == mrsdata.dt, "Recorded dwell time is different from 1/sw."
+    te = dataset[0x5200, 0x9229][0][0x0018, 0x9114][0][0x0018, 0x9082].value
+
+    # The voxel size etc. below could be moved into load_dicom, but looking at the test data acquired with Philips
+    # (dicom/No_Name/Mrs_Dti_Qa/MRSshortTELN_401/IM-0001-0002.dcm), the length of volume localization (0018, 9126)
+    # is 1 as opposed to 3. 
+    volume_localization = dataset[0x0018, 0x9126]
+
+    # TODO: Check against TWIX, would be great if any DICOM and TWIX from same acquisitions!
+    in_plane_rot = dataset[0x5200,0x9229][0][0x00018, 0x9112][0][0x0018,0x1314].value
+    normal_vector = numpy.array(volume_localization[0][0x0018, 0x9105].value)
+    pos_vector = numpy.array(volume_localization[0][0x0018, 0x9106].value)
+    voi_size = [volume_localization[2][0x0018, 0x9104].value,
+                volume_localization[1][0x0018, 0x9104].value,
+                volume_localization[0][0x0018, 0x9104].value]
+
+    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)
+    #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_vector,
+                                      voi_size)
+    return MRSData(numpy.array(mrsdata),
+                   dt,
+                   mrsdata.f0,
+                   te=te,
+                   tr=mrsdata.tr,
+                   transform=transform)
+
 # def anonymize_siemens_dicom(filename, anonymized_filename):
 # TODO: anonymize dicom
 #     """