Basic reconstruction with nabu¶
In this notebook, we see how to use the different building blocks of nabu to parse a dataset and perform a basic reconstruction.
This notebook uses a dataset of a bamboo stick (acquired at ESRF ID19, courtesy Ludovic Broche). The projections were binned by a factor of 4 in each dimension, and the angular range was also subsampled by 4.
Browse the dataset¶
We use nabu utility analyze_dataset
which builds a data structure with all information on the dataset.
[1]:
import numpy as np
from tempfile import mkdtemp
from nabu.resources.dataset_analyzer import analyze_dataset
from nabu.resources.nxflatfield import update_dataset_info_flats_darks
from nabu.testutils import get_file
[2]:
print("Getting dataset (downloading if necessary) ...")
data_path = get_file("bamboo_reduced.nx")
print("... OK")
output_dir = mkdtemp(prefix="nabu_reconstruction")
Getting dataset (downloading if necessary) ...
... OK
[3]:
# Parse the ".nx" file. This NX file is our entry point when it comes to data,
# as it's only the format which is remotely stable
# From this .nx file, build a data structure with readily available information
dataset_info = analyze_dataset(data_path)
Estimate the center of rotation¶
[4]:
from nabu.pipeline.estimators import estimate_cor
[5]:
cor = estimate_cor(
"sliding-window", # estimator name
dataset_info,
do_flatfield=True,
)
print("Estimated center of rotation: %.2f" % cor)
Estimating center of rotation
CenterOfRotationSlidingWindow.find_shift({'side': 'center', 'window_width': None, 'roi_yxhw': None, 'median_filt_shape': None, 'peak_fit_radius': 1, 'high_pass': None, 'low_pass': None, 'return_validity': False, 'return_relative_to_middle': False})
Estimated center of rotation: 338.99
Define how the data should be processed¶
Instantiate the various components of the pipeline.
[6]:
from nabu.io.reader import NXTomoReader
from nabu.preproc.flatfield import FlatField
from nabu.preproc.phase import PaganinPhaseRetrieval
from nabu.reconstruction.fbp import Backprojector
/home/pierre/.venv/py311/lib/python3.11/site-packages/pytools/persistent_dict.py:63: RecommendedHashNotFoundWarning: Unable to import recommended hash 'siphash24.siphash13', falling back to 'hashlib.sha256'. Run 'python3 -m pip install siphash24' to install the recommended hash.
warn("Unable to import recommended hash 'siphash24.siphash13', "
[7]:
# Define the sub-region to read (and reconstruct)
# Read all projections, from row 270 to row 288
sub_region = (
slice(None),
slice(270, 289),
slice(None)
)
[8]:
reader = NXTomoReader(
data_path,
sub_region=sub_region,
)
radios_shape = reader.output_shape
n_angles, n_z, n_x = radios_shape
[9]:
flatfield = FlatField(
radios_shape,
dataset_info.get_reduced_flats(sub_region=sub_region),
dataset_info.get_reduced_darks(sub_region=sub_region),
radios_indices=sorted(list(dataset_info.projections.keys())),
)
[10]:
paganin = PaganinPhaseRetrieval(
radios_shape[1:],
distance=dataset_info.distance,
energy=dataset_info.energy,
delta_beta=250., # should be tuned
pixel_size=dataset_info.pixel_size * 1e-6,
)
[11]:
reconstruction = Backprojector(
(n_angles, n_x),
angles=dataset_info.rotation_angles,
rot_center=cor,
halftomo=dataset_info.is_halftomo,
padding_mode="edges",
scale_factor=1/(dataset_info.pixel_size * 1e-4), # mu/cm
extra_options={"centered_axis": True, "clip_outer_circle": True}
)
/home/pierre/.venv/py311/lib/python3.11/site-packages/pytools/persistent_dict.py:63: RecommendedHashNotFoundWarning: Unable to import recommended hash 'siphash24.siphash13', falling back to 'hashlib.sha256'. Run 'python3 -m pip install siphash24' to install the recommended hash.
warn("Unable to import recommended hash 'siphash24.siphash13', "
/home/pierre/.venv/py311/lib/python3.11/site-packages/pytools/persistent_dict.py:63: RecommendedHashNotFoundWarning: Unable to import recommended hash 'siphash24.siphash13', falling back to 'hashlib.sha256'. Run 'python3 -m pip install siphash24' to install the recommended hash.
warn("Unable to import recommended hash 'siphash24.siphash13', "
/home/pierre/.venv/py311/lib/python3.11/site-packages/skcuda/cublas.py:284: UserWarning: creating CUBLAS context to get version number
warnings.warn('creating CUBLAS context to get version number')
Read data¶
[12]:
radios = reader.load_data()
Pre-processing¶
[13]:
flatfield.normalize_radios(radios) # in-place
print()
[14]:
radios_phase = np.zeros_like(radios)
for i in range(radios.shape[0]):
paganin.retrieve_phase(radios[i], output=radios_phase[i])
Reconstruction¶
[15]:
volume = np.zeros((n_z, n_x, n_x), "f")
for i in range(n_z):
volume[i] = reconstruction.fbp(radios[:, i, :])
[16]:
import matplotlib.pyplot as plt
[17]:
plt.figure()
plt.imshow(volume[0], cmap="gray")
plt.show()
Going further: GPU processing¶
All the above components have a Cuda backend: SinoBuilder
becomes CudaSinoBuilder
, PaganinPhaseRetrieval
becomes CudaPaganinPhaseRetrieval
, and so on. Importantly, the cuda counterpart of these classes have the same API.
[ ]: