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 CORFinder
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-4-e6b532b7cb10> in <module>
----> 1 from nabu.pipeline.estimators import CORFinder

~/.local/lib/python3.8/site-packages/nabu/pipeline/estimators.py in <module>
      5 import inspect
      6 import numpy as np
----> 7 import scipy.fft  #  pylint: disable=E0611
      8 from silx.io import get_data
      9 from ..preproc.flatfield import FlatFieldDataUrls

ModuleNotFoundError: No module named 'scipy.fft'
[5]:
cor_finder = CORFinder(
    "sliding-window",
    dataset_info,
    do_flatfield=True,
)
cor = cor_finder.find_cor()
print("Estimated center of rotation: %.2f" % cor)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-5-517af09c15c6> in <module>
----> 1 cor_finder = CORFinder(
      2     "sliding-window",
      3     dataset_info,
      4     do_flatfield=True,
      5 )

NameError: name 'CORFinder' is not defined

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
---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-6-b70f0fd6c8e2> in <module>
----> 1 from nabu.io.reader import NXTomoReader
      2 from nabu.preproc.flatfield import FlatField
      3 from nabu.preproc.phase import PaganinPhaseRetrieval
      4 from nabu.reconstruction.fbp import Backprojector

ImportError: cannot import name 'NXTomoReader' from 'nabu.io.reader' (/var/lib/gitlab-runner/.local/lib/python3.8/site-packages/nabu/io/reader.py)
[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
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-8-3f7b6efb6350> in <module>
----> 1 reader = NXTomoReader(
      2     data_path,
      3     sub_region=sub_region,
      4 )
      5 radios_shape = reader.output_shape

NameError: name 'NXTomoReader' is not defined
[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())),
)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-9-43a589e0d43e> in <module>
----> 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())),

NameError: name 'FlatField' is not defined
[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,
)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-10-4ae2eb148aef> in <module>
----> 1 paganin = PaganinPhaseRetrieval(
      2     radios_shape[1:],
      3     distance=dataset_info.distance,
      4     energy=dataset_info.energy,
      5     delta_beta=250., # should be tuned

NameError: name 'PaganinPhaseRetrieval' is not defined
[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}
)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-11-c08982cb1f66> in <module>
----> 1 reconstruction = Backprojector(
      2     (n_angles, n_x),
      3     angles=dataset_info.rotation_angles,
      4     rot_center=cor,
      5     halftomo=dataset_info.is_halftomo,

NameError: name 'Backprojector' is not defined

Read data

[12]:
radios = reader.load_data()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-12-ad4922173686> in <module>
----> 1 radios = reader.load_data()

NameError: name 'reader' is not defined

Pre-processing

[13]:
flatfield.normalize_radios(radios) # in-place
print()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-13-b9e4a1f57ebe> in <module>
----> 1 flatfield.normalize_radios(radios) # in-place
      2 print()

NameError: name 'flatfield' is not defined
[14]:
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)
<ipython-input-14-8b41c1a1caa6> in <module>
----> 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

[15]:
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)
<ipython-input-15-00222c7e518f> in <module>
----> 1 volume = np.zeros((n_z, n_x, n_x), "f")
      2
      3 for i in range(n_z):
      4     volume[i] = reconstruction.fbp(radios[:, i, :])

NameError: name 'n_z' is not defined
[16]:
import matplotlib.pyplot as plt
[17]:
plt.figure()
plt.imshow(volume[0], cmap="gray")
plt.show()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-17-394cfad5f2ae> in <module>
      1 plt.figure()
----> 2 plt.imshow(volume[0], cmap="gray")
      3 plt.show()

NameError: name 'volume' is not defined
<Figure size 432x288 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.

[ ]: