Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Maintenance of the Mayo clinic human CT dataset loader #1589

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Maintenance of the Mayo clinic human CT dataset loader #1589

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
f7e34a9
BUG: Load projection data and references from mayo dataset
pierocor Nov 26, 2020
8a5a197
FIX: Add flying focal spot
pierocor Dec 18, 2020
a851738
STY: Add documentation
pierocor Dec 18, 2020
96a7d5e
STY: pull request mods
pierocor Feb 2, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
66 changes: 43 additions & 23 deletions odl/contrib/datasets/ct/examples/mayo_reconstruct.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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 "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 ~48 000 projections in the full dataset.
there are ~32 000 projections per patient 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/'
rec_folder = odl.__path__[0] + '/../../data/LDCT-and-Projection-data/' \
'L004/08-21-2018-84608/1.000000-Full dose images-59704/'
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

  • I suggest to use odl.contrib.datasets.util.get_data_dir() to determine a storage location.
  • What are the reasons to change the dataset directory? It would be good to have a small comment on top that says which data it is.
  • Please use dir instead of folder. The latter is Windows lingo, and Python uses e.g. is_dir() to check if a Path object represents a directory.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The Mayo dataset cannot be downloaded directly from a url but it requires its own software. This prevents to use odl.contrib.datasets.util.get_data(filename, subset, url) for local caching. Would this work better?

Suggested change
# 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 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))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same question as above: why change the slice?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It has been adjusted to the example considered (patient L004). In this case, the projection data counts ~32k images and the above slice choice reconstructs a central part of the volume.

# Load reconstruction data
print("Loading reference data from {:s}".format(rec_folder))
recon_space, volume = mayo.load_reconstruction(rec_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)
Expand All @@ -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)
pierocor marked this conversation as resolved.
Show resolved Hide resolved

# 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 (axial)')
ref.show('Reference (axial)')
diff.show('Diff (axial)')

# 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)
Loading