Skip to content
Snippets Groups Projects
Commit d716673d authored by Tamino Huxohl's avatar Tamino Huxohl
Browse files

formatting

parent 058b2484
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 152 additions and 58 deletions
mu_map/data/prepare.py
+ 152
− 58
View file @ d716673d
...@@ -51,7 +51,6 @@ headers.file_recon_nac_nsc = "file_recon_nac_nsc" ...@@ -51,7 +51,6 @@ headers.file_recon_nac_nsc = "file_recon_nac_nsc"
headers.file_mu_map = "file_mu_map" headers.file_mu_map = "file_mu_map"
def get_protocol(projection: pydicom.dataset.FileDataset) -> str: def get_protocol(projection: pydicom.dataset.FileDataset) -> str:
""" """
Get the protocol (stress, rest) of a projection image by checking if Get the protocol (stress, rest) of a projection image by checking if
...@@ -69,7 +68,9 @@ def get_protocol(projection: pydicom.dataset.FileDataset) -> str: ...@@ -69,7 +68,9 @@ def get_protocol(projection: pydicom.dataset.FileDataset) -> str:
raise ValueError(f"Unkown protocol in projection {projection.SeriesDescription}") raise ValueError(f"Unkown protocol in projection {projection.SeriesDescription}")
def find_projections(dicom_images: List[pydicom.dataset.FileDataset]) -> pydicom.dataset.FileDataset: def find_projections(
dicom_images: List[pydicom.dataset.FileDataset],
) -> pydicom.dataset.FileDataset:
""" """
Find all projections in a list of DICOM images belonging to a study. Find all projections in a list of DICOM images belonging to a study.
...@@ -81,7 +82,9 @@ def find_projections(dicom_images: List[pydicom.dataset.FileDataset]) -> pydicom ...@@ -81,7 +82,9 @@ def find_projections(dicom_images: List[pydicom.dataset.FileDataset]) -> pydicom
for dicom_image in _filter: for dicom_image in _filter:
# filter for allows series descriptions # filter for allows series descriptions
if dicom_image.SeriesDescription not in ["Stress", "Rest", "Stress_2"]: if dicom_image.SeriesDescription not in ["Stress", "Rest", "Stress_2"]:
logger.warning(f"Skip projection with unkown protocol [{dicom_image.SeriesDescription}]") logger.warning(
f"Skip projection with unkown protocol [{dicom_image.SeriesDescription}]"
)
# print(f" - {DICOMTime.Study.to_datetime(di)}, {DICOMTime.Series.to_datetime(di)}, {DICOMTime.Content.to_datetime(di)}, {DICOMTime.Acquisition.to_datetime(di)}") # print(f" - {DICOMTime.Study.to_datetime(di)}, {DICOMTime.Series.to_datetime(di)}, {DICOMTime.Content.to_datetime(di)}, {DICOMTime.Acquisition.to_datetime(di)}")
continue continue
dicom_images.append(dicom_image) dicom_images.append(dicom_image)
...@@ -91,7 +94,10 @@ def find_projections(dicom_images: List[pydicom.dataset.FileDataset]) -> pydicom ...@@ -91,7 +94,10 @@ def find_projections(dicom_images: List[pydicom.dataset.FileDataset]) -> pydicom
return dicom_images return dicom_images
def is_recon_type(scatter_corrected: bool, attenuation_corrected: bool) -> Callable[pydicom.dataset.FileDataset, bool]:
def is_recon_type(
scatter_corrected: bool, attenuation_corrected: bool
) -> Callable[[pydicom.dataset.FileDataset], bool]:
""" """
Get a filter function for reconstructions that are (non-)scatter and/or (non-)attenuation corrected. Get a filter function for reconstructions that are (non-)scatter and/or (non-)attenuation corrected.
...@@ -111,7 +117,12 @@ def is_recon_type(scatter_corrected: bool, attenuation_corrected: bool) -> Calla ...@@ -111,7 +117,12 @@ def is_recon_type(scatter_corrected: bool, attenuation_corrected: bool) -> Calla
return lambda dicom: filter_str in dicom.SeriesDescription return lambda dicom: filter_str in dicom.SeriesDescription
def find_reconstruction(dicom_images: List[pydicom.dataset.FileDataset], projection: pydicom.dataset.FileDataset, scatter_corrected: bool, attenuation_corrected: bool) -> List[pydicom.dataset.FileDataset]: def find_reconstruction(
dicom_images: List[pydicom.dataset.FileDataset],
projection: pydicom.dataset.FileDataset,
scatter_corrected: bool,
attenuation_corrected: bool,
) -> List[pydicom.dataset.FileDataset]:
""" """
Find a reconstruction in a list of dicom images of a study belonging to a projection. Find a reconstruction in a list of dicom images of a study belonging to a projection.
...@@ -125,35 +136,52 @@ def find_reconstruction(dicom_images: List[pydicom.dataset.FileDataset], project ...@@ -125,35 +136,52 @@ def find_reconstruction(dicom_images: List[pydicom.dataset.FileDataset], project
_filter = filter(lambda image: "RECON TOMO" in image.ImageType, dicom_images) _filter = filter(lambda image: "RECON TOMO" in image.ImageType, dicom_images)
_filter = filter(lambda image: protocol in image.SeriesDescription, _filter) _filter = filter(lambda image: protocol in image.SeriesDescription, _filter)
_filter = filter(lambda image: STUDY_DESCRIPTION in image.SeriesDescription, _filter) _filter = filter(
_filter = filter(lambda image: "CT" not in image.SeriesDescription, _filter) # remove µ-maps lambda image: STUDY_DESCRIPTION in image.SeriesDescription, _filter
)
_filter = filter(
lambda image: "CT" not in image.SeriesDescription, _filter
) # remove µ-maps
_filter = filter(is_recon_type(scatter_corrected, attenuation_corrected), _filter) _filter = filter(is_recon_type(scatter_corrected, attenuation_corrected), _filter)
# _filter = list(_filter) # _filter = list(_filter)
# print("DEBUG Reconstructions: ") # print("DEBUG Reconstructions: ")
# for r in _filter: # for r in _filter:
# try: # try:
# print(f" - {r.SeriesDescription:>50} at {DICOMTime.Study.to_datetime(r)}, {DICOMTime.Series.to_datetime(r)}, {DICOMTime.Content.to_datetime(r)}, {DICOMTime.Acquisition.to_datetime(r)}") # print(f" - {r.SeriesDescription:>50} at {DICOMTime.Study.to_datetime(r)}, {DICOMTime.Series.to_datetime(r)}, {DICOMTime.Content.to_datetime(r)}, {DICOMTime.Acquisition.to_datetime(r)}")
# except Exception as e: # except Exception as e:
# print(f"Error {e}") # print(f"Error {e}")
_filter = filter(lambda image: DICOMTime.Acquisition.to_datetime(image) == DICOMTime.Acquisition.to_datetime(projection), _filter) _filter = filter(
lambda image: DICOMTime.Acquisition.to_datetime(image)
== DICOMTime.Acquisition.to_datetime(projection),
_filter,
)
dicom_images = list(_filter) dicom_images = list(_filter)
if len(dicom_images) == 0: if len(dicom_images) == 0:
raise ValueError(f"No reconstruction with SC={scatter_corrected}, AC={attenuation_corrected} available") raise ValueError(
f"No reconstruction with SC={scatter_corrected}, AC={attenuation_corrected} available"
)
# sort oldest to be first # sort oldest to be first
if len(dicom_images) > 1: if len(dicom_images) > 1:
logger.warning(f"Multiple reconstructions ({len(dicom_images)}) with SC={scatter_corrected}, AC={attenuation_corrected} for projection {projection.SeriesDescription} of patient {projection.PatientID}") logger.warning(
dicom_images.sort(key=lambda image: DICOMTime.Series.to_datetime(image), reverse=True) f"Multiple reconstructions ({len(dicom_images)}) with SC={scatter_corrected}, AC={attenuation_corrected} for projection {projection.SeriesDescription} of patient {projection.PatientID}"
)
dicom_images.sort(
key=lambda image: DICOMTime.Series.to_datetime(image), reverse=True
)
return dicom_images[0] return dicom_images[0]
MAX_TIME_DIFF_S = 30 MAX_TIME_DIFF_S = 30
def find_attenuation_map( def find_attenuation_map(
dicom_images: List[pydicom.dataset.FileDataset], projection: pydicom.dataset.FileDataset, reconstructions: List[pydicom.dataset.FileDataset], dicom_images: List[pydicom.dataset.FileDataset],
projection: pydicom.dataset.FileDataset,
reconstructions: List[pydicom.dataset.FileDataset],
) -> pydicom.dataset.FileDataset: ) -> pydicom.dataset.FileDataset:
""" """
Find a reconstruction in a list of dicom images of a study belonging to a projection and reconstructions. Find a reconstruction in a list of dicom images of a study belonging to a projection and reconstructions.
...@@ -164,13 +192,28 @@ def find_attenuation_map( ...@@ -164,13 +192,28 @@ def find_attenuation_map(
:returns: the according attenuation map :returns: the according attenuation map
""" """
protocol = get_protocol(projection) protocol = get_protocol(projection)
recon_times = list(map(lambda recon: DICOMTime.Series.to_datetime(recon), reconstructions)) recon_times = list(
map(lambda recon: DICOMTime.Series.to_datetime(recon), reconstructions)
)
_filter = filter(lambda image: "RECON TOMO" in image.ImageType, dicom_images) _filter = filter(lambda image: "RECON TOMO" in image.ImageType, dicom_images)
_filter = filter(lambda image: protocol in image.SeriesDescription, _filter) _filter = filter(lambda image: protocol in image.SeriesDescription, _filter)
_filter = filter(lambda image: STUDY_DESCRIPTION in image.SeriesDescription, _filter) _filter = filter(
lambda image: STUDY_DESCRIPTION in image.SeriesDescription, _filter
)
_filter = filter(lambda image: " µ-map]" in image.SeriesDescription, _filter) _filter = filter(lambda image: " µ-map]" in image.SeriesDescription, _filter)
_filter = filter(lambda image: any(map(lambda recon_time: (DICOMTime.Series.to_datetime(image) - recon_time).seconds < MAX_TIME_DIFF_S, recon_times)), _filter) _filter = filter(
lambda image: any(
map(
lambda recon_time: (
DICOMTime.Series.to_datetime(image) - recon_time
).seconds
< MAX_TIME_DIFF_S,
recon_times,
)
),
_filter,
)
dicom_images = list(_filter) dicom_images = list(_filter)
if len(dicom_images) == 0: if len(dicom_images) == 0:
...@@ -178,13 +221,19 @@ def find_attenuation_map( ...@@ -178,13 +221,19 @@ def find_attenuation_map(
# sort oldest to be first # sort oldest to be first
if len(dicom_images) > 1: if len(dicom_images) > 1:
logger.warning(f"Multiple attenuation maps ({len(dicom_images)}) for projection {projection.SeriesDescription} of patient {projection.PatientID}") logger.warning(
dicom_images.sort(key=lambda image: DICOMTime.Series.to_datetime(image), reverse=True) f"Multiple attenuation maps ({len(dicom_images)}) for projection {projection.SeriesDescription} of patient {projection.PatientID}"
)
dicom_images.sort(
key=lambda image: DICOMTime.Series.to_datetime(image), reverse=True
)
return dicom_images[0] return dicom_images[0]
def get_relevant_images(patient: pydicom.dataset.FileDataset, dicom_dir: str) -> List[pydicom.dataset.FileDataset]: def get_relevant_images(
patient: pydicom.dataset.FileDataset, dicom_dir: str
) -> List[pydicom.dataset.FileDataset]:
""" """
Get all relevant images of a patient. Get all relevant images of a patient.
...@@ -243,7 +292,6 @@ def get_relevant_images(patient: pydicom.dataset.FileDataset, dicom_dir: str) -> ...@@ -243,7 +292,6 @@ def get_relevant_images(patient: pydicom.dataset.FileDataset, dicom_dir: str) ->
return dicom_images return dicom_images
if __name__ == "__main__": if __name__ == "__main__":
parser = argparse.ArgumentParser( parser = argparse.ArgumentParser(
description="Prepare a dataset from DICOM directories", description="Prepare a dataset from DICOM directories",
...@@ -315,38 +363,62 @@ if __name__ == "__main__": ...@@ -315,38 +363,62 @@ if __name__ == "__main__":
data = pd.DataFrame(dict([(key, []) for key in vars(headers).keys()])) data = pd.DataFrame(dict([(key, []) for key in vars(headers).keys()]))
for i, patient in enumerate(patients, start=1): for i, patient in enumerate(patients, start=1):
logger.debug(f"Process patient {str(i):>3}/{len(patients)} - {patient.PatientName.given_name}, {patient.PatientName.family_name}:") logger.debug(
f"Process patient {str(i):>3}/{len(patients)} - {patient.PatientName.given_name}, {patient.PatientName.family_name}:"
)
dicom_images = get_relevant_images(patient, dicom_dir_by_patient[patient.PatientID]) dicom_images = get_relevant_images(
patient, dicom_dir_by_patient[patient.PatientID]
)
projections = find_projections(dicom_images) projections = find_projections(dicom_images)
logger.debug(f"- Found {len(projections)}: projections with protocols {list(map(lambda p: get_protocol(p), projections))}") logger.debug(
f"- Found {len(projections)}: projections with protocols {list(map(lambda p: get_protocol(p), projections))}"
)
for projection in projections: for projection in projections:
protocol = get_protocol(projection) protocol = get_protocol(projection)
reconstructions = [] reconstructions = []
recon_headers = [headers.file_recon_nac_nsc, headers.file_recon_nac_sc, headers.file_recon_ac_nsc, headers.file_recon_ac_sc] recon_headers = [
recon_postfixes = ["recon_nac_nsc", "recon_nac_sc", "recon_ac_nsc", "recon_ac_sc"] headers.file_recon_nac_nsc,
for ac, sc in [(False, False), (False, True), (True, False), (True, True)]: headers.file_recon_nac_sc,
recon = find_reconstruction(dicom_images, projection, scatter_corrected=sc, attenuation_corrected=ac) headers.file_recon_ac_nsc,
headers.file_recon_ac_sc,
]
recon_postfixes = [
"recon_nac_nsc",
"recon_nac_sc",
"recon_ac_nsc",
"recon_ac_sc",
]
for ac, sc in [
(False, False),
(False, True),
(True, False),
(True, True),
]:
recon = find_reconstruction(
dicom_images,
projection,
scatter_corrected=sc,
attenuation_corrected=ac,
)
reconstructions.append(recon) reconstructions.append(recon)
mu_map = find_attenuation_map(dicom_images, projection, reconstructions) mu_map = find_attenuation_map(dicom_images, projection, reconstructions)
# for i, projection in enumerate(projections):
# for i, projection in enumerate(projections):
# _time = DICOMTime.Series.to_datetime(projection) # _time = DICOMTime.Series.to_datetime(projection)
# print(f" - Projection: {projection.SeriesDescription:>10} at {DICOMTime.Study.to_datetime(projection)}, {DICOMTime.Series.to_datetime(projection)}, {DICOMTime.Content.to_datetime(projection)}, {DICOMTime.Acquisition.to_datetime(projection)}") # print(f" - Projection: {projection.SeriesDescription:>10} at {DICOMTime.Study.to_datetime(projection)}, {DICOMTime.Series.to_datetime(projection)}, {DICOMTime.Content.to_datetime(projection)}, {DICOMTime.Acquisition.to_datetime(projection)}")
# recons = [] # recons = []
# for sc, ac in [(False, False), (False, True), (True, False), (True, True)]: # for sc, ac in [(False, False), (False, True), (True, False), (True, True)]:
# r = find_reconstruction(dicom_images, projection, scatter_corrected=sc, attenuation_corrected=ac) # r = find_reconstruction(dicom_images, projection, scatter_corrected=sc, attenuation_corrected=ac)
# print(f" - {r.SeriesDescription:>50} at {DICOMTime.Study.to_datetime(r)}, {DICOMTime.Series.to_datetime(r)}, {DICOMTime.Content.to_datetime(r)}, {DICOMTime.Acquisition.to_datetime(r)}") # print(f" - {r.SeriesDescription:>50} at {DICOMTime.Study.to_datetime(r)}, {DICOMTime.Series.to_datetime(r)}, {DICOMTime.Content.to_datetime(r)}, {DICOMTime.Acquisition.to_datetime(r)}")
# recons.append(r) # recons.append(r)
# mu_map = find_attenuation_map(dicom_images, projection, recons) # mu_map = find_attenuation_map(dicom_images, projection, recons)
# print(f" - {mu_map.SeriesDescription:>50} at {DICOMTime.Study.to_datetime(mu_map)}, {DICOMTime.Series.to_datetime(mu_map)}, {DICOMTime.Content.to_datetime(mu_map)}") # print(f" - {mu_map.SeriesDescription:>50} at {DICOMTime.Study.to_datetime(mu_map)}, {DICOMTime.Series.to_datetime(mu_map)}, {DICOMTime.Content.to_datetime(mu_map)}")
# print(f" -") # print(f" -")
# extract pixel spacings and assert that they are equal for all reconstruction images # extract pixel spacings and assert that they are equal for all reconstruction images
_map_lists = map( _map_lists = map(
lambda image: [*image.PixelSpacing, image.SliceThickness], lambda image: [*image.PixelSpacing, image.SliceThickness],
...@@ -382,18 +454,14 @@ if __name__ == "__main__": ...@@ -382,18 +454,14 @@ if __name__ == "__main__":
], ],
[*reconstructions, mu_map], [*reconstructions, mu_map],
) )
_map_lists = map( _map_lists = map(lambda shape: list(map(int, shape)), _map_lists)
lambda shape: list(map(int, shape)), _map_lists
)
_map_ndarrays = map(lambda shape: np.array(shape), _map_lists) _map_ndarrays = map(lambda shape: np.array(shape), _map_lists)
shapes = list(_map_ndarrays) shapes = list(_map_ndarrays)
shapes.sort(key=lambda shape: shape[2]) shapes.sort(key=lambda shape: shape[2])
shape = shapes[0] shape = shapes[0]
# extract and sort energy windows # extract and sort energy windows
energy_windows = ( energy_windows = projection.EnergyWindowInformationSequence
projection.EnergyWindowInformationSequence
)
energy_windows = map( energy_windows = map(
lambda ew: ew.EnergyWindowRangeSequence[0], energy_windows lambda ew: ew.EnergyWindowRangeSequence[0], energy_windows
) )
...@@ -414,40 +482,66 @@ if __name__ == "__main__": ...@@ -414,40 +482,66 @@ if __name__ == "__main__":
headers.weight: float(projection.PatientWeight), headers.weight: float(projection.PatientWeight),
headers.size: float(projection.PatientSize), headers.size: float(projection.PatientSize),
headers.protocol: protocol, headers.protocol: protocol,
headers.datetime_acquisition: DICOMTime.Series.to_datetime(projection), headers.datetime_acquisition: DICOMTime.Series.to_datetime(
headers.datetime_reconstruction: DICOMTime.Series.to_datetime(reconstructions[0]), projection
),
headers.datetime_reconstruction: DICOMTime.Series.to_datetime(
reconstructions[0]
),
headers.pixel_spacing_x: pixel_spacing[0], headers.pixel_spacing_x: pixel_spacing[0],
headers.pixel_spacing_y: pixel_spacing[1], headers.pixel_spacing_y: pixel_spacing[1],
headers.pixel_spacing_z: pixel_spacing[2], headers.pixel_spacing_z: pixel_spacing[2],
headers.shape_x: shape[0], headers.shape_x: shape[0],
headers.shape_y: shape[1], headers.shape_y: shape[1],
headers.shape_z: shape[2], headers.shape_z: shape[2],
headers.radiopharmaceutical: projection.RadiopharmaceuticalInformationSequence[0].Radiopharmaceutical, headers.radiopharmaceutical: projection.RadiopharmaceuticalInformationSequence[
headers.radionuclide_dose: projection.RadiopharmaceuticalInformationSequence[0].RadionuclideTotalDose, 0
headers.radionuclide_code: projection.RadiopharmaceuticalInformationSequence[0].RadionuclideCodeSequence[0].CodeValue, ].Radiopharmaceutical,
headers.radionuclide_meaning: projection.RadiopharmaceuticalInformationSequence[0].RadionuclideCodeSequence[0].CodeMeaning, headers.radionuclide_dose: projection.RadiopharmaceuticalInformationSequence[
0
].RadionuclideTotalDose,
headers.radionuclide_code: projection.RadiopharmaceuticalInformationSequence[
0
]
.RadionuclideCodeSequence[0]
.CodeValue,
headers.radionuclide_meaning: projection.RadiopharmaceuticalInformationSequence[
0
]
.RadionuclideCodeSequence[0]
.CodeMeaning,
headers.energy_window_peak_lower: energy_windows[0][0], headers.energy_window_peak_lower: energy_windows[0][0],
headers.energy_window_peak_upper: energy_windows[0][1], headers.energy_window_peak_upper: energy_windows[0][1],
headers.energy_window_scatter_lower: energy_windows[1][0], headers.energy_window_scatter_lower: energy_windows[1][0],
headers.energy_window_scatter_upper: energy_windows[1][1], headers.energy_window_scatter_upper: energy_windows[1][1],
headers.detector_count: len(projection.DetectorInformationSequence), headers.detector_count: len(projection.DetectorInformationSequence),
headers.collimator_type: projection.DetectorInformationSequence[0].CollimatorType, headers.collimator_type: projection.DetectorInformationSequence[
headers.rotation_start: float(projection.RotationInformationSequence[0].StartAngle), 0
headers.rotation_step: float(projection.RotationInformationSequence[0].AngularStep), ].CollimatorType,
headers.rotation_scan_arc: float(projection.RotationInformationSequence[0].ScanArc), headers.rotation_start: float(
projection.RotationInformationSequence[0].StartAngle
),
headers.rotation_step: float(
projection.RotationInformationSequence[0].AngularStep
),
headers.rotation_scan_arc: float(
projection.RotationInformationSequence[0].ScanArc
),
} }
_filename_base = f"{_id:04d}-{protocol.lower()}" _filename_base = f"{_id:04d}-{protocol.lower()}"
_ext = "dcm" _ext = "dcm"
img_headers = [headers.file_projection, *recon_headers, headers.file_mu_map] img_headers = [
headers.file_projection,
*recon_headers,
headers.file_mu_map,
]
img_postfixes = ["projection", *recon_postfixes, "mu_map"] img_postfixes = ["projection", *recon_postfixes, "mu_map"]
images = [projection, *reconstructions, mu_map] images = [projection, *reconstructions, mu_map]
for header, image, postfix in zip(img_headers, images, img_postfixes): for header, image, postfix in zip(img_headers, images, img_postfixes):
_image = pydicom.dcmread(image.filename) _image = pydicom.dcmread(image.filename)
filename = f"{_filename_base}-{postfix}.{_ext}" filename = f"{_filename_base}-{postfix}.{_ext}"
pydicom.dcmwrite( pydicom.dcmwrite(os.path.join(args.images_dir, filename), _image)
os.path.join(args.images_dir, filename), _image
)
row[header] = filename row[header] = filename
_id += 1 _id += 1
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment