From f7e34a9ff008fcf920751977812b578265c80b04 Mon Sep 17 00:00:00 2001 From: pierocor Date: Thu, 26 Nov 2020 20:35:06 +0100 Subject: [PATCH 1/4] BUG: Load projection data and references from mayo dataset --- .../datasets/ct/examples/mayo_reconstruct.py | 64 ++++-- odl/contrib/datasets/ct/mayo.py | 206 +++++++++++------- odl/contrib/datasets/ct/mayo_dicom_dict.py | 6 +- 3 files changed, 171 insertions(+), 105 deletions(-) diff --git a/odl/contrib/datasets/ct/examples/mayo_reconstruct.py b/odl/contrib/datasets/ct/examples/mayo_reconstruct.py index 81f211d47a7..966cf9dc94d 100644 --- a/odl/contrib/datasets/ct/examples/mayo_reconstruct.py +++ b/odl/contrib/datasets/ct/examples/mayo_reconstruct.py @@ -4,26 +4,35 @@ stored in the location indicated by "mayo_dir". In this example we only use a subset of the data for performance reasons, -there are ~48 000 projections in the full dataset. +there are ~32 000 projections in the full dataset. """ - +import numpy as np import odl from odl.contrib.datasets.ct import mayo +from time import perf_counter -mayo_dir = '' # replace with your local folder +# define data folders +proj_folder = odl.__path__[0] + '/../../data/LDCT-and-Projection-data/' \ + 'L004/08-21-2018-10971/1.000000-Full dose projections-24362/' +img_folder = odl.__path__[0] + '/../../data/LDCT-and-Projection-data/' \ + 'L004/08-21-2018-84608/1.000000-Full dose images-59704/' -# Load reference reconstruction -volume_folder = mayo_dir + '/Training Cases/L067/full_1mm_sharp' -partition, volume = mayo.load_reconstruction(volume_folder) +# Load projection data +print("Loading projection data from {:s}".format(proj_folder)) +geometry, proj_data = mayo.load_projections(proj_folder, + indices=slice(16000, 19000)) +# Load reconstruction data +print("Loading reference data from {:s}".format(img_folder)) +recon_space, volume = mayo.load_reconstruction(img_folder) -# Load a subset of the projection data -data_folder = mayo_dir + '/Training Cases/L067/full_DICOM-CT-PD' -geometry, proj_data = mayo.load_projections(data_folder, - indices=slice(20000, 28000)) +# ray transform +ray_trafo = odl.tomo.RayTransform(recon_space, geometry) -# Reconstruction space and ray transform -space = odl.uniform_discr_frompartition(partition, dtype='float32') -ray_trafo = odl.tomo.RayTransform(space, geometry) +# Interpolate projection data for a flat grid +radial_dist = geometry.src_radius + geometry.det_radius +flat_proj_data = mayo.interpolate_flat_grid(proj_data, + ray_trafo.range.grid, + radial_dist) # Define FBP operator fbp = odl.tomo.fbp_op(ray_trafo, padding=True) @@ -32,16 +41,27 @@ td_window = odl.tomo.tam_danielson_window(ray_trafo, n_pi=3) # Calculate FBP reconstruction -fbp_result = fbp(td_window * proj_data) +start = perf_counter() +fbp_result = fbp(td_window * flat_proj_data) +stop = perf_counter() +print('FBP done after {:.3f} seconds'.format(stop-start)) + +fbp_result_HU = (fbp_result-0.0192)/0.0192*1000 + +# Save reconstruction in Numpy format +fbp_filename = proj_folder+'fbp_result.npy' +print("Saving reconstruction data in {:s}".format(fbp_filename)) +np.save(fbp_filename, fbp_result_HU) + +# Compare the computed recon to reference reconstruction (coronal slice) +ref = recon_space.element(volume) +diff = recon_space.element(volume - fbp_result_HU.asarray()) # Compare the computed recon to reference reconstruction (coronal slice) -ref = space.element(volume) -fbp_result.show('Recon (coronal)', clim=[0.7, 1.3]) -ref.show('Reference (coronal)', clim=[0.7, 1.3]) -(ref - fbp_result).show('Diff (coronal)', clim=[-0.1, 0.1]) +fbp_result_HU.show('Recon (coronal)') +ref.show('Reference (coronal)') +diff.show('Diff (coronal)') -# Also visualize sagittal slice (note that we only used a subset) coords = [0, None, None] -fbp_result.show('Recon (sagittal)', clim=[0.7, 1.3], coords=coords) -ref.show('Reference (sagittal)', clim=[0.7, 1.3], coords=coords) -(ref - fbp_result).show('Diff (sagittal)', clim=[-0.1, 0.1], coords=coords) +fbp_result_HU.show('Recon (sagittal)', coords=coords) +ref.show('Reference (sagittal)', coords=coords) \ No newline at end of file diff --git a/odl/contrib/datasets/ct/mayo.py b/odl/contrib/datasets/ct/mayo.py index d74a487f913..4932603fda1 100644 --- a/odl/contrib/datasets/ct/mayo.py +++ b/odl/contrib/datasets/ct/mayo.py @@ -19,17 +19,19 @@ from __future__ import division import numpy as np import os -import dicom +import pydicom import odl import tqdm -from dicom.datadict import DicomDictionary, NameDict, CleanName +from pydicom.datadict import DicomDictionary, keyword_dict from odl.discr.discr_utils import linear_interpolator from odl.contrib.datasets.ct.mayo_dicom_dict import new_dict_items # Update the DICOM dictionary with the extra Mayo tags DicomDictionary.update(new_dict_items) -NameDict.update((CleanName(tag), tag) for tag in new_dict_items) +# Update the reverse mapping from name to tag +new_names_dict = dict([(val[4], tag) for tag, val in new_dict_items.items()]) +keyword_dict.update(new_names_dict) __all__ = ('load_projections', 'load_reconstruction') @@ -52,35 +54,42 @@ def _read_projections(folder, indices): for i, file_name in enumerate(tqdm.tqdm(file_names, 'Loading projection data')): # read the file - dataset = dicom.read_file(folder + '/' + file_name) + dataset = pydicom.read_file(folder + '/' + file_name) + + if data_array is None: + # Get some required data + rows = dataset.NumberofDetectorRows + cols = dataset.NumberofDetectorColumns + rescale_intercept = dataset.RescaleIntercept + rescale_slope = dataset.RescaleSlope + + data_array = np.empty((len(file_names), cols, rows), + dtype='float32') + angles = np.empty(len(file_names), dtype='float32') + + else: + # Sanity checks + assert rows == dataset.NumberofDetectorRows + assert cols == dataset.NumberofDetectorColumns + assert rescale_intercept == dataset.RescaleIntercept + assert rescale_slope == dataset.RescaleSlope - # Get some required data - rows = dataset.NumberofDetectorRows - cols = dataset.NumberofDetectorColumns - hu_factor = dataset.HUCalibrationFactor - rescale_intercept = dataset.RescaleIntercept - rescale_slope = dataset.RescaleSlope # Load the array as bytes proj_array = np.array(np.frombuffer(dataset.PixelData, 'H'), dtype='float32') - proj_array = proj_array.reshape([rows, cols], order='F').T + proj_array = proj_array.reshape([cols, rows]) + # proj_array = proj_array.reshape([rows, cols], order='F').T - # Rescale array + # Rescale array (no HU) proj_array *= rescale_slope proj_array += rescale_intercept - proj_array /= hu_factor - - # Store results - if data_array is None: - # We need to load the first dataset before we know the shape - data_array = np.empty((len(file_names), cols, rows), - dtype='float32') data_array[i] = proj_array[:, ::-1] + angles[i] = dataset.DetectorFocalCenterAngularPosition datasets.append(dataset) - return datasets, data_array + return datasets, data_array, angles def load_projections(folder, indices=None): @@ -93,6 +102,10 @@ def load_projections(folder, indices=None): indices : optional Indices of the projections to load. Accepts advanced indexing such as slice or list of indices. + num_slices: int, optional + Number of slices to consider for the reconstruction volume; + the other parameters are hard-coded. With default value (None) + a *temporary* volume is created. Returns ------- @@ -102,21 +115,20 @@ def load_projections(folder, indices=None): Projection data, given as the line integral of the linear attenuation coefficient (g/cm^3). Its unit is thus g/cm^2. """ - datasets, data_array = _read_projections(folder, indices) + datasets, data_array, angles = _read_projections(folder, indices) - # Get the angles - angles = [d.DetectorFocalCenterAngularPosition for d in datasets] - angles = -np.unwrap(angles) - np.pi # different definition of angles + # Reverse angular axis and set origin at 6 o'clock + angles = -np.unwrap(angles) - np.pi # Set minimum and maximum corners - shape = np.array([datasets[0].NumberofDetectorColumns, - datasets[0].NumberofDetectorRows]) - pixel_size = np.array([datasets[0].DetectorElementTransverseSpacing, - datasets[0].DetectorElementAxialSpacing]) + det_shape = np.array([datasets[0].NumberofDetectorColumns, + datasets[0].NumberofDetectorRows]) + det_pixel_size = np.array([datasets[0].DetectorElementTransverseSpacing, + datasets[0].DetectorElementAxialSpacing]) # Correct from center of pixel to corner of pixel - minp = -(np.array(datasets[0].DetectorCentralElement) - 0.5) * pixel_size - maxp = minp + shape * pixel_size + det_minp = -(np.array(datasets[0].DetectorCentralElement) - 0.5) * det_pixel_size + det_maxp = det_minp + det_shape * det_pixel_size # Select geometry parameters src_radius = datasets[0].DetectorFocalCenterRadialDistance @@ -126,46 +138,48 @@ def load_projections(folder, indices=None): # For unknown reasons, mayo does not include the tag # "TableFeedPerRotation", which is what we want. # Instead we manually compute the pitch - pitch = ((datasets[-1].DetectorFocalCenterAxialPosition - - datasets[0].DetectorFocalCenterAxialPosition) / - ((np.max(angles) - np.min(angles)) / (2 * np.pi))) - - # Get flying focal spot data - offset_axial = np.array([d.SourceAxialPositionShift for d in datasets]) - offset_angular = np.array([d.SourceAngularPositionShift for d in datasets]) - offset_radial = np.array([d.SourceRadialDistanceShift for d in datasets]) - - # TODO(adler-j): Implement proper handling of flying focal spot. - # Currently we do not fully account for it, merely making some "first - # order corrections" to the detector position and radial offset. - - # Update angles with flying focal spot (in plane direction). - # This increases the resolution of the reconstructions. - angles = angles - offset_angular - - # We correct for the mean offset due to the rotated angles, we need to - # shift the detector. - offset_detector_by_angles = det_radius * np.mean(offset_angular) - minp[0] -= offset_detector_by_angles - maxp[0] -= offset_detector_by_angles - - # We currently apply only the mean of the offsets - src_radius = src_radius + np.mean(offset_radial) - - # Partially compensate for a movement of the source by moving the object - # instead. We need to rescale by the magnification to get the correct - # change in the detector. This approximation is only exactly valid on the - # axis of rotation. - mean_offset_along_axis_for_ffz = np.mean(offset_axial) * ( - src_radius / (src_radius + det_radius)) + table_dist = datasets[-1].DetectorFocalCenterAxialPosition - \ + datasets[0].DetectorFocalCenterAxialPosition + num_rot = (angles[-1] - angles[0]) / (2 * np.pi) + pitch = table_dist / num_rot + + # TODO: Understand and re-implement the flying focal spot + # # Get flying focal spot data + # offset_axial = np.array([d.SourceAxialPositionShift for d in datasets]) + # offset_angular = np.array([d.SourceAngularPositionShift for d in datasets]) + # offset_radial = np.array([d.SourceRadialDistanceShift for d in datasets]) + + # # TODO(adler-j): Implement proper handling of flying focal spot. + # # Currently we do not fully account for it, merely making some "first + # # order corrections" to the detector position and radial offset. + + # # Update angles with flying focal spot (in plane direction). + # # This increases the resolution of the reconstructions. + # angles = angles - offset_angular + + # # We correct for the mean offset due to the rotated angles, we need to + # # shift the detector. + # offset_detector_by_angles = det_radius * np.mean(offset_angular) + # det_minp[0] -= offset_detector_by_angles + # det_maxp[0] -= offset_detector_by_angles + + # # We currently apply only the mean of the offsets + # src_radius = src_radius + np.mean(offset_radial) + + # # Partially compensate for a movement of the source by moving the object + # # instead. We need to rescale by the magnification to get the correct + # # change in the detector. This approximation is only exactly valid on the + # # axis of rotation. + # mean_offset_along_axis_for_ffz = np.mean(offset_axial) * ( + # src_radius / (src_radius + det_radius)) # Create partition for detector - detector_partition = odl.uniform_partition(minp, maxp, shape) + detector_partition = odl.uniform_partition(det_minp, det_maxp, det_shape) # Convert offset to odl definitions - offset_along_axis = (mean_offset_along_axis_for_ffz + - datasets[0].DetectorFocalCenterAxialPosition - - angles[0] / (2 * np.pi) * pitch) + # offset_along_axis = (mean_offset_along_axis_for_ffz + + offset_along_axis = datasets[0].DetectorFocalCenterAxialPosition - \ + angles[0] / (2 * np.pi) * pitch # Assemble geometry angle_partition = odl.nonuniform_partition(angles) @@ -176,24 +190,52 @@ def load_projections(folder, indices=None): pitch=pitch, offset_along_axis=offset_along_axis) - # Create a *temporary* ray transform (we need its range) - spc = odl.uniform_discr([-1] * 3, [1] * 3, [32] * 3) - ray_trafo = odl.tomo.RayTransform(spc, geometry, interp='linear') + return geometry, data_array + + +def get_default_recon_space(): + # Create a *temporary* ray transform (we need its range) + num_slices = 97 + pixel_spacing = np.array([0.75,0.75]) + num_pixel = np.array([512,512]) + slice_dist = 5. + origin = np.zeros(3) + mid_table = (datasets[0].DetectorFocalCenterAxialPosition + + datasets[-1].DetectorFocalCenterAxialPosition) / 2 + min_pt = np.copy(origin) + min_pt[:2] -= pixel_spacing * num_pixel / 2 + min_pt[2] += mid_table - num_slices * slice_dist / 2 + + max_pt = np.copy(min_pt) + max_pt[:2] += pixel_spacing * num_pixel + max_pt[2] += num_slices * slice_dist + + recon_dim = np.array([*num_pixel, num_slices], dtype=np.int32) + recon_space = odl.uniform_discr(min_pt, max_pt, + shape=recon_dim, + dtype=np.float32) + return recon_space + + +# ray_trafo = odl.tomo.RayTransform(recon_space, geometry, interp='linear') + +def interpolate_flat_grid(data_array, range_grid, radial_dist): # convert coordinates - theta, up, vp = ray_trafo.range.grid.meshgrid - d = src_radius + det_radius - u = d * np.arctan(up / d) - v = d / np.sqrt(d**2 + up**2) * vp + theta, up, vp = range_grid.meshgrid #ray_trafo.range.grid.meshgrid + # d = src_radius + det_radius + u = radial_dist * np.arctan(up / radial_dist) + v = radial_dist / np.sqrt(radial_dist**2 + up**2) * vp # Calculate projection data in rectangular coordinates since we have no # backend that supports cylindrical interpolator = linear_interpolator( - data_array, ray_trafo.range.coord_vectors + data_array, range_grid.coord_vectors # ray_trafo.range.grid.coord_vectors ) proj_data = interpolator((theta, u, v)) - return geometry, proj_data.asarray() + return proj_data + def load_reconstruction(folder, slice_start=0, slice_end=-1): @@ -227,7 +269,7 @@ def load_reconstruction(folder, slice_start=0, slice_end=-1): This function should handle all of these peculiarities and give a volume with the correct coordinate system attached. """ - file_names = sorted([f for f in os.listdir(folder) if f.endswith(".IMA")]) + file_names = sorted([f for f in os.listdir(folder) if f.endswith(".dcm")]) if len(file_names) == 0: raise ValueError('No DICOM files found in {}'.format(folder)) @@ -239,7 +281,7 @@ def load_reconstruction(folder, slice_start=0, slice_end=-1): for file_name in tqdm.tqdm(file_names, 'loading volume data'): # read the file - dataset = dicom.read_file(folder + '/' + file_name) + dataset = pydicom.read_file(folder + '/' + file_name) # Get parameters pixel_size = np.array(dataset.PixelSpacing) @@ -257,18 +299,17 @@ def load_reconstruction(folder, slice_start=0, slice_end=-1): # TODO: Optimize these computations hu_values = (dataset.RescaleSlope * data_array + dataset.RescaleIntercept) - densities = (hu_values + 1000) / 1000 # Store results - volumes.append(densities) + volumes.append(hu_values) datasets.append(dataset) voxel_size = np.array(list(pixel_size) + [pixel_thickness]) shape = np.array([rows, cols, len(volumes)]) # Compute geometry parameters - mid_pt = (np.array(dataset.ReconstructionTargetCenterPatient) - - np.array(dataset.DataCollectionCenterPatient)) + mid_pt = np.array(dataset.ReconstructionTargetCenterPatient) + mid_pt[1] += datasets[0].TableHeight reconstruction_size = (voxel_size * shape) min_pt = mid_pt - reconstruction_size / 2 max_pt = mid_pt + reconstruction_size / 2 @@ -294,10 +335,11 @@ def load_reconstruction(folder, slice_start=0, slice_end=-1): max_pt[2] += 0.5 * slice_distance partition = odl.uniform_partition(min_pt, max_pt, shape) + recon_space = odl.uniform_discr_frompartition(partition, dtype='float32') volume = np.transpose(np.array(volumes), (1, 2, 0)) - return partition, volume + return recon_space, volume if __name__ == '__main__': diff --git a/odl/contrib/datasets/ct/mayo_dicom_dict.py b/odl/contrib/datasets/ct/mayo_dicom_dict.py index dabd45dd9ee..2b9f6a3e119 100644 --- a/odl/contrib/datasets/ct/mayo_dicom_dict.py +++ b/odl/contrib/datasets/ct/mayo_dicom_dict.py @@ -1,5 +1,8 @@ # _dicom_dict.py -"""DICOM data dictionary auto-generated by DICOM-CT-PD-dict_to_pydicom.py""" +""" +DICOM data dictionary auto-generated by DICOM-CT-PD-dict_to_pydicom.py +Update: Added (7041,1001) "Water Attenuation Coefficient" +""" from __future__ import absolute_import new_dict_items = { @@ -3932,6 +3935,7 @@ 0x70391008: ('CS', '1', 'Scatter Correction Flag', '', 'ScatterCorrectionFlag'), 0x70391009: ('CS', '1', 'Log Flag', '', 'LogFlag'), 0x70410010: ('LO', '1', 'Lesion Information Module', '', 'LesionInformationModule'), + 0x70411001: ('DS', '1', 'Water Attenuation Coefficient', '','WaterAttenuationCoefficient'), 0x70411003: ('US', '1', 'Numberof Lesions', '', 'NumberofLesions'), 0x70411004: ('ST', '1', 'Lesion Pathology Array', '', 'LesionPathologyArray'), 0x70411005: ('FL', '1-n', 'Lesion Angular Position Array', '', 'LesionAngularPositionArray'), From 8a5a19734525e5e7fe4d8fc788e1199db02da815 Mon Sep 17 00:00:00 2001 From: Piero Date: Fri, 18 Dec 2020 13:08:16 +0000 Subject: [PATCH 2/4] FIX: Add flying focal spot --- .../datasets/ct/examples/mayo_reconstruct.py | 8 +- odl/contrib/datasets/ct/mayo.py | 79 +++++++++---------- odl/contrib/datasets/ct/mayo_dicom_dict.py | 4 +- 3 files changed, 43 insertions(+), 48 deletions(-) diff --git a/odl/contrib/datasets/ct/examples/mayo_reconstruct.py b/odl/contrib/datasets/ct/examples/mayo_reconstruct.py index 966cf9dc94d..0b42372c9ec 100644 --- a/odl/contrib/datasets/ct/examples/mayo_reconstruct.py +++ b/odl/contrib/datasets/ct/examples/mayo_reconstruct.py @@ -49,7 +49,7 @@ fbp_result_HU = (fbp_result-0.0192)/0.0192*1000 # Save reconstruction in Numpy format -fbp_filename = proj_folder+'fbp_result.npy' +fbp_filename = proj_folder+'/fbp_result.npy' print("Saving reconstruction data in {:s}".format(fbp_filename)) np.save(fbp_filename, fbp_result_HU) @@ -58,9 +58,9 @@ diff = recon_space.element(volume - fbp_result_HU.asarray()) # Compare the computed recon to reference reconstruction (coronal slice) -fbp_result_HU.show('Recon (coronal)') -ref.show('Reference (coronal)') -diff.show('Diff (coronal)') +fbp_result_HU.show('Recon (axial)') +ref.show('Reference (axial)') +diff.show('Diff (axial)') coords = [0, None, None] fbp_result_HU.show('Recon (sagittal)', coords=coords) diff --git a/odl/contrib/datasets/ct/mayo.py b/odl/contrib/datasets/ct/mayo.py index 4932603fda1..14eb32da398 100644 --- a/odl/contrib/datasets/ct/mayo.py +++ b/odl/contrib/datasets/ct/mayo.py @@ -19,6 +19,7 @@ from __future__ import division import numpy as np import os +import sys import pydicom import odl import tqdm @@ -40,6 +41,8 @@ def _read_projections(folder, indices): """Read mayo projections from a folder.""" datasets = [] + data_array = [] + angles = [] # Get the relevant file names file_names = sorted([f for f in os.listdir(folder) if f.endswith(".dcm")]) @@ -49,24 +52,24 @@ def _read_projections(folder, indices): file_names = file_names[indices] - data_array = None - for i, file_name in enumerate(tqdm.tqdm(file_names, 'Loading projection data')): # read the file - dataset = pydicom.read_file(folder + '/' + file_name) - - if data_array is None: + try: + dataset = pydicom.read_file(folder + os.path.sep + file_name) + except: + print("corrupted file: {}".format(file_name), file=sys.stderr) + print("error:\n{}".format(sys.exc_info()[1]), file=sys.stderr) + print("skipping to next file..", file=sys.stderr) + continue + + if not data_array: # Get some required data rows = dataset.NumberofDetectorRows cols = dataset.NumberofDetectorColumns rescale_intercept = dataset.RescaleIntercept rescale_slope = dataset.RescaleSlope - data_array = np.empty((len(file_names), cols, rows), - dtype='float32') - angles = np.empty(len(file_names), dtype='float32') - else: # Sanity checks assert rows == dataset.NumberofDetectorRows @@ -85,14 +88,16 @@ def _read_projections(folder, indices): proj_array *= rescale_slope proj_array += rescale_intercept - data_array[i] = proj_array[:, ::-1] - angles[i] = dataset.DetectorFocalCenterAngularPosition + data_array.append(proj_array[:, ::-1]) + angles.append(dataset.DetectorFocalCenterAngularPosition) datasets.append(dataset) + data_array = np.array(data_array) + angles = np.array(angles) return datasets, data_array, angles -def load_projections(folder, indices=None): +def load_projections(folder, indices=None, use_ffs=True): """Load geometry and data stored in Mayo format from folder. Parameters @@ -143,35 +148,17 @@ def load_projections(folder, indices=None): num_rot = (angles[-1] - angles[0]) / (2 * np.pi) pitch = table_dist / num_rot - # TODO: Understand and re-implement the flying focal spot - # # Get flying focal spot data - # offset_axial = np.array([d.SourceAxialPositionShift for d in datasets]) - # offset_angular = np.array([d.SourceAngularPositionShift for d in datasets]) - # offset_radial = np.array([d.SourceRadialDistanceShift for d in datasets]) - - # # TODO(adler-j): Implement proper handling of flying focal spot. - # # Currently we do not fully account for it, merely making some "first - # # order corrections" to the detector position and radial offset. - - # # Update angles with flying focal spot (in plane direction). - # # This increases the resolution of the reconstructions. - # angles = angles - offset_angular - - # # We correct for the mean offset due to the rotated angles, we need to - # # shift the detector. - # offset_detector_by_angles = det_radius * np.mean(offset_angular) - # det_minp[0] -= offset_detector_by_angles - # det_maxp[0] -= offset_detector_by_angles - - # # We currently apply only the mean of the offsets - # src_radius = src_radius + np.mean(offset_radial) + # offsets: focal spot -> detector’s focal center + offset_angular = np.array([d.SourceAngularPositionShift for d in datasets]) + offset_radial = np.array([d.SourceRadialDistanceShift for d in datasets]) + offset_axial = np.array([d.SourceAxialPositionShift for d in datasets]) + + # angles have inverse convention + shift_d = np.cos(-offset_angular) * (src_radius + offset_radial) - src_radius + shift_t = np.sin(-offset_angular) * (src_radius + offset_radial) + shift_r = + offset_axial - # # Partially compensate for a movement of the source by moving the object - # # instead. We need to rescale by the magnification to get the correct - # # change in the detector. This approximation is only exactly valid on the - # # axis of rotation. - # mean_offset_along_axis_for_ffz = np.mean(offset_axial) * ( - # src_radius / (src_radius + det_radius)) + shifts = np.transpose(np.vstack([shift_d, shift_t, shift_r])) # Create partition for detector detector_partition = odl.uniform_partition(det_minp, det_maxp, det_shape) @@ -183,12 +170,20 @@ def load_projections(folder, indices=None): # Assemble geometry angle_partition = odl.nonuniform_partition(angles) + + # Flying focal spot + src_shift_func = None + if use_ffs: + src_shift_func = lambda angle: odl.tomo.flying_focal_spot( + angle, apart=angle_partition, shifts=shifts) + geometry = odl.tomo.ConeBeamGeometry(angle_partition, detector_partition, src_radius=src_radius, det_radius=det_radius, pitch=pitch, - offset_along_axis=offset_along_axis) + offset_along_axis=offset_along_axis, + src_shift_func=src_shift_func) return geometry, data_array @@ -281,7 +276,7 @@ def load_reconstruction(folder, slice_start=0, slice_end=-1): for file_name in tqdm.tqdm(file_names, 'loading volume data'): # read the file - dataset = pydicom.read_file(folder + '/' + file_name) + dataset = pydicom.read_file(folder + os.path.sep + file_name) # Get parameters pixel_size = np.array(dataset.PixelSpacing) diff --git a/odl/contrib/datasets/ct/mayo_dicom_dict.py b/odl/contrib/datasets/ct/mayo_dicom_dict.py index 2b9f6a3e119..a2c6256c510 100644 --- a/odl/contrib/datasets/ct/mayo_dicom_dict.py +++ b/odl/contrib/datasets/ct/mayo_dicom_dict.py @@ -3916,8 +3916,8 @@ 0x70311031: ('FL', '1', 'Constant Radial Distance', '', 'ConstantRadialDistance'), 0x70311033: ('FL', '1-n', 'Detector Central Element', '', 'DetectorCentralElement'), 0x70330010: ('LO', '1', 'Source Dynamics Module', '', 'SourceDynamicsModule'), - 0x7033100B: ('FL', '1', 'Source Axial Position Shift', '', 'SourceAxialPositionShift'), - 0x7033100C: ('FL', '1', 'Source Angular Position Shift', '', 'SourceAngularPositionShift'), + 0x7033100B: ('FL', '1', 'Source Angular Position Shift', '', 'SourceAngularPositionShift'), + 0x7033100C: ('FL', '1', 'Source Axial Position Shift', '', 'SourceAxialPositionShift'), 0x7033100D: ('FL', '1', 'Source Radial Distance Shift', '', 'SourceRadialDistanceShift'), 0x7033100E: ('CS', '1', 'Flying Focal Spot Mode', '', 'FlyingFocalSpotMode'), 0x70331013: ('US', '1', 'Numberof Source Angular Steps', '', 'NumberofSourceAngularSteps'), From a85173887ea16345a5a6cc6b6763bfcb3fbf4375 Mon Sep 17 00:00:00 2001 From: Piero Date: Fri, 18 Dec 2020 14:16:33 +0000 Subject: [PATCH 3/4] STY: Add documentation --- .../datasets/ct/examples/mayo_reconstruct.py | 10 +-- odl/contrib/datasets/ct/mayo.py | 78 ++++++++----------- odl/contrib/datasets/ct/mayo_dicom_dict.py | 3 +- 3 files changed, 38 insertions(+), 53 deletions(-) diff --git a/odl/contrib/datasets/ct/examples/mayo_reconstruct.py b/odl/contrib/datasets/ct/examples/mayo_reconstruct.py index 0b42372c9ec..efc47058f43 100644 --- a/odl/contrib/datasets/ct/examples/mayo_reconstruct.py +++ b/odl/contrib/datasets/ct/examples/mayo_reconstruct.py @@ -1,10 +1,10 @@ """Reconstruct Mayo dataset using FBP and compare to reference recon. Note that this example requires that Mayo has been previously downloaded and is -stored in the location indicated by "mayo_dir". +stored in the location indicated by "proj_folder" and "rec_folder". In this example we only use a subset of the data for performance reasons, -there are ~32 000 projections in the full dataset. +there are ~32 000 projections per patient in the full dataset. """ import numpy as np import odl @@ -14,7 +14,7 @@ # define data folders proj_folder = odl.__path__[0] + '/../../data/LDCT-and-Projection-data/' \ 'L004/08-21-2018-10971/1.000000-Full dose projections-24362/' -img_folder = odl.__path__[0] + '/../../data/LDCT-and-Projection-data/' \ +rec_folder = odl.__path__[0] + '/../../data/LDCT-and-Projection-data/' \ 'L004/08-21-2018-84608/1.000000-Full dose images-59704/' # Load projection data @@ -22,8 +22,8 @@ geometry, proj_data = mayo.load_projections(proj_folder, indices=slice(16000, 19000)) # Load reconstruction data -print("Loading reference data from {:s}".format(img_folder)) -recon_space, volume = mayo.load_reconstruction(img_folder) +print("Loading reference data from {:s}".format(rec_folder)) +recon_space, volume = mayo.load_reconstruction(rec_folder) # ray transform ray_trafo = odl.tomo.RayTransform(recon_space, geometry) diff --git a/odl/contrib/datasets/ct/mayo.py b/odl/contrib/datasets/ct/mayo.py index 14eb32da398..4577964143c 100644 --- a/odl/contrib/datasets/ct/mayo.py +++ b/odl/contrib/datasets/ct/mayo.py @@ -42,7 +42,6 @@ def _read_projections(folder, indices): """Read mayo projections from a folder.""" datasets = [] data_array = [] - angles = [] # Get the relevant file names file_names = sorted([f for f in os.listdir(folder) if f.endswith(".dcm")]) @@ -82,19 +81,15 @@ def _read_projections(folder, indices): proj_array = np.array(np.frombuffer(dataset.PixelData, 'H'), dtype='float32') proj_array = proj_array.reshape([cols, rows]) - # proj_array = proj_array.reshape([rows, cols], order='F').T - - # Rescale array (no HU) - proj_array *= rescale_slope - proj_array += rescale_intercept - data_array.append(proj_array[:, ::-1]) - angles.append(dataset.DetectorFocalCenterAngularPosition) datasets.append(dataset) data_array = np.array(data_array) - angles = np.array(angles) - return datasets, data_array, angles + # Rescale array + data_array *= rescale_slope + data_array += rescale_intercept + + return datasets, data_array def load_projections(folder, indices=None, use_ffs=True): @@ -107,10 +102,9 @@ def load_projections(folder, indices=None, use_ffs=True): indices : optional Indices of the projections to load. Accepts advanced indexing such as slice or list of indices. - num_slices: int, optional - Number of slices to consider for the reconstruction volume; - the other parameters are hard-coded. With default value (None) - a *temporary* volume is created. + use_ffs : bool, optional + If ``True``, a source shift is applied to compensate the flying focal spot. + Default: ``True`` Returns ------- @@ -120,10 +114,12 @@ def load_projections(folder, indices=None, use_ffs=True): Projection data, given as the line integral of the linear attenuation coefficient (g/cm^3). Its unit is thus g/cm^2. """ - datasets, data_array, angles = _read_projections(folder, indices) + datasets, data_array = _read_projections(folder, indices) + # Get the angles + angles = np.array([d.DetectorFocalCenterAngularPosition for d in datasets]) # Reverse angular axis and set origin at 6 o'clock - angles = -np.unwrap(angles) - np.pi + angles = -np.unwrap(angles) - np.pi # Set minimum and maximum corners det_shape = np.array([datasets[0].NumberofDetectorColumns, @@ -148,11 +144,11 @@ def load_projections(folder, indices=None, use_ffs=True): num_rot = (angles[-1] - angles[0]) / (2 * np.pi) pitch = table_dist / num_rot - # offsets: focal spot -> detector’s focal center + # offsets: detector’s focal center -> focal spot offset_angular = np.array([d.SourceAngularPositionShift for d in datasets]) offset_radial = np.array([d.SourceRadialDistanceShift for d in datasets]) offset_axial = np.array([d.SourceAxialPositionShift for d in datasets]) - + # angles have inverse convention shift_d = np.cos(-offset_angular) * (src_radius + offset_radial) - src_radius shift_t = np.sin(-offset_angular) * (src_radius + offset_radial) @@ -164,7 +160,6 @@ def load_projections(folder, indices=None, use_ffs=True): detector_partition = odl.uniform_partition(det_minp, det_maxp, det_shape) # Convert offset to odl definitions - # offset_along_axis = (mean_offset_along_axis_for_ffz + offset_along_axis = datasets[0].DetectorFocalCenterAxialPosition - \ angles[0] / (2 * np.pi) * pitch @@ -187,35 +182,26 @@ def load_projections(folder, indices=None, use_ffs=True): return geometry, data_array - -def get_default_recon_space(): - # Create a *temporary* ray transform (we need its range) - num_slices = 97 - pixel_spacing = np.array([0.75,0.75]) - num_pixel = np.array([512,512]) - slice_dist = 5. - origin = np.zeros(3) - mid_table = (datasets[0].DetectorFocalCenterAxialPosition + - datasets[-1].DetectorFocalCenterAxialPosition) / 2 - min_pt = np.copy(origin) - min_pt[:2] -= pixel_spacing * num_pixel / 2 - min_pt[2] += mid_table - num_slices * slice_dist / 2 - - max_pt = np.copy(min_pt) - max_pt[:2] += pixel_spacing * num_pixel - max_pt[2] += num_slices * slice_dist - - recon_dim = np.array([*num_pixel, num_slices], dtype=np.int32) - recon_space = odl.uniform_discr(min_pt, max_pt, - shape=recon_dim, - dtype=np.float32) - return recon_space +def interpolate_flat_grid(data_array, range_grid, radial_dist): + """Return the linear interpolator of the projection data on a flat detector. + Parameters + ---------- + data_array : `numpy.ndarray` + Projection data on the cylindrical detector that should be interpolated. + range_grid : RectGrid + Rectilinear grid on the flat detector. + radial_dist : float + The constant radial distance, that is the distance between the detector’s + focal center and its central element. -# ray_trafo = odl.tomo.RayTransform(recon_space, geometry, interp='linear') + Returns + ------- + proj_data : `numpy.ndarray` + Interpolated projection data on the flat rectilinear grid. + """ -def interpolate_flat_grid(data_array, range_grid, radial_dist): # convert coordinates theta, up, vp = range_grid.meshgrid #ray_trafo.range.grid.meshgrid # d = src_radius + det_radius @@ -225,14 +211,13 @@ def interpolate_flat_grid(data_array, range_grid, radial_dist): # Calculate projection data in rectangular coordinates since we have no # backend that supports cylindrical interpolator = linear_interpolator( - data_array, range_grid.coord_vectors # ray_trafo.range.grid.coord_vectors + data_array, range_grid.coord_vectors ) proj_data = interpolator((theta, u, v)) return proj_data - def load_reconstruction(folder, slice_start=0, slice_end=-1): """Load a volume from folder, also returns the corresponding partition. @@ -291,7 +276,6 @@ def load_reconstruction(folder, slice_start=0, slice_end=-1): data_array = np.rot90(data_array, -1) # Convert from storage type to densities - # TODO: Optimize these computations hu_values = (dataset.RescaleSlope * data_array + dataset.RescaleIntercept) diff --git a/odl/contrib/datasets/ct/mayo_dicom_dict.py b/odl/contrib/datasets/ct/mayo_dicom_dict.py index a2c6256c510..7341021fda8 100644 --- a/odl/contrib/datasets/ct/mayo_dicom_dict.py +++ b/odl/contrib/datasets/ct/mayo_dicom_dict.py @@ -1,7 +1,8 @@ # _dicom_dict.py """ DICOM data dictionary auto-generated by DICOM-CT-PD-dict_to_pydicom.py -Update: Added (7041,1001) "Water Attenuation Coefficient" +Update: Add (7041,1001) "Water Attenuation Coefficient" + Swap (7033,100B) and (7033,100C) """ from __future__ import absolute_import From 96a7d5e9b09b295e56e99f06b2e0edc69a45a0c4 Mon Sep 17 00:00:00 2001 From: pierocor Date: Tue, 2 Feb 2021 16:11:08 +0100 Subject: [PATCH 4/4] STY: pull request mods --- odl/contrib/datasets/README.md | 2 +- .../datasets/ct/examples/mayo_reconstruct.py | 44 +++++++------ odl/contrib/datasets/ct/mayo.py | 65 ++++++++++--------- 3 files changed, 58 insertions(+), 53 deletions(-) diff --git a/odl/contrib/datasets/README.md b/odl/contrib/datasets/README.md index 822084e9255..1c622815c9e 100644 --- a/odl/contrib/datasets/README.md +++ b/odl/contrib/datasets/README.md @@ -22,7 +22,7 @@ Reference datasets with accompanying ODL geometries etc. * `walnut_data` * `lotus_root_data` - CT data as provided by Mayo Clinic. The data is from a human and of high resolution (512x512). To access the data, see [the webpage](https://www.aapm.org/GrandChallenge/LowDoseCT/#registration). Note that downloading this dataset requires signing up and signing a terms of use form. + CT data as provided by Mayo Clinic. The data is from a human and of high resolution (512x512). To access the data, see [the webpage](https://ctcicblog.mayo.edu/hubcap/patient-ct-projection-data-library/). Note that downloading this dataset requires installing the NBIA Data Retriever. * `load_projections` * `load_reconstruction` * `images` diff --git a/odl/contrib/datasets/ct/examples/mayo_reconstruct.py b/odl/contrib/datasets/ct/examples/mayo_reconstruct.py index efc47058f43..ab5a4ac4c2d 100644 --- a/odl/contrib/datasets/ct/examples/mayo_reconstruct.py +++ b/odl/contrib/datasets/ct/examples/mayo_reconstruct.py @@ -1,29 +1,38 @@ """Reconstruct Mayo dataset using FBP and compare to reference recon. -Note that this example requires that Mayo has been previously downloaded and is -stored in the location indicated by "proj_folder" and "rec_folder". +Note that this example requires that projection and reconstruction data +of a CT scan from the Mayo dataset have been previously downloaded (see +[the webpage](https://ctcicblog.mayo.edu/hubcap/patient-ct-projection-data-library/)) +and are stored in the locations indicated by "proj_dir" and "rec_dir". -In this example we only use a subset of the data for performance reasons, -there are ~32 000 projections per patient in the full dataset. + +In this example we only use a subset of the data for performance reasons. +The number of projections per patient varies in the full dataset. To get +a reconstruction of the central part of the volume, modify the indices in +the argument of the `mayo.load_projections` function accordingly. """ import numpy as np +import os import odl from odl.contrib.datasets.ct import mayo from time import perf_counter -# define data folders -proj_folder = odl.__path__[0] + '/../../data/LDCT-and-Projection-data/' \ - 'L004/08-21-2018-10971/1.000000-Full dose projections-24362/' -rec_folder = odl.__path__[0] + '/../../data/LDCT-and-Projection-data/' \ - 'L004/08-21-2018-84608/1.000000-Full dose images-59704/' +# replace with your local directory +mayo_dir = '' +# define projection and reconstruction data directories +# e.g. for patient L004 full dose CT scan: +proj_dir = os.path.join( + mayo_dir, 'L004/08-21-2018-10971/1.000000-Full dose projections-24362/') +rec_dir = os.path.join( + mayo_dir, 'L004/08-21-2018-84608/1.000000-Full dose images-59704/') -# Load projection data -print("Loading projection data from {:s}".format(proj_folder)) -geometry, proj_data = mayo.load_projections(proj_folder, +# Load projection data restricting to a central slice +print("Loading projection data from {:s}".format(proj_dir)) +geometry, proj_data = mayo.load_projections(proj_dir, indices=slice(16000, 19000)) # Load reconstruction data -print("Loading reference data from {:s}".format(rec_folder)) -recon_space, volume = mayo.load_reconstruction(rec_folder) +print("Loading reference data from {:s}".format(rec_dir)) +recon_space, volume = mayo.load_reconstruction(rec_dir) # ray transform ray_trafo = odl.tomo.RayTransform(recon_space, geometry) @@ -48,11 +57,6 @@ fbp_result_HU = (fbp_result-0.0192)/0.0192*1000 -# Save reconstruction in Numpy format -fbp_filename = proj_folder+'/fbp_result.npy' -print("Saving reconstruction data in {:s}".format(fbp_filename)) -np.save(fbp_filename, fbp_result_HU) - # Compare the computed recon to reference reconstruction (coronal slice) ref = recon_space.element(volume) diff = recon_space.element(volume - fbp_result_HU.asarray()) @@ -64,4 +68,4 @@ coords = [0, None, None] fbp_result_HU.show('Recon (sagittal)', coords=coords) -ref.show('Reference (sagittal)', coords=coords) \ No newline at end of file +ref.show('Reference (sagittal)', coords=coords) diff --git a/odl/contrib/datasets/ct/mayo.py b/odl/contrib/datasets/ct/mayo.py index 4577964143c..1ff20aa9eed 100644 --- a/odl/contrib/datasets/ct/mayo.py +++ b/odl/contrib/datasets/ct/mayo.py @@ -11,9 +11,9 @@ In addition to the standard ODL requirements, this library also requires: - tqdm - - dicom - - A copy of the Mayo dataset, see - https://www.aapm.org/GrandChallenge/LowDoseCT/#registration + - pydicom + - Samples from the Mayo dataset, see + https://ctcicblog.mayo.edu/hubcap/patient-ct-projection-data-library/ """ from __future__ import division @@ -23,6 +23,7 @@ import pydicom import odl import tqdm +from functools import partial from pydicom.datadict import DicomDictionary, keyword_dict from odl.discr.discr_utils import linear_interpolator @@ -38,16 +39,16 @@ __all__ = ('load_projections', 'load_reconstruction') -def _read_projections(folder, indices): - """Read mayo projections from a folder.""" +def _read_projections(dir, indices): + """Read mayo projections from a directory.""" datasets = [] data_array = [] # Get the relevant file names - file_names = sorted([f for f in os.listdir(folder) if f.endswith(".dcm")]) + file_names = sorted([f for f in os.listdir(dir) if f.endswith(".dcm")]) if len(file_names) == 0: - raise ValueError('No DICOM files found in {}'.format(folder)) + raise ValueError('No DICOM files found in {}'.format(dir)) file_names = file_names[indices] @@ -55,12 +56,11 @@ def _read_projections(folder, indices): 'Loading projection data')): # read the file try: - dataset = pydicom.read_file(folder + os.path.sep + file_name) + dataset = pydicom.read_file(os.path.join(dir, file_name)) except: print("corrupted file: {}".format(file_name), file=sys.stderr) print("error:\n{}".format(sys.exc_info()[1]), file=sys.stderr) - print("skipping to next file..", file=sys.stderr) - continue + raise if not data_array: # Get some required data @@ -76,7 +76,6 @@ def _read_projections(folder, indices): assert rescale_intercept == dataset.RescaleIntercept assert rescale_slope == dataset.RescaleSlope - # Load the array as bytes proj_array = np.array(np.frombuffer(dataset.PixelData, 'H'), dtype='float32') @@ -84,7 +83,7 @@ def _read_projections(folder, indices): data_array.append(proj_array[:, ::-1]) datasets.append(dataset) - data_array = np.array(data_array) + data_array = np.stack(data_array) # Rescale array data_array *= rescale_slope data_array += rescale_intercept @@ -92,13 +91,13 @@ def _read_projections(folder, indices): return datasets, data_array -def load_projections(folder, indices=None, use_ffs=True): - """Load geometry and data stored in Mayo format from folder. +def load_projections(dir, indices=None, use_ffs=True): + """Load geometry and data stored in Mayo format from dir. Parameters ---------- - folder : str - Path to the folder where the Mayo DICOM files are stored. + dir : str + Path to the directory where the Mayo DICOM files are stored. indices : optional Indices of the projections to load. Accepts advanced indexing such as slice or list of indices. @@ -114,7 +113,7 @@ def load_projections(folder, indices=None, use_ffs=True): Projection data, given as the line integral of the linear attenuation coefficient (g/cm^3). Its unit is thus g/cm^2. """ - datasets, data_array = _read_projections(folder, indices) + datasets, data_array = _read_projections(dir, indices) # Get the angles angles = np.array([d.DetectorFocalCenterAngularPosition for d in datasets]) @@ -139,8 +138,8 @@ def load_projections(folder, indices=None, use_ffs=True): # For unknown reasons, mayo does not include the tag # "TableFeedPerRotation", which is what we want. # Instead we manually compute the pitch - table_dist = datasets[-1].DetectorFocalCenterAxialPosition - \ - datasets[0].DetectorFocalCenterAxialPosition + table_dist = (datasets[-1].DetectorFocalCenterAxialPosition - + datasets[0].DetectorFocalCenterAxialPosition) num_rot = (angles[-1] - angles[0]) / (2 * np.pi) pitch = table_dist / num_rot @@ -160,8 +159,8 @@ def load_projections(folder, indices=None, use_ffs=True): detector_partition = odl.uniform_partition(det_minp, det_maxp, det_shape) # Convert offset to odl definitions - offset_along_axis = datasets[0].DetectorFocalCenterAxialPosition - \ - angles[0] / (2 * np.pi) * pitch + offset_along_axis = (datasets[0].DetectorFocalCenterAxialPosition - + angles[0] / (2 * np.pi) * pitch) # Assemble geometry angle_partition = odl.nonuniform_partition(angles) @@ -169,8 +168,10 @@ def load_projections(folder, indices=None, use_ffs=True): # Flying focal spot src_shift_func = None if use_ffs: - src_shift_func = lambda angle: odl.tomo.flying_focal_spot( - angle, apart=angle_partition, shifts=shifts) + src_shift_func = partial( + odl.tomo.flying_focal_spot, apart=angle_partition, shifts=shifts) + else: + src_shift_func = None geometry = odl.tomo.ConeBeamGeometry(angle_partition, detector_partition, @@ -203,8 +204,8 @@ def interpolate_flat_grid(data_array, range_grid, radial_dist): """ # convert coordinates - theta, up, vp = range_grid.meshgrid #ray_trafo.range.grid.meshgrid - # d = src_radius + det_radius + theta, up, vp = range_grid.meshgrid + u = radial_dist * np.arctan(up / radial_dist) v = radial_dist / np.sqrt(radial_dist**2 + up**2) * vp @@ -218,13 +219,13 @@ def interpolate_flat_grid(data_array, range_grid, radial_dist): return proj_data -def load_reconstruction(folder, slice_start=0, slice_end=-1): - """Load a volume from folder, also returns the corresponding partition. +def load_reconstruction(dir, slice_start=0, slice_end=-1): + """Load a volume from dir, also returns the corresponding partition. Parameters ---------- - folder : str - Path to the folder where the DICOM files are stored. + dir : str + Path to the directory where the DICOM files are stored. slice_start : int Index of the first slice to use. Used for subsampling. slice_end : int @@ -249,10 +250,10 @@ def load_reconstruction(folder, slice_start=0, slice_end=-1): This function should handle all of these peculiarities and give a volume with the correct coordinate system attached. """ - file_names = sorted([f for f in os.listdir(folder) if f.endswith(".dcm")]) + file_names = sorted([f for f in os.listdir(dir) if f.endswith(".dcm")]) if len(file_names) == 0: - raise ValueError('No DICOM files found in {}'.format(folder)) + raise ValueError('No DICOM files found in {}'.format(dir)) volumes = [] datasets = [] @@ -261,7 +262,7 @@ def load_reconstruction(folder, slice_start=0, slice_end=-1): for file_name in tqdm.tqdm(file_names, 'loading volume data'): # read the file - dataset = pydicom.read_file(folder + os.path.sep + file_name) + dataset = pydicom.read_file(os.path.join(dir, file_name)) # Get parameters pixel_size = np.array(dataset.PixelSpacing)