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.

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

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 3
      1 import numpy as np
      2 from tempfile import mkdtemp
----> 3 from nabu.resources.dataset_analyzer import analyze_dataset
      4 from nabu.resources.nxflatfield import update_dataset_info_flats_darks
      5 from nabu.testutils import get_file

ModuleNotFoundError: No module named 'nabu'

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) ...

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[2], line 2
      1 print("Getting dataset (downloading if necessary) ...")
----> 2 data_path = get_file("bamboo_reduced.nx")
      3 print("... OK")
      4 output_dir = mkdtemp(prefix="nabu_reconstruction")

NameError: name 'get_file' is not defined

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[3], line 4
      1 # Parse the ".nx" file. This NX file is our entry point when it comes to data,
      2 # as it's only the format which is remotely stable
      3 # From this .nx file, build a data structure with readily available information
----> 4 dataset_info = analyze_dataset(data_path)

NameError: name 'analyze_dataset' is not defined

Estimate the center of rotation

from nabu.pipeline.estimators import estimate_cor

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[4], line 1
----> 1 from nabu.pipeline.estimators import estimate_cor

ModuleNotFoundError: No module named 'nabu'

cor = estimate_cor(
    "sliding-window", # estimator name
    dataset_info,
    do_flatfield=True,
)
print("Estimated center of rotation: %.2f" % cor)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[5], line 1
----> 1 cor = estimate_cor(
      2     "sliding-window", # estimator name
      3     dataset_info,
      4     do_flatfield=True,
      5 )
      6 print("Estimated center of rotation: %.2f" % cor)

NameError: name 'estimate_cor' is not defined

Define how the data should be processed

Instantiate the various components of the pipeline.

from nabu.io.reader import NXTomoReader
from nabu.preproc.flatfield import FlatField
from nabu.preproc.phase import PaganinPhaseRetrieval
from nabu.reconstruction.fbp import Backprojector

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[6], line 1
----> 1 from nabu.io.reader import NXTomoReader
      2 from nabu.preproc.flatfield import FlatField
      3 from nabu.preproc.phase import PaganinPhaseRetrieval

ModuleNotFoundError: No module named 'nabu'

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

reader = NXTomoReader(
    data_path,
    sub_region=sub_region,
)
radios_shape = reader.output_shape
n_angles, n_z, n_x = radios_shape

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[8], line 1
----> 1 reader = NXTomoReader(
      2     data_path,
      3     sub_region=sub_region,
      4 )
      5 radios_shape = reader.output_shape
      6 n_angles, n_z, n_x = radios_shape

NameError: name 'NXTomoReader' is not defined

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())),
)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[9], line 1
----> 1 flatfield = FlatField(
      2     radios_shape,
      3     dataset_info.get_reduced_flats(sub_region=sub_region),
      4     dataset_info.get_reduced_darks(sub_region=sub_region),
      5     radios_indices=sorted(list(dataset_info.projections.keys())),
      6 )

NameError: name 'FlatField' is not defined

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,
)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[10], line 1
----> 1 paganin = PaganinPhaseRetrieval(
      2     radios_shape[1:],
      3     distance=dataset_info.distance,
      4     energy=dataset_info.energy,
      5     delta_beta=250., # should be tuned
      6     pixel_size=dataset_info.pixel_size * 1e-6,
      7 )

NameError: name 'PaganinPhaseRetrieval' is not defined

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}
)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[11], line 1
----> 1 reconstruction = Backprojector(
      2     (n_angles, n_x),
      3     angles=dataset_info.rotation_angles,
      4     rot_center=cor,
      5     halftomo=dataset_info.is_halftomo,
      6     padding_mode="edges",
      7     scale_factor=1/(dataset_info.pixel_size * 1e-4), # mu/cm
      8     extra_options={"centered_axis": True, "clip_outer_circle": True}
      9 )

NameError: name 'Backprojector' is not defined

Read data

radios = reader.load_data()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[12], line 1
----> 1 radios = reader.load_data()

NameError: name 'reader' is not defined

Pre-processing

flatfield.normalize_radios(radios) # in-place
print()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[13], line 1
----> 1 flatfield.normalize_radios(radios) # in-place
      2 print()

NameError: name 'flatfield' is not defined

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[14], line 1
----> 1 radios_phase = np.zeros_like(radios)
      2 for i in range(radios.shape[0]):
      3     paganin.retrieve_phase(radios[i], output=radios_phase[i])

NameError: name 'radios' is not defined

Reconstruction

volume = np.zeros((n_z, n_x, n_x), "f")

for i in range(n_z):
    volume[i] = reconstruction.fbp(radios[:, i, :])

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[15], line 1
----> 1 volume = np.zeros((n_z, n_x, n_x), "f")
      3 for i in range(n_z):
      4     volume[i] = reconstruction.fbp(radios[:, i, :])

NameError: name 'n_z' is not defined

import matplotlib.pyplot as plt

plt.figure()
plt.imshow(volume[0], cmap="gray")
plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[17], line 2
      1 plt.figure()
----> 2 plt.imshow(volume[0], cmap="gray")
      3 plt.show()

NameError: name 'volume' is not defined

<Figure size 640x480 with 0 Axes>

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.