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.

In [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

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

In [2]:

print("Getting dataset (downloading if necessary) ...")
data_path = get_file("bamboo_reduced.nx")
print("... OK")
output_dir = mkdtemp(prefix="nabu_reconstruction")

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

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

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

In [4]:

from nabu.pipeline.estimators import estimate_cor

from nabu.pipeline.estimators import estimate_cor

In [5]:

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

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.

In [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

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

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

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

In [8]:

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

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

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

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

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

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

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

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

Could not get the current device. Was a CUDA context created ?

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

File ~/.venv/py313/lib/python3.13/site-packages/nabu/reconstruction/fbp_base.py:128, in BackprojectorBase.__init__(self, sino_shape, slice_shape, angles, rot_center, padding_mode, halftomo, filter_name, slice_roi, scale_factor, extra_options, backend_options)
    126 self._check_textures_availability()
    127 self._init_geometry(sino_shape, slice_shape, angles, rot_center, halftomo, slice_roi)
--> 128 self._init_filter(filter_name)
    129 self._allocate_memory()
    130 self._compute_angles()

File ~/.venv/py313/lib/python3.13/site-packages/nabu/reconstruction/fbp_base.py:274, in BackprojectorBase._init_filter(self, filter_name)
    272 #
    273 sinofilter_other_kwargs = self._get_filter_init_extra_options()
--> 274 self.sino_filter = self.SinoFilterClass(
    275     self.sino_shape,
    276     filter_name=self.filter_name,
    277     padding_mode=self.padding_mode,
    278     extra_options={"cutoff": self.extra_options.get("filter_cutoff", 1.0)},
    279     **sinofilter_other_kwargs,
    280 )
    281 if self.halftomo:
    282     # When doing half-tomography, each projections is seen "twice".
    283     # SinoFilter normalizes with pi/n_angles, but in half-tomography here n_angles is somehow halved.
    284     # TODO it should even be "n_turns", where n_turns can be computed from the angles
    285     self.sino_filter.set_filter(self.sino_filter.filter_f * (self.n_angles / np.pi * 2))

File ~/.venv/py313/lib/python3.13/site-packages/nabu/reconstruction/filtering_cuda.py:24, in CudaSinoFilter.__init__(self, sino_shape, filter_name, padding_mode, crop_filtered_data, extra_options, cuda_options)
     22 self._cuda_options = cuda_options or {}
     23 self.cuda = CudaProcessing(**self._cuda_options)
---> 24 super().__init__(
     25     sino_shape,
     26     filter_name=filter_name,
     27     padding_mode=padding_mode,
     28     crop_filtered_data=crop_filtered_data,
     29     extra_options=extra_options,
     30 )
     31 self._init_kernels()

File ~/.venv/py313/lib/python3.13/site-packages/nabu/reconstruction/filtering.py:76, in SinoFilter.__init__(self, sino_shape, filter_name, padding_mode, crop_filtered_data, extra_options)
     74 self._set_padding_mode(padding_mode)
     75 self._calculate_shapes(sino_shape, crop_filtered_data)
---> 76 self._init_fft()
     77 self._allocate_memory()
     78 self._compute_filter(filter_name)

File ~/.venv/py313/lib/python3.13/site-packages/nabu/reconstruction/filtering_cuda.py:35, in CudaSinoFilter._init_fft(self)
     33 def _init_fft(self):
     34     fft_cls = get_fft_class(self.extra_options["fft_backend"])
---> 35     self.fft = fft_cls(
     36         self.sino_padded_shape,
     37         dtype=np.float32,
     38         axes=(-1,),
     39     )

File ~/.venv/py313/lib/python3.13/site-packages/nabu/processing/fft_base.py:48, in _BaseFFT.__init__(self, shape, dtype, r2c, axes, normalize, **backend_options)
     46 self._configure_batched_transform()
     47 self._configure_normalization(normalize)
---> 48 self._compute_fft_plans()

File ~/.venv/py313/lib/python3.13/site-packages/nabu/processing/fft_base.py:135, in _BaseVKFFT._compute_fft_plans(self)
    134 def _compute_fft_plans(self):
--> 135     self._vkfft_plan = self.get_fft_obj(
    136         self.shape,
    137         self.dtype,
    138         ndim=self._vkfft_ndim,
    139         inplace=False,
    140         norm=self._vkfft_norm,
    141         r2c=self.r2c,
    142         dct=False,
    143         axes=self.axes,
    144         strides=None,
    145         **self._vkfft_other_init_kwargs,
    146     )

File ~/.venv/py313/lib/python3.13/site-packages/nabu/processing/fft_cuda.py:30, in get_vkfft_cuda(slf, *args, **kwargs)
     29 def get_vkfft_cuda(slf, *args, **kwargs):
---> 30     return _get_vkfft_cuda(*args, **kwargs)

File ~/.venv/py313/lib/python3.13/site-packages/nabu/processing/fft_cuda.py:26, in _get_vkfft_cuda(*args, **kwargs)
     24 @maybe_cached
     25 def _get_vkfft_cuda(*args, **kwargs):
---> 26     return CudaVkFFTApp(*args, **kwargs)

File ~/.venv/py313/lib/python3.13/site-packages/pyvkfft/cuda.py:185, in VkFFTApp.__init__(self, shape, dtype, ndim, inplace, stream, norm, r2c, dct, dst, axes, strides, tune_config, r2c_odd, verbose, **kwargs)
    183 self.config = self._make_config()
    184 if self.config is None:
--> 185     raise RuntimeError("Error creating VkFFTConfiguration. Was the CUDA context properly initialised ?")
    187 res = ctypes.c_int(0)
    188 # Size of tmp buffer allocated by VkFFT - if any

RuntimeError: Error creating VkFFTConfiguration. Was the CUDA context properly initialised ?

Read data¶

In [12]:

radios = reader.load_data()

radios = reader.load_data()

Pre-processing¶

In [13]:

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

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

In [14]:

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

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

Reconstruction¶

In [15]:

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

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

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 4
      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 'reconstruction' is not defined

In [16]:

import matplotlib.pyplot as plt

import matplotlib.pyplot as plt

In [17]:

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

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.

In [ ]: