ENH/svs_slaser_dkd (#136)

* Handling for slaser_dkd sequences DICOM format

* Add twix handling of slaser_dkd sequences

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

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:

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

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

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)

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]