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

@@ -722,7 +722,7 @@ def identify_integrated_references(img, inst_num):
     # Three cases as explained by DKD.
     '''
     dkd_slaserVOI_dkd
-        ScanMode = 8 
+        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

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

tests/test_siemens_special_data.py

@@ -63,7 +63,7 @@ def test_dicom_anon(tmp_path):
     assert hdr_ext['OriginalFile'][0] == str(data_path_anon)
 
 
-data_path_slaser_dkd = siemens_path / 'special_cases_slaser_dkd' / 'DICOM'
+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},
@@ -77,16 +77,46 @@ def test_dicom_anon(tmp_path):
 
 
 def test_dicom_slaser_dkd(tmp_path):
-    for idx, key in enumerate(slaser_dkd_options):
-        print(key)
+    for key in slaser_dkd_options:
         subprocess.run(
             ['spec2nii', 'dicom',
              '-f', key,
              '-o', tmp_path,
-             '-j', data_path_slaser_dkd / key])
+             '-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]