From 53eac54224b81ddf7b4a0afe99a9615bcfe09a2e Mon Sep 17 00:00:00 2001 From: constracti Date: Wed, 8 May 2019 17:34:58 +0300 Subject: [PATCH] export siemens csa header --- csa2.py | 85 ------------------------------------------- dicomtools.py | 98 ++++++++++++++++++++++++++++++++++++++++++++++++-- nifti2dicom.py | 22 ++++++++++-- 3 files changed, 114 insertions(+), 91 deletions(-) delete mode 100644 csa2.py diff --git a/csa2.py b/csa2.py deleted file mode 100644 index e4cf627..0000000 --- a/csa2.py +++ /dev/null @@ -1,85 +0,0 @@ -# http://nipy.org/nibabel/dicom/siemens_csa.html - -import math - -def decode(arr): - assert arr[0:4] == b"SV10" - # arr[4:8] # unused - ntags = int.from_bytes(arr[8:12], "little") - # int.from_bytes(arr[8:12], "little") # unused 77 - i = 16 - hdr = {} - for ctags in range(0, ntags): - key = arr[i:i+64].split(bytes(1), 1)[0].decode() - hdr[key] = {} - hdr[key]["VM"] = int.from_bytes(arr[i+64:i+68], "little") - hdr[key]["VR"] = arr[i+68:i+72].split(bytes(1), 1)[0].decode() - hdr[key]["SyngoDT"] = int.from_bytes(arr[i+72:i+76], "little") - nitems = int.from_bytes(arr[i+76:i+80], "little") - # int.from_bytes(arr[80:84], "little") # unused 77 or 205 - i += 84 - hdr[key]["Data"] = [] - for citems in range(0, nitems): - item_len = int.from_bytes(arr[i:i+4], "little") - int.from_bytes(arr[i+4:i+8], "little") # unused - int.from_bytes(arr[i+8:i+12], "little") # unused - int.from_bytes(arr[i+12:i+16], "little") # unused - if item_len: - hdr[key]["Data"].append(arr[i+16:i+16+item_len].split(bytes(1), 1)[0].decode()) - else: - hdr[key]["Data"].append(None) - i += 16 + math.ceil(item_len / 4) * 4 - return hdr - -def encode(hdr): - arr = b"" - arr += b"SV10" - arr += bytes(4) # unused - arr += len(hdr).to_bytes(4, "little") - arr += bytes(4) # unused 77 - for key, val in hdr.items(): - arr += key.encode().ljust(64, bytes(1)) - arr += val["VM"].to_bytes(4, "little") - arr += val["VR"].encode().ljust(4, bytes(1)) - arr += val["SyngoDT"].to_bytes(4, "little") - arr += len(val["Data"]).to_bytes(4, "little") - arr += bytes(4) # unused 77 or 205 - for dat in val["Data"]: - item_len = len(dat) - arr += item_len.to_bytes(4, "little") - arr += bytes(4) # unused - arr += bytes(4) # unused - arr += bytes(4) # unused - if item_len: - arr += dat.encode().ljust(math.ceil(item_len / 4) * 4, bytes(1)) - arr += bytes(4) # append one more zero - return arr - -def diff(hdr1, hdr2): - keys1 = iter(sorted(hdr1.keys())) - keys2 = iter(sorted(hdr2.keys())) - key1 = next(keys1, None) - key2 = next(keys2, None) - while key1 and key2: - val1 = hdr1[key1] - val2 = hdr2[key2] - if key1 < key2: - print("< {}: {}".format(key1, val1)) - key1 = next(keys1, None) - elif key2 < key1: - print("> {}: {}".format(key2, val2)) - key2 = next(keys2, None) - else: - if val1 != val2: - print("< {}: {}".format(key1, val1)) - print("> {}: {}".format(key2, val2)) - key1 = next(keys1, None) - key2 = next(keys2, None) - while key1: - val1 = hdr1[key1] - print("< {}: {}".format(key1, val1)) - key1 = next(keys1, None) - while key2: - val2 = hdr2[key2] - print("> {}: {}".format(key2, val2)) - key2 = next(keys2, None) diff --git a/dicomtools.py b/dicomtools.py index 75e1231..01fa52f 100644 --- a/dicomtools.py +++ b/dicomtools.py @@ -2,6 +2,7 @@ import re import math import datetime +from collections import OrderedDict import numpy import pydicom @@ -10,8 +11,6 @@ except ImportError: pass -import csa2 - def file_is_valid(filename): try: pydicom.dcmread(filename, specific_tags=[]) @@ -48,7 +47,7 @@ def get_affine(dataset1, dataset2): # number of slices nslices = 1 if (0x0029, 0x0010) in dataset1 and dataset1[0x0029, 0x0010].value == "SIEMENS CSA HEADER": - csa_image_header_info = csa2.decode(dataset1[0x0029, 0x1010].value) + csa_image_header_info = csa2_decode(dataset1[0x0029, 0x1010].value) if csa_image_header_info["NumberOfImagesInMosaic"]["Data"]: nslices = int(csa_image_header_info["NumberOfImagesInMosaic"]["Data"][0]) # shape, zooms & affine @@ -89,6 +88,99 @@ def get_affine(dataset1, dataset2): # return return shape, zooms, affine + +######## +# csa2 # +######## + +# http://nipy.org/nibabel/dicom/siemens_csa.html + +def csa2_decode(arr): + assert arr[0:4] == b"SV10" + # arr[4:8] # unused + ntags = int.from_bytes(arr[8:12], "little") + # int.from_bytes(arr[12:16], "little") # unused 77 + i = 16 + hdr = OrderedDict() + for ctags in range(0, ntags): + key = arr[i:i+64].split(bytes(1), 1)[0].decode() + hdr[key] = {} + hdr[key]["VM"] = int.from_bytes(arr[i+64:i+68], "little") + hdr[key]["VR"] = arr[i+68:i+72].split(bytes(1), 1)[0].decode() + hdr[key]["SyngoDT"] = int.from_bytes(arr[i+72:i+76], "little") + nitems = int.from_bytes(arr[i+76:i+80], "little") + # int.from_bytes(arr[80:84], "little") # unused 77 or 205 + i += 84 + hdr[key]["Data"] = [] + for citems in range(nitems): + item_len = int.from_bytes(arr[i:i+4], "little") + # int.from_bytes(arr[i+4:i+8], "little") # unused + # int.from_bytes(arr[i+8:i+12], "little") # unused + # int.from_bytes(arr[i+12:i+16], "little") # unused + if item_len: + hdr[key]["Data"].append(arr[i+16:i+16+item_len].split(bytes(1), 1)[0].decode()) + else: + hdr[key]["Data"].append(None) + i += 16 + math.ceil(item_len / 4) * 4 + return hdr + +def csa2_encode(hdr): + arr = b"" + arr += b"SV10" + arr += bytes(4) # unused + arr += len(hdr).to_bytes(4, "little") + arr += bytes(4) # unused 77 + for key, val in hdr.items(): + arr += key.encode().ljust(64, bytes(1)) + arr += val["VM"].to_bytes(4, "little") + arr += val["VR"].encode().ljust(4, bytes(1)) + arr += val["SyngoDT"].to_bytes(4, "little") + arr += len(val["Data"]).to_bytes(4, "little") + arr += bytes(4) # unused 77 or 205 + for dat in val["Data"]: + item_len = len(dat) + 1 if dat is not None else 0 + arr += item_len.to_bytes(4, "little") + arr += bytes(4) # unused + arr += bytes(4) # unused + arr += bytes(4) # unused + if item_len: + arr += dat.encode().ljust(math.ceil(item_len / 4) * 4, bytes(1)) + arr += bytes(4) # append one more zero + return arr + +def csa2_diff(hdr1, hdr2): + keys1 = iter(sorted(hdr1.keys())) + keys2 = iter(sorted(hdr2.keys())) + key1 = next(keys1, None) + key2 = next(keys2, None) + while key1 and key2: + val1 = hdr1[key1] + val2 = hdr2[key2] + if key1 < key2: + print("< {}: {}".format(key1, val1)) + key1 = next(keys1, None) + elif key2 < key1: + print("> {}: {}".format(key2, val2)) + key2 = next(keys2, None) + else: + if val1 != val2: + print("< {}: {}".format(key1, val1)) + print("> {}: {}".format(key2, val2)) + key1 = next(keys1, None) + key2 = next(keys2, None) + while key1: + val1 = hdr1[key1] + print("< {}: {}".format(key1, val1)) + key1 = next(keys1, None) + while key2: + val2 = hdr2[key2] + print("> {}: {}".format(key2, val2)) + key2 = next(keys2, None) + +########## +# linear # +########## + def _prop(dataset, tag): if type(tag) is str: return dataset.data_element(tag) diff --git a/nifti2dicom.py b/nifti2dicom.py index fbd450a..14f8773 100755 --- a/nifti2dicom.py +++ b/nifti2dicom.py @@ -43,8 +43,6 @@ def nifti2dicom(path): if (0x0019, 0x0010) in dataset and dataset[0x0019, 0x0010].value == "SIEMENS MR HEADER": # (0x0019, 0x100b) SliceMeasurementDuration dataset.pop((0x0019, 0x100b), None) - # (0x0019, 0x1016) TimeAfterStart - dataset.pop((0x0019, 0x1016), None) # (0x0019, 0x1029) MosaicRefAcqTimes dataset.pop((0x0019, 0x1029), None) # Window Explanation Algorithm not specified @@ -56,7 +54,7 @@ def nifti2dicom(path): dataset.pop((0x0028, 0x1055), None) if (0x0029, 0x0010) in dataset and dataset[0x0029, 0x0010].value == "SIEMENS CSA HEADER": # (0x0029, 0x1010) CSA Image Header Info - dataset.pop((0x0029, 0x1010), None) + csa_image_header_info = dicomtools.csa2_decode(dataset[0x0029, 0x1010].value) if (0x0043, 0x0010) in dataset and dataset[0x0043, 0x0010].value == "GEMS_PARM_01": # (0x0043, 0x1028) Unique image iden dataset.pop((0x0043, 0x1028), None) @@ -117,6 +115,9 @@ def nifti2dicom(path): # (0x0019, 0x1015) SlicePosition_PCS if (0x0019, 0x1015) in dataset: dicomtools.linear_float_array((0x0019, 0x1015), dataset, dataset1, dataset2) + # (0x0019, 0x1016) TimeAfterStart + if (0x0019, 0x1016) in dataset: + dicomtools.linear_float((0x0019, 0x1016), dataset, dataset1, dataset2) elif (0x0019, 0x0010) in dataset and dataset[0x0019, 0x0010].value in ["GEMS_ACQU_01", "GEMS_IDEN_01"]: # (0x0019, 0x10a2) Raw data run number if (0x0019, 0x10a2) in dataset: @@ -139,6 +140,21 @@ def nifti2dicom(path): # (0x0028, 0x0107) Largest Image Pixel Value if (0x0028, 0x0107) in dataset: dataset[0x0028, 0x0107].value = data_slice.max() + if (0x0029, 0x0010) in dataset and dataset[0x0029, 0x0010].value == "SIEMENS CSA HEADER": + # (0x0029, 0x1010) CSA Image Header Info + if csa_image_header_info["Actual3DImaPartNumber"]["Data"]: + csa_image_header_info["Actual3DImaPartNumber"]["Data"][0] = str(f).ljust(8) + elif not csa_image_header_info["MosaicRefAcqTimes"]["Data"]: # TODO linear int on csa_image_header_info + csa_image_header_info["ProtocolSliceNumber"]["Data"][0] = str(f).ljust(8) + # csa_image_header_info["GSWDDataType"] CORONAL + # csa_image_header_info["RFSWDDataType"] CORONAL + # csa_image_header_info["ICE_Dims"]["Data"][0] * + # csa_image_header_info["MosaicRefAcqTimes"]["Data"] FMRI + # csa_image_header_info["SliceMeasurementDuration"]["Data"][0] CORONAL + csa_image_header_info["SlicePosition_PCS"]["Data"][0:3] = ["{:.8f}".format(x) for x in dataset.ImagePositionPatient] + if csa_image_header_info["TimeAfterStart"]["Data"]: + csa_image_header_info["TimeAfterStart"]["Data"][0] = "{:.8f}".format(dataset[0x0019, 0x1016].value) + dataset[0x0029, 0x1010].value = dicomtools.csa2_encode(csa_image_header_info) if (0x2001, 0x0010) in dataset and dataset[0x2001, 0x0010].value == "Philips Imaging DD 001": # (0x2001, 0x100a) Slice Number MR if (0x2001, 0x100a) in dataset: