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
WARNING:silx.DEPRECATION:Module silx.third_party is deprecated since silx version 2.0.0. Use 'fabio' instead.
  File "/home/pierre/.venv/py311/lib/python3.11/site-packages/silx/third_party/__init__.py", line 36, in <module>
    deprecated_warning(
WARNING:silx.DEPRECATION:Module silx.third_party.EdfFile is deprecated since silx version 2.0.0. Use 'fabio.open and fabio.edfimage.EdfImage' instead.
  File "/home/pierre/.venv/py311/lib/python3.11/site-packages/silx/third_party/EdfFile.py", line 118, in <module>
    deprecated_warning(
[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)

# A tomography scan usually contains a series of darks/flats.
# These darks/flats are then averaged (the legacy pipeline called them 'refHST' and 'darkend')
# The flat-field normalization should be done only with reduced flats/darks.
update_dataset_info_flats_darks(dataset_info, True, darks_flats_dir=output_dir)
Could not load darks from /tmp/nabu_reconstructionstm_lna6/bamboo_reduced_darks.hdf5
Could not load flats from /tmp/nabu_reconstructionstm_lna6/bamboo_reduced_flats.hdf5
Loaded darks from /tmp/nabu_testdata_pierre/bamboo_reduced_darks.hdf5
Loaded flats from /tmp/nabu_testdata_pierre/bamboo_reduced_flats.hdf5

Estimate the center of rotation

[4]:
from nabu.pipeline.estimators import CORFinder
[5]:
cor_finder = CORFinder(
    "sliding-window",
    dataset_info,
    do_flatfield=True,
)
cor = cor_finder.find_cor()
print("Estimated center of rotation: %.2f" % cor)
Estimating center of rotation
CenterOfRotationSlidingWindow.find_shift({'window_width': None, 'roi_yxhw': None, 'median_filt_shape': None, 'padding_mode': None, 'peak_fit_radius': 1, 'high_pass': None, 'low_pass': None, 'return_validity': False})
Estimated center of rotation: 339.49

Define how the data should be processed

Instantiate the various components of the pipeline.

[6]:
from nabu.io.reader import ChunkReader
from nabu.preproc.flatfield import FlatFieldDataUrls
from nabu.preproc.phase import PaganinPhaseRetrieval
from nabu.reconstruction.sinogram import SinoBuilder
from nabu.reconstruction.fbp import Backprojector
[7]:
# Define the sub-region to read (and reconstruct)
# In each radio: will read (start_x, end_x, start_y, end_y)
sub_region = (None, None, 270, 289)
[8]:
reader = ChunkReader(
    dataset_info.projections,
    sub_region=sub_region,
    convert_float=True,
)
radios = reader.data
[9]:
flatfield = FlatFieldDataUrls(
    radios.shape,
    dataset_info.flats,
    dataset_info.darks,
    radios_indices=sorted(list(dataset_info.projections.keys())),
    sub_region=sub_region
)
[10]:
paganin = PaganinPhaseRetrieval(
    radios[0].shape,
    distance=dataset_info.distance,
    energy=dataset_info.energy,
    delta_beta=250., # should be tuned
    pixel_size=dataset_info.pixel_size * 1e-6,
)
[11]:
sino_builder = SinoBuilder(
    radios_shape=radios.shape,
    rot_center=cor,
    halftomo=dataset_info.is_halftomo
)
sino_shape = sino_builder.output_shape[1:]
[12]:
reconstruction = Backprojector(
    sino_shape,
    angles=dataset_info.rotation_angles,
    rot_center=cor,
    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/skcuda/cublas.py:284: UserWarning: creating CUBLAS context to get version number
  warnings.warn('creating CUBLAS context to get version number')

Read data

[13]:
reader.load_data(overwrite=True)

Pre-processing

[14]:
flatfield.normalize_radios(radios) # in-place
print()

[15]:
radios_phase = np.zeros_like(radios)
for i in range(radios.shape[0]):
    paganin.retrieve_phase(radios[i], output=radios_phase[i])

Build sinogram

[16]:
# For 180 degrees scan, it's simply a transposition.
# For 360 degrees scan + half-acquisition, this duplicates data - not elegant but fine in this case
sinos = sino_builder.get_sinos(radios_phase)

Reconstruction

[17]:
volume = np.zeros((sinos.shape[0], sinos.shape[-1], sinos.shape[-1]), "f")

for i in range(sinos.shape[0]):
    volume[i] = reconstruction.fbp(sinos[0])
[18]:
import matplotlib.pyplot as plt
[19]:
plt.figure()
plt.imshow(volume[0], cmap="gray")
plt.show()
../_images/Tutorials_nabu_basic_reconstruction_26_0.png

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.

[ ]:

[ ]: