Add handling and tests for philips mega sequence. (#131)

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

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`

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

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,38 +211,79 @@ 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:
         raise IncorrectDataDimensionsError(
             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
 

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

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'