From 5765d0a1f0c8cd122281490f92d962782bd7a713 Mon Sep 17 00:00:00 2001 From: pierocor Date: Tue, 2 Feb 2021 16:11:08 +0100 Subject: [PATCH] STY: pull request mods --- odl/contrib/datasets/README.md | 2 +- .../datasets/ct/examples/mayo_reconstruct.py | 36 +++++----- odl/contrib/datasets/ct/mayo.py | 65 ++++++++++--------- 3 files changed, 53 insertions(+), 50 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..a2c93fb3746 100644 --- a/odl/contrib/datasets/ct/examples/mayo_reconstruct.py +++ b/odl/contrib/datasets/ct/examples/mayo_reconstruct.py @@ -1,29 +1,36 @@ """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. """ 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, +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 +55,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 +66,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..5bff56fedcd 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)