From 82e483d352980fbddf6acfe69854fe4c18af45df Mon Sep 17 00:00:00 2001 From: wtclarke Date: Wed, 11 Oct 2023 10:56:59 +0100 Subject: [PATCH] Add handling and tests for CSI data stored in Simens enhanced DICOM --- CHANGELOG.md | 1 + spec2nii/Siemens/dicomfunctions.py | 115 ++++++++++++++++++++++++++--- tests/spec2nii_test_data | 2 +- tests/test_siemens_dicom.py | 49 ++++++++++++ 4 files changed, 156 insertions(+), 11 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a2cd509..8e75b23 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,7 @@ This document contains the Spec2nii release history in reverse chronological ord 0.7.1 (WIP) -------------------------------- - The --anon flag can be passed with any call to anonymise after writing files. +- The Siemens enhanced dicom filetype pathway now handles CSI data. 0.7.0 (Saturday 5th August 2023) -------------------------------- diff --git a/spec2nii/Siemens/dicomfunctions.py b/spec2nii/Siemens/dicomfunctions.py index f211d04..8781d33 100644 --- a/spec2nii/Siemens/dicomfunctions.py +++ b/spec2nii/Siemens/dicomfunctions.py @@ -241,6 +241,25 @@ def process_siemens_svs(img, verbose): return process_siemens_svs_xa(img, verbose) +def _get_enhanced_dcm_img_orientation_pat(img): + """Extract the image orientation patient field from an enhanced dicom scan + + :param img: DICOM image object + :type img: nibabel.nicom.dicomwrappers.SiemensWrapper + :return: 3x2 Numpy array containing the orientation info + :rtype: np.ndarray + """ + dcm_hdrs = img.dcm_data.PerFrameFunctionalGroupsSequence[0] + dcm_hdrs1 = img.dcm_data.SharedFunctionalGroupsSequence[0] + # 1) Extract dicom parameters + try: + return np.array(dcm_hdrs.PlaneOrientationSequence[0].ImageOrientationPatient).reshape(2, 3) + except AttributeError: + # For VB19 data formatted as MRSpectroscopyStorage this seems to be different to + # the XA data I've seen. But this could also be a sequence specific thing! + return np.array(dcm_hdrs1.PlaneOrientationSequence[0].ImageOrientationPatient).reshape(2, 3) + + def process_siemens_svs_xa(img, verbose): """Process Siemens DICOM SVS data acquired on a NumarisX system.""" specData = np.frombuffer(img.dcm_data[('5600', '0020')].value, dtype=np.single) @@ -250,13 +269,7 @@ def process_siemens_svs_xa(img, verbose): dcm_hdrs = img.dcm_data.PerFrameFunctionalGroupsSequence[0] dcm_hdrs1 = img.dcm_data.SharedFunctionalGroupsSequence[0] # 1) Extract dicom parameters - try: - imageOrientationPatient = np.array(dcm_hdrs.PlaneOrientationSequence[0].ImageOrientationPatient).reshape(2, 3) - except AttributeError: - # For VB19 data formatted as MRSpectroscopyStorage this seems to be different to - # the XA data I've seen. But this could also be a sequence specific thing! - imageOrientationPatient = np.array(dcm_hdrs1.PlaneOrientationSequence[0].ImageOrientationPatient).reshape(2, 3) - # VoiPosition - this does not have the FOV shift that imagePositionPatient has + imageOrientationPatient = _get_enhanced_dcm_img_orientation_pat(img) imagePositionPatient = dcm_hdrs.PlanePositionSequence[0].ImagePositionPatient # The first two elements could easily be the reverse of this. @@ -309,9 +322,55 @@ def process_siemens_csi(img, verbose): def process_siemens_csi_xa(img, verbose): - """Process Siemens DICOM CSI data acquired on a NumarisX system.""" - # NOTE: WTC expect to need to reverse the phase convention as above for svs. - raise NotImplementedError('Method process_siemens_csi_xa not implemented, example data needed!') + """Process Siemens DICOM CSI data acquired on a NumarisX system or Numaris4 Vx using enhanced DCM export""" + + warnings.warn( + 'The orientation of CSI stored in enhanced DICOM remains weakly tested.' + 'Please provide more test data.') + + specData = np.frombuffer(img.dcm_data[('5600', '0020')].value, dtype=np.single) + # It appears the phase convention is reversed in the NumarisX systems. + specDataCmplx = specData[0::2] - 1j * specData[1::2] + + # The data in img.dcm_data.SharedFunctionalGroupsSequence[0].MRSpectroscopyFOVGeometrySequence[0] + # is the uninterpolated (lower) resolution that doesn't match the data size + rows = img.dcm_data.Rows + cols = img.dcm_data.Columns + slices = int(img.dcm_data.NumberOfFrames) + spectral_points = img.dcm_data.DataPointColumns + + specDataCmplx = specDataCmplx.reshape((slices, rows, cols, spectral_points)) + specDataCmplx = np.moveaxis(specDataCmplx, (0, 1, 2), (2, 1, 0)) + + dcm_hdrs = img.dcm_data.PerFrameFunctionalGroupsSequence[0] + dcm_hdrs1 = img.dcm_data.SharedFunctionalGroupsSequence[0] + + # 1) Extract dicom parameters + imageOrientationPatient = _get_enhanced_dcm_img_orientation_pat(img) + # VoiPosition - this does not have the FOV shift that imagePositionPatient has + imagePositionPatient = dcm_hdrs.PlanePositionSequence[0].ImagePositionPatient + + # The first two elements could easily be the reverse of this. + xyzMM = np.array([float(dcm_hdrs1.PixelMeasuresSequence[0].PixelSpacing[0]), + float(dcm_hdrs1.PixelMeasuresSequence[0].PixelSpacing[1]), + float(dcm_hdrs1.PixelMeasuresSequence[0].SliceThickness)]) + + # Note that half_shift = False. This is the opposite as for classic/original DCM. + # THIS NEEDS TESTING AND DOCUMENTING. For an explanation of original shift see spec2nii/notes/siemens.md + currNiftiOrientation = dcm_to_nifti_orientation(imageOrientationPatient, + imagePositionPatient, + xyzMM, + specDataCmplx.shape[:3], + half_shift=False, + verbose=verbose) + + dwelltime = 1 / img.dcm_data.SpectralWidth + meta = extractDicomMetadata_xa(img) + + # Look for a VOI in the data + meta = _detect_and_fill_voi_enh(img, meta) + + return specDataCmplx, currNiftiOrientation, dwelltime, meta def process_siemens_csi_vx(img, verbose): @@ -388,6 +447,39 @@ def _detect_and_fill_voi(img, meta): return meta +def _detect_and_fill_voi_enh(img, meta): + """Look for a VOI volume in the headers of a CSI scan stored in enhanced dicom. + If present populate the VOI header field. + + :param img: Nibable DICOM image object + :type img: SiemensWrapper + :param meta: Existing NIfTI-MRS meta object + :type meta: hdr_ext + :return: Modified meta hdr_ext object + :rtype: hdr_ext + """ + + imageOrientationPatient = _get_enhanced_dcm_img_orientation_pat(img) + + # VoiPosition - this does not have the FOV shift that imagePositionPatient has + vol_loc_seq = img.dcm_data.VolumeLocalizationSequence + imagePositionPatient = vol_loc_seq[0].MidSlabPosition + + # This order needs validating!! + xyzMM = np.array([vol_loc_seq[1].SlabThickness, + vol_loc_seq[2].SlabThickness, + vol_loc_seq[0].SlabThickness]) + + voiNiftiOrientation = dcm_to_nifti_orientation(imageOrientationPatient, + imagePositionPatient, + xyzMM, + (1, 1, 1)) + + meta.set_standard_def('VOI', voiNiftiOrientation.Q44.tolist()) + + return meta + + def extractDicomMetadata_xa(dcmdata): """ Extract information from the nibabel DICOM object to insert into the json header ext. For Numarix X systems @@ -557,6 +649,9 @@ def set_standard_def(nifti_mrs_key, location, key, cast=None): # 'InstitutionAddress' set_standard_def('InstitutionAddress', dcmdata.dcm_data, 'InstitutionAddress') # 'TxCoil' + obj.set_standard_def( + 'TxCoil', + dcmdata.csa_header['tags']['TransmittingCoil']['items'][0]) # 'RxCoil' if len(dcmdata.csa_header['tags']['ReceivingCoil']['items']) > 0: obj.set_standard_def('RxCoil', diff --git a/tests/spec2nii_test_data b/tests/spec2nii_test_data index fd8886c..6f205d3 160000 --- a/tests/spec2nii_test_data +++ b/tests/spec2nii_test_data @@ -1 +1 @@ -Subproject commit fd8886c05f1c684c6e6c6b0769bae30d4ff1316a +Subproject commit 6f205d389085a420b38689f573045166517bb445 diff --git a/tests/test_siemens_dicom.py b/tests/test_siemens_dicom.py index 4308de9..b009a15 100644 --- a/tests/test_siemens_dicom.py +++ b/tests/test_siemens_dicom.py @@ -26,6 +26,12 @@ tim_trio_non_image = siemens_path / 'VBData' / 'tim_trio_sop_class_uid' / 'timtrio_mr_vb19a.dcm' tim_trio_mrspecstorage = siemens_path / 'VBData' / 'tim_trio_sop_class_uid' / 'timtrio_mrdcm_vb13a_2_0.dcm' +enhanced_csi_1 = siemens_path / 'enhanced_dcm_csi' / 'rk_enhanced' / '1.dcm' +enhanced_csi_2 = siemens_path / 'enhanced_dcm_csi'\ + / 'sm_enhanced' / 'Phantom.Psychiatrie^Studien^aslac.Sr27.I1.20231006084036808.dcm' +enhanced_csi_2_comp = siemens_path / 'enhanced_dcm_csi'\ + / 'sm_classic' / 'PHANTOM.MR.PSYCHIATRIE_STUDIEN_ASLAC.0027.0001.2023.10.06.08.40.21.250955.407258759.IMA' + def test_VB_svs(tmp_path): @@ -161,3 +167,46 @@ def test_sop_class_vb(tmp_path): assert img_2.shape == (1, 1, 1, 2048) assert np.iscomplexobj(img_2.dataobj) + + +def test_enhanced_csi(tmp_path): + + # First test look for right shape (16, 16, 1, 1024) + subprocess.run(['spec2nii', 'dicom', + '-f', 'csi_1', + '-o', tmp_path, + '-j', enhanced_csi_1]) + + img_1 = read_nifti_mrs(tmp_path / 'csi_1.nii.gz') + assert img_1.shape == (16, 16, 1, 1024) + assert np.iscomplexobj(img_1.dataobj) + + # Second test: compare with classic dicom export + subprocess.run(['spec2nii', 'dicom', + '-f', 'csi_2', + '-o', tmp_path, + '-j', enhanced_csi_2]) + subprocess.run(['spec2nii', 'dicom', + '-f', 'csi_2_comp', + '-o', tmp_path, + '-j', enhanced_csi_2_comp]) + + img_2 = read_nifti_mrs(tmp_path / 'csi_2.nii.gz') + img_2_comp = read_nifti_mrs(tmp_path / 'csi_2_comp.nii.gz') + + hdr_ext_codes = img_2.header.extensions.get_codes() + hdr_ext = json.loads(img_2.header.extensions[hdr_ext_codes.index(44)].get_content()) + hdr_ext_codes = img_2_comp.header.extensions.get_codes() + hdr_ext_comp = json.loads(img_2_comp.header.extensions[hdr_ext_codes.index(44)].get_content()) + + assert img_2.shape == img_2_comp.shape + assert np.allclose(img_2.dataobj[:], img_2_comp.dataobj[:]) + for key in hdr_ext: + print(key) + # Skip keys with known differences. + if key in ('RxCoil', 'SequenceName', 'PatientID', 'ConversionTime', 'OriginalFile'): + continue + if isinstance(hdr_ext[key], (float, complex)): + assert np.isclose(hdr_ext[key], hdr_ext_comp[key]) + else: + assert hdr_ext[key] == hdr_ext_comp[key]