ENH: Add handling and tests for CSI data stored in Siemens enhanced DICOM

CHANGELOG.md

@@ -3,6 +3,7 @@ This document contains the Spec2nii release history in reverse chronological ord
 0.7.1 (WIP)
 --------------------------------
 - The --anon flag can be passed with any call to anonymise after writing files.
+- The Siemens enhanced dicom filetype pathway now handles CSI data.
 
 0.7.0 (Saturday 5th August 2023)
 --------------------------------

spec2nii/Siemens/dicomfunctions.py

@@ -241,6 +241,25 @@ def process_siemens_svs(img, verbose):
         return process_siemens_svs_xa(img, verbose)
 
 
+def _get_enhanced_dcm_img_orientation_pat(img):
+    """Extract the image orientation patient field from an enhanced dicom scan
+
+    :param img: DICOM image object
+    :type img: nibabel.nicom.dicomwrappers.SiemensWrapper
+    :return: 3x2 Numpy array containing the orientation info
+    :rtype: np.ndarray
+    """
+    dcm_hdrs = img.dcm_data.PerFrameFunctionalGroupsSequence[0]
+    dcm_hdrs1 = img.dcm_data.SharedFunctionalGroupsSequence[0]
+    # 1) Extract dicom parameters
+    try:
+        return np.array(dcm_hdrs.PlaneOrientationSequence[0].ImageOrientationPatient).reshape(2, 3)
+    except AttributeError:
+        # For VB19 data formatted as MRSpectroscopyStorage this seems to be different to
+        # the XA data I've seen. But this could also be a sequence specific thing!
+        return np.array(dcm_hdrs1.PlaneOrientationSequence[0].ImageOrientationPatient).reshape(2, 3)
+
+
 def process_siemens_svs_xa(img, verbose):
     """Process Siemens DICOM SVS data acquired on a NumarisX system."""
     specData = np.frombuffer(img.dcm_data[('5600', '0020')].value, dtype=np.single)
@@ -250,13 +269,7 @@ def process_siemens_svs_xa(img, verbose):
     dcm_hdrs = img.dcm_data.PerFrameFunctionalGroupsSequence[0]
     dcm_hdrs1 = img.dcm_data.SharedFunctionalGroupsSequence[0]
     # 1) Extract dicom parameters
-    try:
-        imageOrientationPatient = np.array(dcm_hdrs.PlaneOrientationSequence[0].ImageOrientationPatient).reshape(2, 3)
-    except AttributeError:
-        # For VB19 data formatted as MRSpectroscopyStorage this seems to be different to
-        # the XA data I've seen. But this could also be a sequence specific thing!
-        imageOrientationPatient = np.array(dcm_hdrs1.PlaneOrientationSequence[0].ImageOrientationPatient).reshape(2, 3)
-    # VoiPosition - this does not have the FOV shift that imagePositionPatient has
+    imageOrientationPatient = _get_enhanced_dcm_img_orientation_pat(img)
     imagePositionPatient = dcm_hdrs.PlanePositionSequence[0].ImagePositionPatient
 
     # The first two elements could easily be the reverse of this.
@@ -309,9 +322,55 @@ def process_siemens_csi(img, verbose):
 
 
 def process_siemens_csi_xa(img, verbose):
-    """Process Siemens DICOM CSI data acquired on a NumarisX system."""
-    # NOTE: WTC expect to need to reverse the phase convention as above for svs.
-    raise NotImplementedError('Method process_siemens_csi_xa not implemented, example data needed!')
+    """Process Siemens DICOM CSI data acquired on a NumarisX system or Numaris4 Vx using enhanced DCM export"""
+
+    warnings.warn(
+        'The orientation of CSI stored in enhanced DICOM remains weakly tested.'
+        'Please provide more test data.')
+
+    specData = np.frombuffer(img.dcm_data[('5600', '0020')].value, dtype=np.single)
+    # It appears the phase convention is reversed in the NumarisX systems.
+    specDataCmplx = specData[0::2] - 1j * specData[1::2]
+
+    # The data in img.dcm_data.SharedFunctionalGroupsSequence[0].MRSpectroscopyFOVGeometrySequence[0]
+    # is the uninterpolated (lower) resolution that doesn't match the data size
+    rows = img.dcm_data.Rows
+    cols = img.dcm_data.Columns
+    slices = int(img.dcm_data.NumberOfFrames)
+    spectral_points = img.dcm_data.DataPointColumns
+
+    specDataCmplx = specDataCmplx.reshape((slices, rows, cols, spectral_points))
+    specDataCmplx = np.moveaxis(specDataCmplx, (0, 1, 2), (2, 1, 0))
+
+    dcm_hdrs = img.dcm_data.PerFrameFunctionalGroupsSequence[0]
+    dcm_hdrs1 = img.dcm_data.SharedFunctionalGroupsSequence[0]
+
+    # 1) Extract dicom parameters
+    imageOrientationPatient = _get_enhanced_dcm_img_orientation_pat(img)
+    # VoiPosition - this does not have the FOV shift that imagePositionPatient has
+    imagePositionPatient = dcm_hdrs.PlanePositionSequence[0].ImagePositionPatient
+
+    # The first two elements could easily be the reverse of this.
+    xyzMM = np.array([float(dcm_hdrs1.PixelMeasuresSequence[0].PixelSpacing[0]),
+                      float(dcm_hdrs1.PixelMeasuresSequence[0].PixelSpacing[1]),
+                      float(dcm_hdrs1.PixelMeasuresSequence[0].SliceThickness)])
+
+    # Note that half_shift = False. This is the opposite as for classic/original DCM.
+    # THIS NEEDS TESTING AND DOCUMENTING. For an explanation of original shift see spec2nii/notes/siemens.md
+    currNiftiOrientation = dcm_to_nifti_orientation(imageOrientationPatient,
+                                                    imagePositionPatient,
+                                                    xyzMM,
+                                                    specDataCmplx.shape[:3],
+                                                    half_shift=False,
+                                                    verbose=verbose)
+
+    dwelltime = 1 / img.dcm_data.SpectralWidth
+    meta = extractDicomMetadata_xa(img)
+
+    # Look for a VOI in the data
+    meta = _detect_and_fill_voi_enh(img, meta)
+
+    return specDataCmplx, currNiftiOrientation, dwelltime, meta
 
 
 def process_siemens_csi_vx(img, verbose):
@@ -388,6 +447,39 @@ def _detect_and_fill_voi(img, meta):
     return meta
 
 
+def _detect_and_fill_voi_enh(img, meta):
+    """Look for a VOI volume in the headers of a CSI scan stored in enhanced dicom.
+    If present populate the VOI header field.
+
+    :param img: Nibable DICOM image object
+    :type img: SiemensWrapper
+    :param meta: Existing NIfTI-MRS meta object
+    :type meta: hdr_ext
+    :return: Modified meta hdr_ext object
+    :rtype: hdr_ext
+    """
+
+    imageOrientationPatient = _get_enhanced_dcm_img_orientation_pat(img)
+
+    # VoiPosition - this does not have the FOV shift that imagePositionPatient has
+    vol_loc_seq = img.dcm_data.VolumeLocalizationSequence
+    imagePositionPatient = vol_loc_seq[0].MidSlabPosition
+
+    # This order needs validating!!
+    xyzMM = np.array([vol_loc_seq[1].SlabThickness,
+                      vol_loc_seq[2].SlabThickness,
+                      vol_loc_seq[0].SlabThickness])
+
+    voiNiftiOrientation = dcm_to_nifti_orientation(imageOrientationPatient,
+                                                   imagePositionPatient,
+                                                   xyzMM,
+                                                   (1, 1, 1))
+
+    meta.set_standard_def('VOI', voiNiftiOrientation.Q44.tolist())
+
+    return meta
+
+
 def extractDicomMetadata_xa(dcmdata):
     """ Extract information from the nibabel DICOM object to insert into the json header ext.
     For Numarix X systems
@@ -557,6 +649,9 @@ def set_standard_def(nifti_mrs_key, location, key, cast=None):
     # 'InstitutionAddress'
     set_standard_def('InstitutionAddress', dcmdata.dcm_data, 'InstitutionAddress')
     # 'TxCoil'
+    obj.set_standard_def(
+        'TxCoil',
+        dcmdata.csa_header['tags']['TransmittingCoil']['items'][0])
     # 'RxCoil'
     if len(dcmdata.csa_header['tags']['ReceivingCoil']['items']) > 0:
         obj.set_standard_def('RxCoil',

tests/test_siemens_dicom.py

@@ -26,6 +26,12 @@
 tim_trio_non_image = siemens_path / 'VBData' / 'tim_trio_sop_class_uid' / 'timtrio_mr_vb19a.dcm'
 tim_trio_mrspecstorage = siemens_path / 'VBData' / 'tim_trio_sop_class_uid' / 'timtrio_mrdcm_vb13a_2_0.dcm'
 
+enhanced_csi_1 = siemens_path / 'enhanced_dcm_csi' / 'rk_enhanced' / '1.dcm'
+enhanced_csi_2 = siemens_path / 'enhanced_dcm_csi'\
+    / 'sm_enhanced' / 'Phantom.Psychiatrie^Studien^aslac.Sr27.I1.20231006084036808.dcm'
+enhanced_csi_2_comp = siemens_path / 'enhanced_dcm_csi'\
+    / 'sm_classic' / 'PHANTOM.MR.PSYCHIATRIE_STUDIEN_ASLAC.0027.0001.2023.10.06.08.40.21.250955.407258759.IMA'
+
 
 def test_VB_svs(tmp_path):
 
@@ -161,3 +167,46 @@ def test_sop_class_vb(tmp_path):
 
     assert img_2.shape == (1, 1, 1, 2048)
     assert np.iscomplexobj(img_2.dataobj)
+
+
+def test_enhanced_csi(tmp_path):
+
+    # First test look for right shape (16, 16, 1, 1024)
+    subprocess.run(['spec2nii', 'dicom',
+                    '-f', 'csi_1',
+                    '-o', tmp_path,
+                    '-j', enhanced_csi_1])
+
+    img_1 = read_nifti_mrs(tmp_path / 'csi_1.nii.gz')
+    assert img_1.shape == (16, 16, 1, 1024)
+    assert np.iscomplexobj(img_1.dataobj)
+
+    # Second test: compare with classic dicom export
+    subprocess.run(['spec2nii', 'dicom',
+                    '-f', 'csi_2',
+                    '-o', tmp_path,
+                    '-j', enhanced_csi_2])
+    subprocess.run(['spec2nii', 'dicom',
+                    '-f', 'csi_2_comp',
+                    '-o', tmp_path,
+                    '-j', enhanced_csi_2_comp])
+
+    img_2 = read_nifti_mrs(tmp_path / 'csi_2.nii.gz')
+    img_2_comp = read_nifti_mrs(tmp_path / 'csi_2_comp.nii.gz')
+
+    hdr_ext_codes = img_2.header.extensions.get_codes()
+    hdr_ext = json.loads(img_2.header.extensions[hdr_ext_codes.index(44)].get_content())
+    hdr_ext_codes = img_2_comp.header.extensions.get_codes()
+    hdr_ext_comp = json.loads(img_2_comp.header.extensions[hdr_ext_codes.index(44)].get_content())
+
+    assert img_2.shape == img_2_comp.shape
+    assert np.allclose(img_2.dataobj[:], img_2_comp.dataobj[:])
+    for key in hdr_ext:
+        print(key)
+        # Skip keys with known differences.
+        if key in ('RxCoil', 'SequenceName', 'PatientID', 'ConversionTime', 'OriginalFile'):
+            continue
+        if isinstance(hdr_ext[key], (float, complex)):
+            assert np.isclose(hdr_ext[key], hdr_ext_comp[key])
+        else:
+            assert hdr_ext[key] == hdr_ext_comp[key]