From d3add6804aa0ef3264bf77042f88a0720badd821 Mon Sep 17 00:00:00 2001 From: William T Clarke Date: Thu, 18 Apr 2024 16:13:30 +0100 Subject: [PATCH] ENH/svs_slaser_dkd (#136) * Handling for slaser_dkd sequences DICOM format * Add twix handling of slaser_dkd sequences --- CHANGELOG.md | 5 +- spec2nii/Siemens/dicomfunctions.py | 44 ++++++++++++-- spec2nii/Siemens/twix_special_case.py | 87 +++++++++++++++++++++++++++ spec2nii/Siemens/twixfunctions.py | 42 +++++++++---- spec2nii/spec2nii.py | 9 ++- tests/spec2nii_test_data | 2 +- tests/test_siemens_special_data.py | 59 ++++++++++++++++++ 7 files changed, 225 insertions(+), 23 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d50924c..41feb85 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,10 +1,11 @@ This document contains the Spec2nii release history in reverse chronological order. -0.7.4 (WIP) ---------------------------------- +0.7.4 (Thursday 18th April 2024) +-------------------------------- - Refinements and improvements to the GE SVS pipeline from Mark Mikkelsen. - Add support for older jMRUI text formats which have a slightly different syntax. With thanks to Donnie Cameron. - Handle odd case of XA like .twix headers in a VX baseline scan +- Improved (and validated) handling of Dinesh Deelchand's slaser sequences with integrated references (`svs_slaser(voi)_dkd(2)`) 0.7.3 (Tuesday 12th March 2024) ------------------------------- diff --git a/spec2nii/Siemens/dicomfunctions.py b/spec2nii/Siemens/dicomfunctions.py index 89d2817..7e3ad91 100644 --- a/spec2nii/Siemens/dicomfunctions.py +++ b/spec2nii/Siemens/dicomfunctions.py @@ -719,24 +719,56 @@ def identify_integrated_references(img, inst_num): # Handle CMRR DKD sequence # https://www.cmrr.umn.edu/spectro/ # SEMI-LASER (MRM 2011, NMB 2019) Release 2016-12 + # Three cases as explained by DKD. + ''' + dkd_slaserVOI_dkd + ScanMode = 8 + (NOTE DKD originally said 2, but from testing the sLASER_VE11C_Sept2016 version 8 is the only option + xprot[('sSpecPara', 'lAutoRefScanMode')] = 1 when set to 'off' and 8 when set to 'on') + 4 water ref at start and 4 at the end where the first 2 water ref in + each are acquired with VAPOR RF off but OVS on while last 2 ref are with VAPOR and OVS completely off. + + dkd_slaserVOI_dkd2 + ScanMode =8 + does the same as ScanMode =8 above + + ScanMode =2 with ScanNo=2 + this means 2 water ref at start and 2 at end; + this acq is done with VAPOR and OVS RF off, but all gradients still ON. + ''' seq_file_name = xprot[('tSequenceFileName',)].strip('"').lower() - match = re.search(r'svs_slaservoi_dkd', seq_file_name) - if match and xprot[('sSpecPara', 'lAutoRefScanMode')] == 8.0: + if (re.search(r'svs_slaser(voi)?_dkd2$', seq_file_name) + and xprot[('sSpecPara', 'lAutoRefScanMode')] == 8.0)\ + or (re.search(r'svs_slaser(voi)?_dkd$', seq_file_name) + and xprot[('sSpecPara', 'lAutoRefScanMode')] == 8.0): num_ref = int(xprot[('sSpecPara', 'lAutoRefScanNo')]) num_dyn = int(xprot[('lAverages',)]) total_dyn = num_dyn + (num_ref * 4) if inst_num <= num_ref: # First ecc calibration references - return 1, '_ecc' + return 1, '_rf_off' elif inst_num <= (num_ref * 2): # First quantitation calibration references - return 2, '_quant' + return 2, '_rf_grads_ovs_off' elif (total_dyn - (2 * num_ref)) < inst_num <= (total_dyn - num_ref): # Second ecc calibration references - return 1, '_ecc' + return 1, '_rf_off' elif (total_dyn - num_ref) < inst_num <= total_dyn: # Second quantitation calibration references - return 2, '_quant' + return 2, '_rf_grads_ovs_off' + else: + return 0, '' + elif (re.search(r'svs_slaser(voi)?_dkd2', seq_file_name) + and xprot[('sSpecPara', 'lAutoRefScanMode')] == 2.0): + num_ref = int(xprot[('sSpecPara', 'lAutoRefScanNo')]) + num_dyn = int(xprot[('lAverages',)]) + total_dyn = num_dyn + (num_ref * 2) + if inst_num < num_ref: + # First WS and OVS RF off + return 1, '_vapor_ovs_rfoff' + elif inst_num >= (total_dyn - num_ref): + # SecondWS and OVS RF off + return 1, '_vapor_ovs_rfoff' else: return 0, '' else: diff --git a/spec2nii/Siemens/twix_special_case.py b/spec2nii/Siemens/twix_special_case.py index 90a1e42..7f96e22 100644 --- a/spec2nii/Siemens/twix_special_case.py +++ b/spec2nii/Siemens/twix_special_case.py @@ -4,6 +4,7 @@ Will Clarke, University of Oxford, 2022 """ import numpy as np +import re def smm_svs_herc_hyper(twixObj, reord_data, meta_obj, dim_tags, subseq, subseq_name): @@ -204,3 +205,89 @@ def mgs_svs_ed_twix(twixObj, reord_data, meta_obj, dim_tags): meta_obj.set_standard_def("EditPulse", edit_pulse_val) return reord_data, meta_obj, dim_tags + + +def slaser_dkd(twixObj, reord_data, meta_obj, dim_tags): + + scan_mode = twixObj.hdr['MeasYaps'][('sSpecPara', 'lAutoRefScanMode')] + seq_file_name = twixObj['hdr']['Meas'][('tSequenceFileName')].lower() + + group_out = [] + + if scan_mode == 1.0: + for idx, dt in enumerate(dim_tags): + meta_obj.set_dim_info(idx, dt) + group_out.append( + (reord_data.reshape((1, 1, 1) + reord_data.shape), + meta_obj, + '')) + + elif (re.search(r'svs_slaser(voi)?_dkd2$', seq_file_name) and scan_mode == 8.0)\ + or (re.search(r'svs_slaser(voi)?_dkd$', seq_file_name) and scan_mode == 8.0): + num_ref = int(twixObj.hdr['MeasYaps'][('sSpecPara', 'lAutoRefScanNo')]) + num_dyn = int(twixObj.hdr['MeasYaps'][('lAverages',)]) + + base_shape = (1, 1, 1) + reord_data.shape[:2] + (-1, ) + + for idx, dt in enumerate(dim_tags): + meta_obj.set_dim_info(idx, dt) + + # Main acq + start = (num_ref * 2) + stop = start + num_dyn + group_out.append( + (reord_data[..., start:stop].reshape(base_shape), + meta_obj.copy(), + '')) + + # _rf_off + indices = np.concatenate( + (np.arange(num_ref), np.arange(-2 * num_ref, -num_ref))) + group_out.append( + (reord_data[..., indices].reshape(base_shape), + meta_obj.copy(), + '_rf_off')) + + # _rf_grads_ovs_off + indices = np.concatenate( + (np.arange(num_ref, 2 * num_ref), np.arange(-num_ref, 0))) + group_out.append( + (reord_data[..., indices].reshape(base_shape), + meta_obj.copy(), + '_rf_grads_ovs_off')) + + elif (re.search(r'svs_slaser(voi)?_dkd2', seq_file_name) and scan_mode == 2.0): + num_ref = int(twixObj.hdr['MeasYaps'][('sSpecPara', 'lAutoRefScanNo')]) + num_dyn = int(twixObj.hdr['MeasYaps'][('lAverages',)]) + + base_shape = (1, 1, 1) + reord_data.shape[:2] + (-1, ) + + for idx, dt in enumerate(dim_tags): + meta_obj.set_dim_info(idx, dt) + + # Main acq + start = num_ref + stop = start + num_dyn + group_out.append( + (reord_data[..., start:stop].reshape(base_shape), + meta_obj.copy(), + '')) + + # _vapor_ovs_rfoff + indices = np.concatenate( + (np.arange(num_ref), np.arange(-num_ref, 0))) + group_out.append( + (reord_data[..., indices].reshape(base_shape), + meta_obj.copy(), + '_vapor_ovs_rfoff')) + + else: + print(f'slaser_dkd special case processing: unrecognised case for {seq_file_name} with scan mode {scan_mode}.') + for idx, dt in enumerate(dim_tags): + meta_obj.set_dim_info(idx, dt) + group_out.append( + (reord_data.reshape((1, 1, 1) + reord_data.shape), + meta_obj, + '')) + + return group_out diff --git a/spec2nii/Siemens/twixfunctions.py b/spec2nii/Siemens/twixfunctions.py index e5c3b0c..dcef172 100644 --- a/spec2nii/Siemens/twixfunctions.py +++ b/spec2nii/Siemens/twixfunctions.py @@ -305,7 +305,7 @@ def process_svs(twixObj, base_name_out, name_in, dataKey, dim_overrides, remove_ # AG 03/22/2023 Moved the NIFTI/Out lists and if statement up to work with new Hyper function below. # AG 03/22/2023 Create Lists for assembled data - Starts with XA reference.. - nifit_mrs_out = [] + nifti_mrs_out = [] filename_out = [] if xa_ref_scans is not None: @@ -319,7 +319,7 @@ def process_svs(twixObj, base_name_out, name_in, dataKey, dim_overrides, remove_ meta_obj_ref = extractTwixMetadata(twixObj['hdr'], basename(twixObj[dataKey].filename)) meta_obj_ref.set_standard_def('WaterSuppressed', True) - nifit_mrs_out.append( + nifti_mrs_out.append( assemble_nifti_mrs( xa_ref_scans.reshape(newshape), dwellTime, @@ -395,19 +395,37 @@ def process_svs(twixObj, base_name_out, name_in, dataKey, dim_overrides, remove_ if reord_data_.ndim <= 4: # 4 or less Dims newshape = (1, 1, 1) + reord_data_.shape # Pad 3 singleton dims - nifit_mrs_out.append(assemble_nifti_mrs(reord_data_.reshape(newshape), + nifti_mrs_out.append(assemble_nifti_mrs(reord_data_.reshape(newshape), dwellTime, orientation, meta_obj_)) filename_out.append(f'{mainStr}_{hyp_suffix[ii]}') - return nifit_mrs_out, filename_out + return nifti_mrs_out, filename_out # HERCULES Data elif (xa_or_vx(twixObj['hdr']) == 'xa' and 'smm_svs_herc' in twixObj['hdr']['Meas'][('tSequenceFileName')]): from spec2nii.Siemens.twix_special_case import mgs_svs_ed_twix reord_data, meta_obj, dim_tags = mgs_svs_ed_twix(twixObj, reord_data, meta_obj, dim_tags) + elif re.search( + r'svs_slaser(voi)?_dkd', + twixObj['hdr']['Meas'][('tSequenceFileName')], + re.IGNORECASE): + from spec2nii.Siemens.twix_special_case import slaser_dkd + for data, meta, name in slaser_dkd(twixObj, reord_data, meta_obj, dim_tags): + + nifti_mrs_out.append( + gen_nifti_mrs_hdr_ext( + data, + dwellTime, + meta, + orientation.Q44, + no_conj=True)) + + filename_out.append(f'{mainStr}{name}') + + return nifti_mrs_out, filename_out else: # Set dim tags in meta now as no additional info for idx, dt in enumerate(dim_tags): @@ -417,7 +435,7 @@ def process_svs(twixObj, base_name_out, name_in, dataKey, dim_overrides, remove_ # Pad with three singleton dimensions (x,y,z) newshape = (1, 1, 1) + reord_data.shape - nifit_mrs_out.append(assemble_nifti_mrs(reord_data.reshape(newshape), + nifti_mrs_out.append(assemble_nifti_mrs(reord_data.reshape(newshape), dwellTime, orientation, meta_obj)) @@ -432,7 +450,7 @@ def process_svs(twixObj, base_name_out, name_in, dataKey, dim_overrides, remove_ # Pad with three singleton dimensions (x,y,z) newshape = (1, 1, 1) + reord_data[modIndex].shape - nifit_mrs_out.append( + nifti_mrs_out.append( assemble_nifti_mrs(reord_data[modIndex].reshape(newshape), dwellTime, orientation, @@ -446,7 +464,7 @@ def process_svs(twixObj, base_name_out, name_in, dataKey, dim_overrides, remove_ filename_out.append(out_name) - return nifit_mrs_out, filename_out + return nifti_mrs_out, filename_out def process_fid(twixObj, base_name_out, name_in, dataKey, dim_overrides, remove_os, quiet=False, verbose=False): @@ -580,13 +598,13 @@ def process_fid(twixObj, base_name_out, name_in, dataKey, dim_overrides, remove_ reord_data = np.moveaxis(squeezedData, original, new) # Now assemble data - nifit_mrs_out = [] + nifti_mrs_out = [] filename_out = [] if reord_data.ndim <= 4: # Pad with three singleton dimensions (x,y,z) newshape = (1, 1, 1) + reord_data.shape - nifit_mrs_out.append(assemble_nifti_mrs(reord_data.reshape(newshape), + nifti_mrs_out.append(assemble_nifti_mrs(reord_data.reshape(newshape), dwellTime, orientation, meta_obj, @@ -602,7 +620,7 @@ def process_fid(twixObj, base_name_out, name_in, dataKey, dim_overrides, remove_ # Pad with three singleton dimensions (x,y,z) newshape = (1, 1, 1) + reord_data[modIndex].shape - nifit_mrs_out.append( + nifti_mrs_out.append( assemble_nifti_mrs(reord_data[modIndex].reshape(newshape), dwellTime, orientation, @@ -628,7 +646,7 @@ def process_fid(twixObj, base_name_out, name_in, dataKey, dim_overrides, remove_ meta_obj_ref = extractTwixMetadata(twixObj['hdr'], basename(twixObj[dataKey].filename)) meta_obj_ref.set_standard_def('WaterSuppressed', True) - nifit_mrs_out.append( + nifti_mrs_out.append( assemble_nifti_mrs( xa_ref_scans.reshape(newshape), dwellTime, @@ -638,7 +656,7 @@ def process_fid(twixObj, base_name_out, name_in, dataKey, dim_overrides, remove_ filename_out.append(mainStr + '_ref') - return nifit_mrs_out, filename_out + return nifti_mrs_out, filename_out def assemble_nifti_mrs(data, dwellTime, orientation, meta_obj, dim_tags=None): diff --git a/spec2nii/spec2nii.py b/spec2nii/spec2nii.py index 5e34094..d06c93e 100644 --- a/spec2nii/spec2nii.py +++ b/spec2nii/spec2nii.py @@ -80,8 +80,13 @@ def add_common_parameters(subparser): # Handle dicom subcommand parser_dicom = subparsers.add_parser('dicom', help='Convert from Siemens DICOM format.') parser_dicom.add_argument('file', help='file or directory to convert', type=str) - parser_dicom.add_argument("-t", "--tag", type=str, help="Specify NIfTI MRS tag used for 5th " - "dimension if multiple files are passed.") + parser_dicom.add_argument( + "-t", + "--tag", + type=str, + default='DIM_DYN', + help="Specify NIfTI MRS tag used for 5th dimension if multiple files are passed. " + "Defaults to DIM_DYN.") parser_dicom.add_argument('--voi', action='store_true', help='Output VOI as single voxel NIfTI mask.') parser_dicom = add_common_parameters(parser_dicom) parser_dicom.set_defaults(func=self.dicom) diff --git a/tests/spec2nii_test_data b/tests/spec2nii_test_data index 12973b9..8631cdc 160000 --- a/tests/spec2nii_test_data +++ b/tests/spec2nii_test_data @@ -1 +1 @@ -Subproject commit 12973b964ce06faf4f1a1df1a522d87e2a407067 +Subproject commit 8631cdc0f22abead50a22f202ca6a290629de569 diff --git a/tests/test_siemens_special_data.py b/tests/test_siemens_special_data.py index d7769e8..64edd8b 100644 --- a/tests/test_siemens_special_data.py +++ b/tests/test_siemens_special_data.py @@ -61,3 +61,62 @@ def test_dicom_anon(tmp_path): assert hdr_ext['SpectrometerFrequency'][0] == 123.255582 assert hdr_ext['ResonantNucleus'][0] == "1H" assert hdr_ext['OriginalFile'][0] == str(data_path_anon) + + +data_path_slaser_dkd_dicom = siemens_path / 'special_cases_slaser_dkd' / 'DICOM' +slaser_dkd_options = { + "svs_slaser_dkd_von_wrsoff_13_None": {'': 4, '_rf_off': 0, '_rf_grads_ovs_off': 0}, + "svs_slaser_dkd_von_wrs1_15_None": {'': 4, '_rf_off': 2, '_rf_grads_ovs_off': 2}, + "svs_slaser_dkd_von_wrs2_17_None": {'': 4, '_rf_off': 4, '_rf_grads_ovs_off': 4}, + "svs_slaserVOI_dkd2_von_wrsoff_19_None": {'': 4, '_vapor_ovs_rfoff': 0}, + "svs_slaserVOI_dkd2_von_wrsw1pw3_1_21_None": {'': 4, '_rf_off': 2, '_rf_grads_ovs_off': 2}, + "svs_slaserVOI_dkd2_von_wrsw1pw3_2_23_None": {'': 4, '_rf_off': 4, '_rf_grads_ovs_off': 4}, + "svs_slaserVOI_dkd2_von_wrsw4_1_25_None": {'': 4, '_vapor_ovs_rfoff': 2}, + "svs_slaserVOI_dkd2_von_wrsw4_2_27_None": {'': 4, '_vapor_ovs_rfoff': 4} +} + + +def test_dicom_slaser_dkd(tmp_path): + for key in slaser_dkd_options: + subprocess.run( + ['spec2nii', 'dicom', + '-f', key, + '-o', tmp_path, + '-j', data_path_slaser_dkd_dicom / key]) + + for subkey in slaser_dkd_options[key]: + if slaser_dkd_options[key][subkey] > 0: + print(f'{key}, {subkey}, {slaser_dkd_options[key][subkey]}') + img = read_nifti_mrs(tmp_path / f'{key}{subkey}.nii.gz') + assert img.shape == (1, 1, 1, 2048, slaser_dkd_options[key][subkey]) + + +data_path_slaser_dkd_twix = siemens_path / 'special_cases_slaser_dkd' / 'twix' +slaser_dkd_options_twix = { + "meas_MID00053_FID16506_svs_slaser_dkd_von_wrsoff.dat": {'': 4, '_rf_off': 0, '_rf_grads_ovs_off': 0}, + "meas_MID00054_FID16507_svs_slaser_dkd_von_wrs1.dat": {'': 4, '_rf_off': 2, '_rf_grads_ovs_off': 2}, + "meas_MID00055_FID16508_svs_slaser_dkd_von_wrs2.dat": {'': 4, '_rf_off': 4, '_rf_grads_ovs_off': 4}, + "meas_MID00056_FID16509_svs_slaserVOI_dkd2_von_wrsoff.dat": {'': 4, '_vapor_ovs_rfoff': 0}, + "meas_MID00057_FID16510_svs_slaserVOI_dkd2_von_wrsw1pw3_1.dat": {'': 4, '_rf_off': 2, '_rf_grads_ovs_off': 2}, + "meas_MID00058_FID16511_svs_slaserVOI_dkd2_von_wrsw1pw3_2.dat": {'': 4, '_rf_off': 4, '_rf_grads_ovs_off': 4}, + "meas_MID00059_FID16512_svs_slaserVOI_dkd2_von_wrsw4_1.dat": {'': 4, '_vapor_ovs_rfoff': 2}, + "meas_MID00060_FID16513_svs_slaserVOI_dkd2_von_wrsw4_2.dat": {'': 4, '_vapor_ovs_rfoff': 4} +} + + +def test_twix_slaser_dkd(tmp_path): + for key in slaser_dkd_options_twix: + print(key) + subprocess.run( + ['spec2nii', 'twix', + '-e', 'image', + '-f', key.split('.')[0], + '-o', tmp_path, + '-j', data_path_slaser_dkd_twix / key]) + + for subkey in slaser_dkd_options_twix[key]: + if slaser_dkd_options_twix[key][subkey] > 0: + print(f'{key}, {subkey}, {slaser_dkd_options_twix[key][subkey]}') + img = read_nifti_mrs(tmp_path / f'{key.split(".")[0]}{subkey}.nii.gz') + assert img.shape[:4] == (1, 1, 1, 4096) + assert img.shape[-1] == slaser_dkd_options_twix[key][subkey]