Skip to content

Commit

Permalink
ENH: Validated rda orientation and CSI handling (#128)
Browse files Browse the repository at this point in the history 
* Validated svs rda orientations.

* Enable rda CSI

* Update readme

* Handle lack of VOIPositionSag in XA20 data.
  • Loading branch information
@wtclarke
wtclarke authored Mar 2, 2024
1 parent b915930 commit f5a0c20
Show file tree
Hide file tree
Showing 5 changed files with 148 additions and 31 deletions.
5 changes: 5 additions & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,10 @@
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.

0.7.2 (Thursday 7th December 2023)
----------------------------------
- SpectralWidth now added to header extension automatically to match bids specification.
Expand Down
3 changes: 1 addition & 2 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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`
Expand Down
79 changes: 51 additions & 28 deletions spec2nii/Siemens/rda.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,6 @@
"""
import re
from datetime import datetime
import warnings

import numpy as np

Expand Down Expand Up @@ -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]'])],
Expand All @@ -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])
Expand Down Expand Up @@ -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'
Expand Down Expand Up @@ -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}')
Expand Down
2 changes: 1 addition & 1 deletion tests/spec2nii_test_data
Submodule spec2nii_test_data updated from ca9850 to 9322e3
90 changes: 90 additions & 0 deletions tests/test_siemens_rda.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,96 @@

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):

Expand Down

0 comments on commit f5a0c20

Please sign in to comment.