Loader for XA DICOM (#178)

* Loader for XA DICOM

* Comment: In-plane rot angle not found

* Subtract VOI Position to Table Position to match TWIX

* Extract TR and added phantom XA60

suspect/io/siemens.py

@@ -1,11 +1,14 @@
 import pydicom.tag
 import pydicom
 import numpy
+import re
 import struct
 import warnings
 
 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 +97,37 @@ 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.
+    with open(filename, "rb") as f:
+        dataset = pydicom.dcmread(f)
+        software_version = dataset[0x0018, 0x1020].value
+        if "XA" in software_version:
+            f.seek(0)
+            other_ds = pydicom.filereader.read_dataset(f, True, True)
+            return _load_siemens_dicom_xa(dataset, other_ds)
+        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 +219,75 @@ def load_siemens_dicom(filename):
                    metadata=metadata)
 
 
+def _load_siemens_dicom_xa(dataset, other_ds):
+    """Imports a file in the Siemens .IMA format for XA version.
+
+    Parameters
+    ----------
+    dataset : pydicom.dataset.FileDataset
+        Loaded DICOM dataset object
+    other_ds : pydicom.dataset.FileDataset
+        Other dataset that contains meas headers 
+
+    """
+    # 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
+    tr = float(dataset[0x5200, 0x9229][0][0x0018, 0x9112][0][0x0018, 0x0080].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]
+
+    # Extract in-plane rotation angle from meas header as it is not parsed into DICOM headers
+    # Per testing, when in-plane rotation is 0, sSpecPara.sVoI.dInPlaneRot is simply missing!
+    rgx = r"sSpecPara\.sVoI\.dInPlaneRot\s*=\s*(-?[[0-9]*[.]?[0-9]*]{0,})\s*"
+    in_plane_rot_matches = re.findall(rgx, other_ds[0x0004, 0x00c0].value.decode("latin-1"))
+    if in_plane_rot_matches:
+        in_plane_rot = float(in_plane_rot_matches[0])
+    else:
+        in_plane_rot = 0
+    normal_vector = numpy.array(volume_localization[0][0x0018, 0x9105].value)
+
+    # In the older DICOM, VOI position info is in CSA header, so the position is the same as TWIX. 
+    # In XA DICOM, it's not the same as TWIX. Per testing we can 'reverse' the Position by subtracting 
+    # with TablePosition.
+    # This, of course, assumes the TWIX is the correct one to reconstruct the voxel.
+    pos_vector = numpy.array(volume_localization[0][0x0018, 0x9106].value)
+    table_position = numpy.array(
+        dataset[0x5200, 0x9230].value[0][0x0021,0x10fe][0][0x0021, 0x1059].value
+    )
+    pos_vector = pos_vector - table_position
+
+    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=tr,
+                   transform=transform)
+
 # def anonymize_siemens_dicom(filename, anonymized_filename):
 # TODO: anonymize dicom
 #     """

tests/test_mrs/test_ima.py → tests/test_mrs/test_siemens.py

@@ -12,6 +12,13 @@ def test_svs_30():
     assert data.tr == 2000
 
 
+def test_svs_xa60():
+    # Not .ima, but might be best to leave it here for Siemens spectro DICOM
+    data = suspect.io.siemens.load_siemens_dicom("tests/test_data/siemens/SVS_XA60.dcm")
+    assert data.shape == (1024,)
+    assert data.te == 30
+    assert data.tr == 2000
+
 # def test_svs_30():
 #     data = suspect.io.siemens.load_siemens_dicom("suspect/tests/test_data/siemens_svs_30.IMA")
 #     assert data.shape == (1024,)