Reconstructing a complete scan using nabu and config file, from python¶

In previous notebooks nabu_basic_reconstruction.ipynb and nabu_from_python_with_gpu.ipynb, we saw how to manipulate the individual processing classes for step-by-step processing (eg. FlatField, PaganinPhaseRetrieval ; and their cuda counterparts).

In this notebook, we see how to use a facility class FullFieldReconstructor which assembles all these steps, does the "boring work" for you (memory allocation, data transfers, configuration ingestion, ...). It reconstructs the input volume by parts, handling the necessary overlap between sub-volumes if needed.

1. Get the scan (.nx file) and processing configuration (.conf file)¶

In [1]:

import os
from nabu.testutils import utilstest, get_file
from nabu.pipeline.fullfield.processconfig import ProcessConfig

import os from nabu.testutils import utilstest, get_file from nabu.pipeline.fullfield.processconfig import ProcessConfig

In [2]:

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

# Get the configuration file of this dataset
conf_fname = get_file("bamboo_reduced.conf")

# Change directory to the path where the data is located (only useful for this tutorial)
os.chdir(utilstest.data_home)

# Parse this configuration file
conf = ProcessConfig(conf_fname)

print("Getting dataset (downloading if necessary) ...") data_path = get_file("bamboo_reduced.nx") print("... OK") # Get the configuration file of this dataset conf_fname = get_file("bamboo_reduced.conf") # Change directory to the path where the data is located (only useful for this tutorial) os.chdir(utilstest.data_home) # Parse this configuration file conf = ProcessConfig(conf_fname)

Getting dataset (downloading if necessary) ...

... OK
Option 'double_flatfield_enabled' has been renamed 'double_flatfield' in [preproc]
This is deprecated since version 2025.1.0 and will result in an error in futures versions
Browsing dataset
Updating dataset information with user configuration
Loaded darks from /tmp/nabu_testdata_pierre/bamboo_reduced_darks.hdf5
Loaded flats from /tmp/nabu_testdata_pierre/bamboo_reduced_flats.hdf5
Doing dataset estimations
Could not get an initial estimate for center of rotation in data file
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.986
Doing coupled validation
Cannot do SRCurrent normalization: missing flats and/or projections SRCurrent

In [ ]:

2. Running volume reconstruction¶

You can specify which reconstruction backend to use ("numpy" or "cuda")

In [3]:

from nabu.pipeline.fullfield.reconstruction import FullFieldReconstructor 

from nabu.pipeline.fullfield.reconstruction import FullFieldReconstructor

In [4]:

reconstructor = FullFieldReconstructor(conf, backend="cuda")

reconstructor = FullFieldReconstructor(conf, backend="cuda")

Estimated margin: 8 pixels

In [5]:

reconstructor.reconstruct()

reconstructor.reconstruct()

Processing sub-volume ((0, 1000), (0, 242), (0, 640))

[0] Creating a new pipeline object

[0] Set sub-region to ((0, 1000), (149, 391), (0, 640))

[0] Allocating recs: (234, np.int64(640), np.int64(640))

[0] Allocating tmp_sino: (1000, 640)

[0] Allocating radios: (1000, 242, 640)

[0] Set sub-region to ((0, 1000), (0, 242), (0, 640))

[0] File already exists: /tmp/nabu_testdata_pierre/bamboo_reduced_rec/bamboo_reduced_rec_00000.hdf5. It will be overwritten as requested in configuration

[0] Reading data

[0] Region = ((0, 1000), (0, 242), (0, 640))

[0] Read subvolume (1000, 242, 640) in 0.13 s

[0] Transfering radios to GPU

[0] Applying flat-field

[0] Performing phase retrieval

[0] Taking logarithm

[0] Cropping radios from (1000, 242, 640) to (1000, 234, 640)

[0] Removing rings on sinograms

[0] Reconstruction

[0] Saving data

[0] Wrote /tmp/nabu_testdata_pierre/bamboo_reduced_rec/bamboo_reduced_rec_00000.hdf5

[0] Processing sub-volume ((0, 1000), (226, 476), (0, 640))

[0] Destroying pipeline instance and releasing memory

[0] Creating a new pipeline object

[0] Set sub-region to ((0, 1000), (145, 395), (0, 640))

[0] Allocating recs: (234, np.int64(640), np.int64(640))

[0] Allocating tmp_sino: (1000, 640)

[0] Allocating radios: (1000, 250, 640)

[0] Set sub-region to ((0, 1000), (226, 476), (0, 640))

[0] Reading data

[0] Region = ((0, 1000), (226, 476), (0, 640))

[0] Read subvolume (1000, 250, 640) in 0.14 s

[0] Transfering radios to GPU

[0] Applying flat-field

[0] Performing phase retrieval

[0] Taking logarithm

[0] Cropping radios from (1000, 250, 640) to (1000, 234, 640)

[0] Removing rings on sinograms

[0] Reconstruction

[0] Saving data

[0] Wrote /tmp/nabu_testdata_pierre/bamboo_reduced_rec/bamboo_reduced_rec_00234.hdf5

[0] Processing sub-volume ((0, 1000), (460, 540), (0, 640))

[0] Destroying pipeline instance and releasing memory

[0] Creating a new pipeline object

[0] Set sub-region to ((0, 1000), (230, 310), (0, 640))

[0] Allocating recs: (72, np.int64(640), np.int64(640))

[0] Allocating tmp_sino: (1000, 640)

[0] Allocating radios: (1000, 80, 640)

[0] Set sub-region to ((0, 1000), (460, 540), (0, 640))

[0] Reading data

[0] Region = ((0, 1000), (460, 540), (0, 640))

[0] Read subvolume (1000, 80, 640) in 0.05 s

[0] Transfering radios to GPU

[0] Applying flat-field

[0] Performing phase retrieval

[0] Taking logarithm

[0] Cropping radios from (1000, 80, 640) to (1000, 72, 640)

[0] Removing rings on sinograms

[0] Reconstruction

[0] Saving data

[0] Wrote /tmp/nabu_testdata_pierre/bamboo_reduced_rec/bamboo_reduced_rec_00468.hdf5

In [6]:

output_file = reconstructor.merge_hdf5_reconstructions()

output_file = reconstructor.merge_hdf5_reconstructions()

File /tmp/nabu_testdata_pierre/bamboo_reduced_rec.hdf5 already exists. Overwriting as requested in configuration file
Merging reconstructions to /tmp/nabu_testdata_pierre/bamboo_reduced_rec.hdf5

3. Visualization¶

In [7]:

from tomoscan.io import HDF5File

from tomoscan.io import HDF5File

In [8]:

with HDF5File(output_file, "r") as f:
    data_ptr = f["entry0000/reconstruction/results/data"]
    middle_slice_horiz = data_ptr[reconstructor.delta_z//2]
    middle_slice_vertic = data_ptr[:, :, reconstructor.n_x//2]
    

with HDF5File(output_file, "r") as f: data_ptr = f["entry0000/reconstruction/results/data"] middle_slice_horiz = data_ptr[reconstructor.delta_z//2] middle_slice_vertic = data_ptr[:, :, reconstructor.n_x//2]

In [9]:

import matplotlib.pyplot as plt
plt.figure()
plt.subplot(121)
plt.imshow(middle_slice_horiz, cmap="gray")
plt.xlabel("Horizontal slice, middle")
plt.subplot(122)
plt.imshow(middle_slice_vertic, cmap="gray")
plt.xlabel("Vertical slice, middle")
plt.show()

import matplotlib.pyplot as plt plt.figure() plt.subplot(121) plt.imshow(middle_slice_horiz, cmap="gray") plt.xlabel("Horizontal slice, middle") plt.subplot(122) plt.imshow(middle_slice_vertic, cmap="gray") plt.xlabel("Vertical slice, middle") plt.show()

In [ ]: