STY: pull request mods

Author: pierocor

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`

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)
+ref.show('Reference (sagittal)', coords=coords)

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,29 +39,28 @@
 __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]
 
     for i, file_name in enumerate(tqdm.tqdm(file_names,
                                             '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,29 +76,28 @@ 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')
         proj_array = proj_array.reshape([cols, rows])
         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
 
     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,17 +159,19 @@ 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)
 
     # 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)