BF: Tests on XA50 twix and rda, fixes for mgs_svs_ed

CHANGELOG.md

@@ -1,5 +1,12 @@
 This document contains the Spec2nii release history in reverse chronological order.
 
+0.7.3 (WIP)
+-----------
+- Siemens .rda format now had corrected and validated orientations (tested on VE11 baseline).
+- Siemens .rda format now handles MRSI/CSI data and matches DICOM output. Validated on VE11 baseline data.
+- Fixes in Siemens Twix special case for universal editing sequence (HERMES conditions).
+- Added handling of custom Bruker sequences `mt_sLASER`, `mt_MEGA_sLASER_V35` and `cl_STELASER_PA360_b`.
+
 0.7.2 (Thursday 7th December 2023)
 ----------------------------------
 - SpectralWidth now added to header extension automatically to match bids specification.

README.md

@@ -91,8 +91,7 @@ NIfTI MRS dimension tags (e.g. `DIM_COIL`) can be specified using the `-t` comma
 ### Siemens RDA
 `spec2nii rda RDA_FILE`
 
-Only supports SVS currently. Please contact the developers with examples for MRSI capability.
-Orientation calculation is a WIP, test data for different voxel orientations would be greatly appreciated!
+Compatible with CSI and SVS data. Validated to be the same data and orientation information as DICOM output on VE baselines.
 
 ### UIH DICOM
 `spec2nii uih DCM_FILE_or_DIR`

spec2nii/Siemens/rda.py

@@ -4,7 +4,6 @@
 """
 import re
 from datetime import datetime
-import warnings
 
 import numpy as np
 
@@ -76,25 +75,47 @@ def convert_rda(rda_path, fname_out, verbose):
 
         data = np.fromfile(fp)
 
+    data_cmplx = data[0::2] + 1j * data[1::2]
+
+    # MRSI
     if (int(hdr['CSIMatrixSize[0]'])
             * int(hdr['CSIMatrixSize[1]'])
             * int(hdr['CSIMatrixSize[2]']))\
             > 1:
-        raise MRSINotHandledError('MRSI is currently not handled in the RDA format. Test data needed.')
+        data_cmplx = data_cmplx.reshape((
+            int(hdr['CSIMatrixSize[2]']),
+            int(hdr['CSIMatrixSize[1]']),
+            int(hdr['CSIMatrixSize[0]']),
+            int(hdr['VectorSize'])))
+        data_cmplx = np.moveaxis(data_cmplx, (0, 1, 2), (2, 1, 0))
+        data_shape = data_cmplx.shape[:3]
+
+        imagePositionPatient = np.asarray([
+            _locale_float(hdr['PositionVector[0]']),
+            _locale_float(hdr['PositionVector[1]']),
+            _locale_float(hdr['PositionVector[2]'])])
+
+        half_shift = True
+    # SVS
+    else:
+        data_cmplx = data_cmplx.reshape((1, 1, 1) + data_cmplx.shape)
+        data_shape = (1, 1, 1)
+
+        try:
+            imagePositionPatient = np.asarray([
+                _locale_float(hdr['VOIPositionSag']),
+                _locale_float(hdr['VOIPositionCor']),
+                _locale_float(hdr['VOIPositionTra'])])
+        except KeyError:
+            imagePositionPatient = np.asarray([
+                _locale_float(hdr['PositionVector[0]']),
+                _locale_float(hdr['PositionVector[1]']),
+                _locale_float(hdr['PositionVector[2]'])])
+
+        half_shift = False
 
-    data_cmplx = data[0::2] + 1j * data[1::2]
-    data_cmplx = data_cmplx.reshape((1, 1, 1) + data_cmplx.shape)
     dwelltime = _locale_float(hdr['DwellTime']) / 1E6
 
-    warnings.warn(
-        'The orientation calculations for rda data is mostly untested.'
-        ' Please contribute test data if you can!')
-
-    imagePositionPatient = np.asarray([
-        _locale_float(hdr['PositionVector[0]']),
-        _locale_float(hdr['PositionVector[1]']),
-        _locale_float(hdr['PositionVector[2]'])])
-
     imageOrientationPatient = np.asarray([
         [_locale_float(hdr['RowVector[0]']), _locale_float(hdr['ColumnVector[0]'])],
         [_locale_float(hdr['RowVector[1]']), _locale_float(hdr['ColumnVector[1]'])],
@@ -108,8 +129,9 @@ def convert_rda(rda_path, fname_out, verbose):
     currNiftiOrientation = dcm_to_nifti_orientation(imageOrientationPatient,
                                                     imagePositionPatient,
                                                     xyzMM,
-                                                    (1, 1, 1),
-                                                    verbose=verbose)
+                                                    data_shape,
+                                                    verbose=verbose,
+                                                    half_shift=half_shift)
 
     meta = extractRdaMetadata(hdr)
     meta.set_standard_def('OriginalFile', [rda_path.name])
@@ -137,13 +159,14 @@ def extractRdaMetadata(hdr):
         hdr['Nucleus'])
 
     # Standard defined metadata
-    def set_standard_def(nifti_mrs_key, location, key, cast=None):
-        if key in location\
-                and location[key] is not None:
+    def set_standard_def(nifti_mrs_key, key, cast=None):
+        if key in hdr\
+                and hdr[key] is not None\
+                and hdr[key]:
             if cast is not None:
-                obj.set_standard_def(nifti_mrs_key, cast(location[key]))
+                obj.set_standard_def(nifti_mrs_key, cast(hdr[key]))
             else:
-                obj.set_standard_def(nifti_mrs_key, location[key])
+                obj.set_standard_def(nifti_mrs_key, hdr[key])
 
     # # 5.1 MRS specific Tags
     # 'EchoTime'
@@ -180,22 +203,22 @@ def set_standard_def(nifti_mrs_key, location, key, cast=None):
 
     # # 5.3 Sequence information
     # 'SequenceName'
-    obj.set_standard_def('SequenceName', hdr['SequenceName'])
+    set_standard_def('SequenceName', 'SequenceName')
     # 'ProtocolName'
-    obj.set_standard_def('ProtocolName', hdr['ProtocolName'])
+    set_standard_def('ProtocolName', 'ProtocolName')
     # # 5.4 Sequence information
     # 'PatientPosition'
-    obj.set_standard_def('PatientPosition', hdr['PatientPosition'])
+    set_standard_def('PatientPosition', 'PatientPosition')
     # 'PatientName'
-    obj.set_standard_def('PatientName', hdr['PatientName'])
+    set_standard_def('PatientName', 'PatientName')
     # 'PatientID'
-    obj.set_standard_def('PatientID', hdr['PatientID'])
+    set_standard_def('PatientID', 'PatientID')
     # 'PatientWeight'
-    obj.set_standard_def('PatientWeight', _locale_float(hdr['PatientWeight']))
+    set_standard_def('PatientWeight', 'PatientWeight', cast=_locale_float)
     # 'PatientDoB'
-    obj.set_standard_def('PatientDoB', hdr['PatientBirthDate'])
+    set_standard_def('PatientDoB', 'PatientBirthDate')
     # 'PatientSex'
-    obj.set_standard_def('PatientSex', hdr['PatientSex'])
+    set_standard_def('PatientSex', 'PatientSex')
 
     # # 5.5 Provenance and conversion metadata
     obj.set_standard_def('ConversionMethod', f'spec2nii v{spec2nii_ver}')

spec2nii/Siemens/twix_special_case.py

@@ -142,24 +142,24 @@ def mgs_svs_ed_twix(twixObj, reord_data, meta_obj, dim_tags):
             "OFF": {"PulseOffset": edit_pulse_off, "PulseDuration": pulse_length}}
     elif seq_mode == 1.0:
         # HERMES GABA GSH (3 edit, 1 ctrl condition)
-        edit_cases = 3
+        edit_cases = 4
         dim_info = "HERMES j-difference editing, GABA GSH, four conditions"
         dim_header = {"EditCondition": ["A", "B", "C", "D"]}
         edit_pulse_val = {
             "A": {"PulseOffset": edit_pulse_1, "PulseDuration": 0.02},
-            "B": {"PulseOffset": None, "PulseDuration": None},
-            "C": {"PulseOffset": edit_pulse_off, "PulseDuration": 0.02},
-            "D": {"PulseOffset": [edit_pulse_1, edit_pulse_off], "PulseDuration": 0.02}}
+            "B": {"PulseOffset": edit_pulse_off, "PulseDuration": 0.02},
+            "C": {"PulseOffset": edit_pulse_2, "PulseDuration": 0.02},
+            "D": {"PulseOffset": [edit_pulse_1, edit_pulse_2], "PulseDuration": 0.02}}
     elif seq_mode == 2.0:
         # HERMES GABA GSH EtOH (3 edit, 1 ctrl condition)
         edit_cases = 4
         dim_info = "HERMES j-difference editing, GABA GSH EtOH, four conditions"
         dim_header = {"EditCondition": ["A", "B", "C", "D"]}
         edit_pulse_val = {
-            "A": {"PulseOffset": [edit_pulse_1, edit_pulse_2], "PulseDuration": 0.02},
-            "B": {"PulseOffset": [edit_pulse_3, edit_pulse_2], "PulseDuration": None},
-            "C": {"PulseOffset": [edit_pulse_1, edit_pulse_3], "PulseDuration": 0.02},
-            "D": {"PulseOffset": None, "PulseDuration": None}}
+            "A": {"PulseOffset": [edit_pulse_1, edit_pulse_3], "PulseDuration": 0.02},
+            "B": {"PulseOffset": edit_pulse_off, "PulseDuration": 0.02},
+            "C": {"PulseOffset": [edit_pulse_2, edit_pulse_3], "PulseDuration": 0.02},
+            "D": {"PulseOffset": [edit_pulse_1, edit_pulse_2], "PulseDuration": 0.02}}
     elif seq_mode == 3.0:
         # HERCULES (4 edit conditions)
         edit_cases = 4
@@ -170,16 +170,17 @@ def mgs_svs_ed_twix(twixObj, reord_data, meta_obj, dim_tags):
             "B": {"PulseOffset": [edit_pulse_off, edit_pulse_2], "PulseDuration": 0.02},
             "C": {"PulseOffset": edit_pulse_1, "PulseDuration": 0.02},
             "D": {"PulseOffset": edit_pulse_off, "PulseDuration": 0.02}}
-    elif seq_mode == 3.0:
+    elif seq_mode == 4.0:
+        # This is possibly only a condition for smm_svs_herc as not present in VE11c version.
         # HERMES GABA LAC (3 edit 1 ctrl conditions)
         edit_cases = 4
         dim_info = "HERMES j-difference editing, GABA LAC, four conditions"
         dim_header = {"EditCondition": ["A", "B", "C", "D"]}
         edit_pulse_val = {
             "A": {"PulseOffset": edit_pulse_1, "PulseDuration": 0.02},
-            "B": {"PulseOffset": None, "PulseDuration": None},
-            "C": {"PulseOffset": edit_pulse_off, "PulseDuration": 0.02},
-            "D": {"PulseOffset": [edit_pulse_1, edit_pulse_off], "PulseDuration": 0.02}}
+            "B": {"PulseOffset": edit_pulse_off, "PulseDuration": 0.02},
+            "C": {"PulseOffset": edit_pulse_2, "PulseDuration": 0.02},
+            "D": {"PulseOffset": [edit_pulse_1, edit_pulse_2], "PulseDuration": 0.02}}
     else:
         raise ValueError('Unknown sequence mode in mgs_svs_ed sequence.')
 

tests/test_siemens_rda.py

@@ -17,9 +17,100 @@
 
 xa20_svs_path = siemens_path / 'XAData' / 'XA20/rda/spct_002.MR.MRI-LAB Test_Dir.5.1.114540.rda'
 xa31_locale_svs_path = siemens_path / 'XAData' / 'XA31/rda/locale_XA31.rda'
+xa50_path = siemens_path / 'XAData' / 'XA50' / 'Phantom_20240129.MR.10.1.155227.rda'
 
 latin1_encoding = siemens_path / 'rda' / 'latin1.rda'
 
+ve_path = siemens_path / 'VEData' / 'rda'
+
+paired_ve_data = {
+    'svs_se_c_t15_s10_r10': 'svs_se_c>t15>s10_R10_12_1',
+    'svs_se_iso_tra_sat': 'svs_se_iso_tra_sat_13_1',
+    'svs_se_s_c10_t5_r10': 'svs_se_s>c10>t5_R10_11_1',
+    'svs_se_t_c15_s10_r10': 'svs_se_t>c15>s10_R10_10_1',
+    'csi_se_3d_c_s235_t203_r10': 'csi_se_3D_c>s23.5>t20.3_R10_9_1',
+    'csi_se_3d_s_t235_c203_r10': 'csi_se_3D_s>t23.5>c20.3_R10_8_1',
+    'csi_se_3d_t_c235_s203_r10': 'csi_se_3D_t>c23.5>s20.3_R10_7_1',
+    'csi_se_c_s235_t203_r10': 'csi_se_c>s23.5>t20.3_R10_6_1',
+    'csi_se_iso_tra_sat': 'csi_se_iso_tra_sat_14_1',
+    'csi_se_s_t235_c203_r10': 'csi_se_s>t23.5>c20.3_R10_5_1',
+    'csi_se_t_c235_s203_r10': 'csi_se_t>c23.5>s20.3_R10_4_1'}
+
+
+def test_svs_ve_against_dicom(tmp_path):
+    for rda_file in ve_path.glob("svs_se*.rda"):
+        print(rda_file.stem)
+        dcm_dir = rda_file.parent.parent / 'DICOM' / paired_ve_data[rda_file.stem]
+
+        subprocess.run([
+            'spec2nii', 'rda',
+            '-f', rda_file.stem + '_rda',
+            '-o', tmp_path,
+            '-j', rda_file])
+
+        subprocess.run([
+            'spec2nii', 'dicom',
+            '-f', rda_file.stem + '_dcm',
+            '-o', tmp_path,
+            '-j', dcm_dir])
+
+        rda_out = tmp_path / (rda_file.stem + '_rda.nii.gz')
+        dcm_out = tmp_path / (rda_file.stem + '_dcm.nii.gz')
+        assert rda_out.exists()
+        assert dcm_out.exists()
+
+        dcm = read_nifti_mrs(dcm_out)
+        rda = read_nifti_mrs(rda_out)
+        assert dcm.shape == rda.shape
+        assert np.allclose(dcm.get_fdata(dtype=complex), rda.get_fdata(dtype=complex))
+        assert np.allclose(dcm.get_sform(), rda.get_sform())
+
+        hdr_ext_codes = dcm.header.extensions.get_codes()
+        hdr_ext_dcm = json.loads(dcm.header.extensions[hdr_ext_codes.index(44)].get_content())
+        hdr_ext_codes = rda.header.extensions.get_codes()
+        hdr_ext_rda = json.loads(rda.header.extensions[hdr_ext_codes.index(44)].get_content())
+
+        for key in ['SpectrometerFrequency', 'ResonantNucleus', 'EchoTime', 'RepetitionTime']:
+            assert hdr_ext_dcm[key] == hdr_ext_rda[key]
+
+
+def test_csi_ve_against_dicom(tmp_path):
+    for rda_file in ve_path.glob("csi_se*.rda"):
+        print(rda_file.stem)
+        dcm_dir = rda_file.parent.parent / 'DICOM' / paired_ve_data[rda_file.stem]
+
+        subprocess.run([
+            'spec2nii', 'rda',
+            '-f', rda_file.stem + '_rda',
+            '-o', tmp_path,
+            '-j', rda_file])
+
+        subprocess.run([
+            'spec2nii', 'dicom',
+            '-f', rda_file.stem + '_dcm',
+            '-o', tmp_path,
+            '-j', dcm_dir])
+
+        rda_out = tmp_path / (rda_file.stem + '_rda.nii.gz')
+        dcm_out = tmp_path / (rda_file.stem + '_dcm.nii.gz')
+        assert rda_out.exists()
+        assert dcm_out.exists()
+
+        dcm = read_nifti_mrs(dcm_out)
+        rda = read_nifti_mrs(rda_out)
+        assert dcm.shape == rda.shape
+
+        assert np.allclose(dcm.get_fdata(dtype=complex), rda.get_fdata(dtype=complex))
+        assert np.allclose(dcm.get_sform(), rda.get_sform(), atol=1E-4)
+
+        hdr_ext_codes = dcm.header.extensions.get_codes()
+        hdr_ext_dcm = json.loads(dcm.header.extensions[hdr_ext_codes.index(44)].get_content())
+        hdr_ext_codes = rda.header.extensions.get_codes()
+        hdr_ext_rda = json.loads(rda.header.extensions[hdr_ext_codes.index(44)].get_content())
+
+        for key in ['SpectrometerFrequency', 'ResonantNucleus', 'EchoTime', 'RepetitionTime']:
+            assert hdr_ext_dcm[key] == hdr_ext_rda[key]
+
 
 def test_xa20_svs(tmp_path):
 
@@ -57,6 +148,25 @@ def test_xa31_locale_svs(tmp_path):
     assert np.iscomplexobj(img_t.dataobj)
 
 
+def test_xa50_svs(tmp_path):
+
+    subprocess.run([
+        'spec2nii', 'rda',
+        '-f', 'xa50_svs',
+        '-o', tmp_path,
+        '-j', xa50_path])
+
+    img_t = read_nifti_mrs(tmp_path / 'xa50_svs.nii.gz')
+
+    hdr_ext_codes = img_t.header.extensions.get_codes()
+    hdr_ext = json.loads(img_t.header.extensions[hdr_ext_codes.index(44)].get_content())
+
+    hdr_ext['ResonantNucleus'] = ['1H', ]
+
+    assert img_t.shape == (1, 1, 1, 2048)
+    assert np.iscomplexobj(img_t.dataobj)
+
+
 def test_latin_encoding(tmp_path):
 
     subprocess.check_call(['spec2nii', 'rda',

tests/test_twix.py → tests/test_siemens_twix.py

@@ -23,6 +23,7 @@
 # Special cased data
 hercules_ve = siemens_path / 'HERCULES' / 'Siemens_TIEMO_HERC.dat'
 hercules_xa30 = siemens_path / 'HERCULES' / 'meas_MID02595_FID60346_HERC.dat'
+hermes_xa50 = siemens_path / 'XAData' / 'XA50' / 'smm_svs_herc_v2_hermes.dat'
 
 
 def test_VB(tmp_path):
@@ -202,6 +203,27 @@ def test_XA_HERCULES(tmp_path):
     assert hdr_ext['dim_7'] == 'DIM_DYN'
 
 
+def test_XA_HERMES(tmp_path):
+
+    subprocess.check_call(['spec2nii', 'twix',
+                           '-e', 'image',
+                           '-f', 'hermes_xa',
+                           '-o', tmp_path,
+                           '-j', str(hermes_xa50)])
+
+    img_t = read_nifti_mrs(tmp_path / 'hermes_xa.nii.gz')
+
+    hdr_ext_codes = img_t.header.extensions.get_codes()
+    hdr_ext = json.loads(img_t.header.extensions[hdr_ext_codes.index(44)].get_content())
+
+    assert img_t.shape == (1, 1, 1, 4096, 42, 4, 80)
+    assert np.iscomplexobj(img_t.dataobj)
+
+    assert hdr_ext['dim_5'] == 'DIM_COIL'
+    assert hdr_ext['dim_6'] == 'DIM_EDIT'
+    assert hdr_ext['dim_7'] == 'DIM_DYN'
+
+
 def test_twix_mrsi_orientation(tmp_path):
     '''Test that the (empty) mrsi has information matching the twix equivalent.'''