diff --git a/CHANGELOG.md b/CHANGELOG.md index 23edd76..6e9ca60 100644 --- a/CHANGELOG.md +++ b/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) ---------------------------------- diff --git a/spec2nii/Siemens/rda.py b/spec2nii/Siemens/rda.py index 11ed44d..75b2247 100644 --- a/spec2nii/Siemens/rda.py +++ b/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]) diff --git a/tests/test_siemens_rda.py b/tests/test_siemens_rda.py index 648e5d1..e336a5f 100644 --- a/tests/test_siemens_rda.py +++ b/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',