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

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

# 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_reconstructionpe9t7o54/bamboo_reduced_darks.hdf5
Could not load flats from /tmp/nabu_reconstructionpe9t7o54/bamboo_reduced_flats.hdf5
Could not load darks from /tmp/nabu_testdata_gitlab-runner/bamboo_reduced_darks.hdf5
Could not load flats from /tmp/nabu_testdata_gitlab-runner/bamboo_reduced_flats.hdf5

Saved reduced darks to /tmp/nabu_reconstructionpe9t7o54/bamboo_reduced_darks.hdf5
Saved reduced flats to /tmp/nabu_reconstructionpe9t7o54/bamboo_reduced_flats.hdf5

Estimate the center of rotation

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'

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.

from nabu.io.reader import ChunkReader
from nabu.preproc.flatfield import FlatFieldDataUrls
from nabu.preproc.phase import PaganinPhaseRetrieval
from nabu.reconstruction.fbp import Backprojector

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-6-aeb31a287050> in <module>
      1 from nabu.io.reader import ChunkReader
----> 2 from nabu.preproc.flatfield import FlatFieldDataUrls
      3 from nabu.preproc.phase import PaganinPhaseRetrieval
      4 from nabu.reconstruction.fbp import Backprojector

~/.local/lib/python3.8/site-packages/nabu/preproc/__init__.py in <module>
      1 from .ccd import CCDFilter, Log
      2 from .ctf import CTFPhaseRetrieval
----> 3 from .distortion import DistortionCorrection
      4 from .double_flatfield import DoubleFlatField
      5 from .flatfield import FlatField, FlatFieldDataUrls

~/.local/lib/python3.8/site-packages/nabu/preproc/distortion.py in <module>
      2 from scipy.interpolate import RegularGridInterpolator
      3 from ..utils import check_supported
----> 4 from ..estimation.distortion import estimate_flat_distortion
      5
      6

~/.local/lib/python3.8/site-packages/nabu/estimation/__init__.py in <module>
----> 1 from .alignment import AlignmentBase
      2 from .cor import (
      3     CenterOfRotation,
      4     CenterOfRotationSlidingWindow,
      5     CenterOfRotationGrowingWindow,

~/.local/lib/python3.8/site-packages/nabu/estimation/alignment.py in <module>
      4 from numpy.polynomial.polynomial import Polynomial
      5 from silx.math.medianfilter import medfilt2d
----> 6 import scipy.fft  # pylint: disable=E0611
      7 from ..utils import previouspow2
      8 from ..misc import fourier_filters

ModuleNotFoundError: No module named 'scipy.fft'

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

reader = ChunkReader(
    dataset_info.projections,
    sub_region=sub_region,
    convert_float=True,
)
radios = reader.data
n_angles, n_z, n_x = radios.shape

flatfield = FlatFieldDataUrls(
    radios.shape,
    dataset_info.flats,
    dataset_info.darks,
    radios_indices=sorted(list(dataset_info.projections.keys())),
    sub_region=sub_region
)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-9-3f6404a148bb> in <module>
----> 1 flatfield = FlatFieldDataUrls(
      2     radios.shape,
      3     dataset_info.flats,
      4     dataset_info.darks,
      5     radios_indices=sorted(list(dataset_info.projections.keys())),

NameError: name 'FlatFieldDataUrls' is not defined

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-10-b0121f7c548f> in <module>
----> 1 paganin = PaganinPhaseRetrieval(
      2     radios[0].shape,
      3     distance=dataset_info.distance,
      4     energy=dataset_info.energy,
      5     delta_beta=250., # should be tuned

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

reader.load_data(overwrite=True)

Pre-processing

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

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

NameError: name 'reconstruction' is not defined

import matplotlib.pyplot as plt

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.