Enable rda CSI

CHANGELOG.md

@@ -3,7 +3,7 @@ This document contains the Spec2nii release history in reverse chronological ord
 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.
 
 0.7.2 (Thursday 7th December 2023)
 ----------------------------------

spec2nii/Siemens/rda.py

@@ -76,24 +76,40 @@ 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)
 
-    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
+        imagePositionPatient = np.asarray([
+            _locale_float(hdr['VOIPositionSag']),
+            _locale_float(hdr['VOIPositionCor']),
+            _locale_float(hdr['VOIPositionTra'])])
 
-    warnings.warn(
-        'The orientation calculations for rda data is mostly untested.'
-        ' Please contribute test data if you can!')
+        half_shift = False
 
-    imagePositionPatient = np.asarray([
-        _locale_float(hdr['VOIPositionSag']),
-        _locale_float(hdr['VOIPositionCor']),
-        _locale_float(hdr['VOIPositionTra'])])
+    dwelltime = _locale_float(hdr['DwellTime']) / 1E6
 
     imageOrientationPatient = np.asarray([
         [_locale_float(hdr['RowVector[0]']), _locale_float(hdr['ColumnVector[0]'])],
@@ -108,8 +124,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])

tests/test_siemens_rda.py

@@ -21,14 +21,22 @@
 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'}
+    '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_ve_against_dicom(tmp_path):
+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]
@@ -65,6 +73,44 @@ def test_ve_against_dicom(tmp_path):
             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):
 
     subprocess.check_call(['spec2nii', 'rda',