From f63a41871f421a10484c1fbf85b7caeb97062c03 Mon Sep 17 00:00:00 2001 From: William T Clarke Date: Tue, 12 Mar 2024 15:08:41 +0000 Subject: [PATCH] Add handling and tests for philips mega sequence. (#131) --- CHANGELOG.md | 5 ++- README.md | 4 +- notes/philips.md | 21 +++++++++ spec2nii/Philips/philips_dcm.py | 78 +++++++++++++++++++++++++-------- spec2nii/spec2nii.py | 2 +- tests/spec2nii_test_data | 2 +- tests/test_philips_dicom.py | 65 ++++++++++++++++++++++++++- 7 files changed, 152 insertions(+), 25 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a25cc4e..53d8458 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,11 +1,12 @@ This document contains the Spec2nii release history in reverse chronological order. -0.7.3 (WIP) ------------ +0.7.3 (Tuesday 12th March 2024) +------------------------------- - 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`. +- Philips vendor MEGA-PRESS handled through DICOM pathway. Thanks to Sandeep Ganji and Yansong Zhao for their help. 0.7.2 (Thursday 7th December 2023) ---------------------------------- diff --git a/README.md b/README.md index bee54f8..24f94ab 100644 --- a/README.md +++ b/README.md @@ -43,7 +43,7 @@ This table lists the currently supported formats. I have very limited experience |---------------|----------------|-----|-----|-----------------------| | Siemens Twix | .dat | Yes | † | Yes | | Siemens DICOM | .ima / .dcm | Yes | Yes | Yes | -| Siemens RDA | .rda | Yes | No | Yes (WIP) | +| Siemens RDA | .rda | Yes | Yes | Yes | | Philips | .SPAR/.SDAT | Yes | No | Yes | | Philips | .data/.list | Yes | No | Yes | | Philips DICOM | .dcm | Yes | No | Yes | @@ -119,7 +119,7 @@ NIfTI MRS dimension tags (e.g. `DIM_COIL`) can be specified using the `-t` comma Generates separate reference file if present. -Both classic and enhanced DICOM is handled. +Both classic and enhanced DICOM is handled. Well tested on the vendors' own PRESS and MEGA-PRESS sequence. ### Bruker (2dseq/fid) `spec2nii bruker -m 2DSEQ 2DSEQ_FILE_or_DIR` diff --git a/notes/philips.md b/notes/philips.md index fd47016..7457083 100644 --- a/notes/philips.md +++ b/notes/philips.md @@ -12,3 +12,24 @@ The newer format `MR Spectroscopy Storage` can be identified by the `MediaStorag There are many differences between the format of the two types, not least data storage location, header structure, and the use of multiple files versus one to store the complete data set in. For DICOM tag interpretation in the old format see https://github.com/malaterre/dicom-private-dicts/blob/master/PMS-R32-dict.txt + +### Header description for enhanced + - (2005,1597) – tells if it’s an editing sequence + - (2005,1598) – describes the editing type of the spectrum (0 indicates OFF, 1 indicates ON) + - (2005,1304) – spectrum_mix_no (tells if its water data) + - (2005,1357) – Sample Frequency (BW / spectral width) + - (2001,1083) – Imaging Frequency (main field resonance frequency in MHz) + - (2005,1310) – Spectrum Echo time + - (0018,0080) – Repetition Time + + - (2005,1054) – volume_angulation_ap + - (2005,1055) – volume_angulation_fh + - (2005,1056) – volume_angulation_rl + - (2005,105A) – volume_offcentre_ap + - (2005,105B) – volume_offcentre_fh + - (2005,105C) – volume_offcentre_rl + +Under (2005,1085) you can find the FOV size + - (2005,1057) – FOV in AP + - (2005,1058) – FOV in FH + - (2005,1059) – FOV in RL \ No newline at end of file diff --git a/spec2nii/Philips/philips_dcm.py b/spec2nii/Philips/philips_dcm.py index 1952fd9..79127a2 100644 --- a/spec2nii/Philips/philips_dcm.py +++ b/spec2nii/Philips/philips_dcm.py @@ -84,6 +84,7 @@ def multi_file_dicom(files_in, fname_out, tag, verbose): newshape = (1, 1, 1) + specDataCmplx.shape spec_data = specDataCmplx.reshape(newshape) if spec_ref is not None: + newshape = (1, 1, 1) + spec_ref.shape spec_ref = spec_ref.reshape(newshape) # Data appears to require conjugation to meet standard's conventions. @@ -210,12 +211,9 @@ def _process_philips_svs_new(img, verbose): specData = np.frombuffer(img.dcm_data[('5600', '0020')].value, dtype=np.single) specDataCmplx = specData[0::2] + 1j * specData[1::2] - # In the one piece of data I have been provided the data is twice as long as indicated (1 avg) - # and the second half is a water reference. # (0028,0008) IS [32] # 2,1 Number of Frames nframes = int(img.dcm_data[0x0028, 0x0008].value) spec_points = img.dcm_data.SpectroscopyAcquisitionDataColumns - reference_inlcuded = img.dcm_data[0x0018, 0x9199].value == 'YES' npoints_expected = spec_points * nframes if npoints_expected != specDataCmplx.size: @@ -223,25 +221,69 @@ def _process_philips_svs_new(img, verbose): f'Number of points {specDataCmplx.size}, does not equal frames ({nframes}) ' f'* spectral points ({spec_points}).') - if reference_inlcuded: - spec_data_main = specDataCmplx[:spec_points] - spec_data_ref = specDataCmplx[spec_points:] - spec_data_main = spec_data_main.reshape(nframes - 1, spec_points).T - spec_data_main = spec_data_main.squeeze() - else: - spec_data_main = specDataCmplx.reshape(nframes, spec_points).T - spec_data_main = spec_data_main.squeeze() - spec_data_ref = None + # Only two conditions are known for the Spectroscopy scans 'normal/SV_PRESS' and 'SV_MEGA' + # So special case the MEGA and everything else goes through the PRESS pipeline. + try: + is_edited = img.dcm_data[0x2005, 0x1597].value == 'Y' + except KeyError: + is_edited = False + + if img.dcm_data[0x0008, 0x9209].value == 'SPECTROSCOPY'\ + and is_edited: + # Data is MEGA + on_idx = [] + off_idx = [] + ref_idx = [] + for idx, instance in enumerate(img.dcm_data[0x2005, 0x140f]): + if instance[0x2005, 0x1304].value: + ref_idx.append(idx) + else: + if instance[0x2005, 0x1598].value == '1': + on_idx.append(idx) + elif instance[0x2005, 0x1598].value == '0': + off_idx.append(idx) - # 1) Extract dicom parameters - currNiftiOrientation = _enhanced_dcm_svs_to_orientation(img, verbose) + specDataCmplx = specDataCmplx.reshape((nframes, spec_points)) - dwelltime = 1.0 / img.dcm_data.SpectralWidth - meta = _extractDicomMetadata_new(img) - if reference_inlcuded: + spec_data_main = np.stack( + (specDataCmplx[on_idx, :].T, + specDataCmplx[off_idx, :].T), axis=-1) + spec_data_ref = specDataCmplx[ref_idx, :].T + + meta = _extractDicomMetadata_new(img) meta_r = _extractDicomMetadata_new(img, water_suppressed=False) + + meta.set_dim_info("5th", "DIM_DYN") + meta.set_dim_info("6th", "DIM_EDIT", hdr={"EditCondition": ["ON", "OFF"]}) + + meta_r.set_dim_info("5th", "DIM_DYN") + else: - meta_r = None + # In the one piece of data I have been provided the data is twice as long as indicated (1 avg) + # and the second half is a water reference. + + # This actually tells you if it is selected in the UI, but seems to be a surrogate + # for acquired in the PRESS scan + reference_inlcuded = img.dcm_data[0x0018, 0x9199].value == 'YES' + + if reference_inlcuded: + spec_data_main = specDataCmplx[:spec_points] + spec_data_ref = specDataCmplx[spec_points:] + spec_data_main = spec_data_main.reshape(nframes - 1, spec_points).T + spec_data_main = spec_data_main.squeeze() + else: + spec_data_main = specDataCmplx.reshape(nframes, spec_points).T + spec_data_main = spec_data_main.squeeze() + spec_data_ref = None + + meta = _extractDicomMetadata_new(img) + if reference_inlcuded: + meta_r = _extractDicomMetadata_new(img, water_suppressed=False) + else: + meta_r = None + + currNiftiOrientation = _enhanced_dcm_svs_to_orientation(img, verbose) + dwelltime = 1.0 / img.dcm_data.SpectralWidth return spec_data_main, spec_data_ref, currNiftiOrientation, dwelltime, meta, meta_r diff --git a/spec2nii/spec2nii.py b/spec2nii/spec2nii.py index 86f78d9..5e34094 100644 --- a/spec2nii/spec2nii.py +++ b/spec2nii/spec2nii.py @@ -620,7 +620,7 @@ def philips_dicom(self, args): """Philips DICOM format handler.""" from warnings import warn from spec2nii.Philips.philips_dcm import multi_file_dicom - warn('This Philips DICOM conversion routine is experimental and poorly tested.' + warn('This Philips DICOM conversion routine has limited testing outside vendor PRESS and MEGA sequences.' ' Please get in contact with test data to help improve it.') path_in = Path(args.file) if path_in.is_dir(): diff --git a/tests/spec2nii_test_data b/tests/spec2nii_test_data index d300d59..74415ef 160000 --- a/tests/spec2nii_test_data +++ b/tests/spec2nii_test_data @@ -1 +1 @@ -Subproject commit d300d59ce34eac29d6c24499f87f640f9848864d +Subproject commit 74415ef6bb9556384cd685e766c0e4b47998902b diff --git a/tests/test_philips_dicom.py b/tests/test_philips_dicom.py index 6015a06..608f167 100644 --- a/tests/test_philips_dicom.py +++ b/tests/test_philips_dicom.py @@ -17,6 +17,9 @@ data_path = philips_path / 'DICOM' / 'SV_phantom_H15mm' / 'IM-0020-0002-0001.dcm' data_path_2 = philips_path / 'DICOM_enhanced_multi_dynamic' / 'svsWSAntCing_S002' / 'XX_0004' +data_path_press = philips_path / 'DICOM_enhanced_multi_dynamic' / 'press_mega' / 'XX_0006' +data_path_mpress = philips_path / 'DICOM_enhanced_multi_dynamic' / 'press_mega' / 'XX_0004' + def test_svs(tmp_path): @@ -42,7 +45,7 @@ def test_svs(tmp_path): assert (tmp_path / 'svs_ref.nii.gz').is_file() -def test_svs2(tmp_path): +def test_svs_enhanced(tmp_path): subprocess.check_call(['spec2nii', 'philips_dcm', '-f', 'svs', @@ -62,3 +65,63 @@ def test_svs2(tmp_path): assert hdr_ext['SpectrometerFrequency'][0] == 127.765536 assert hdr_ext['ResonantNucleus'][0] == '1H' assert hdr_ext['OriginalFile'][0] == data_path_2.name + + # Similar data from another test set + subprocess.check_call(['spec2nii', 'philips_dcm', + '-f', 'svs2', + '-o', tmp_path, + '-j', + str(data_path_press)]) + + img_metab = read_nifti_mrs(tmp_path / 'svs2.nii.gz') + img_ref = read_nifti_mrs(tmp_path / 'svs2_ref.nii.gz') + + assert img_metab.shape == (1, 1, 1, 1024) + assert np.iscomplexobj(img_metab.dataobj) + assert np.isclose(1 / img_metab.header['pixdim'][4], 2000.0) + + hdr_ext_codes = img_metab.header.extensions.get_codes() + hdr_ext = json.loads(img_metab.header.extensions[hdr_ext_codes.index(44)].get_content()) + + assert hdr_ext['SpectrometerFrequency'][0] == 127.770768 + assert hdr_ext['ResonantNucleus'][0] == '1H' + assert hdr_ext['OriginalFile'][0] == data_path_press.name + assert hdr_ext['WaterSuppressed'] is True + + assert img_ref.shape == (1, 1, 1, 1024) + hdr_ext_codes = img_ref.header.extensions.get_codes() + hdr_ext = json.loads(img_ref.header.extensions[hdr_ext_codes.index(44)].get_content()) + assert hdr_ext['WaterSuppressed'] is False + + +def test_svs_mega(tmp_path): + subprocess.run([ + 'spec2nii', 'philips_dcm', + '-f', 'svs_mega', + '-o', tmp_path, + '-j', + data_path_mpress]) + + img_metab = read_nifti_mrs(tmp_path / 'svs_mega.nii.gz') + img_ref = read_nifti_mrs(tmp_path / 'svs_mega_ref.nii.gz') + + assert img_metab.shape == (1, 1, 1, 1024, 144, 2) + assert np.iscomplexobj(img_metab.dataobj) + assert np.isclose(1 / img_metab.header['pixdim'][4], 2000.0) + + hdr_ext_codes = img_metab.header.extensions.get_codes() + hdr_ext = json.loads(img_metab.header.extensions[hdr_ext_codes.index(44)].get_content()) + + assert hdr_ext['SpectrometerFrequency'][0] == 127.770768 + assert hdr_ext['ResonantNucleus'][0] == '1H' + assert hdr_ext['OriginalFile'][0] == data_path_mpress.name + assert hdr_ext['WaterSuppressed'] is True + assert hdr_ext['dim_5'] == 'DIM_DYN' + assert hdr_ext['dim_6'] == 'DIM_EDIT' + assert hdr_ext['dim_6_header'] == {'EditCondition': ['ON', 'OFF']} + + assert img_ref.shape == (1, 1, 1, 1024, 9) + hdr_ext_codes = img_ref.header.extensions.get_codes() + hdr_ext = json.loads(img_ref.header.extensions[hdr_ext_codes.index(44)].get_content()) + assert hdr_ext['WaterSuppressed'] is False + assert hdr_ext['dim_5'] == 'DIM_DYN'