Using nabu from python to reconstruct a dataset with GPU

This notebook shows how to use the Nabu software for performing a basic reconstruction of a tomography dataset.
The computations are done on a local machine with a GPU and Cuda available.

This tutorial goes a bit further than nabu_basic_reconstruction.ipynb: - GPU implementation of each component is used - We see how to start from a configuration file and devise a simple processing chain accordingly

The same dataset is used (binned scan of a bamboo stick, thanks Ludovic Broche, ESRF ID19).

1 - Load the dataset informations

We must provide nabu with the the configuration file (nabu.conf), describing the path to the dataset and the processing steps. This is the equivalent of the .par file in PyHST2. In this file, no information is given on the detector size, energy, distance, etc: these informations are extracted from the dataset metadata.

[1]:

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

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-1-cd785c1233ec> in <module>
      1 import os
      2 from nabu.testutils import utilstest, get_file
----> 3 from nabu.pipeline.fullfield.processconfig import ProcessConfig

~/.local/lib/python3.8/site-packages/nabu/pipeline/fullfield/processconfig.py in <module>
     11 from ...resources.utils import get_quantities_and_units
     12 from ...reconstruction.sinogram import get_extended_sinogram_width
---> 13 from ..estimators import estimate_cor
     14 from ..processconfig import ProcessConfigBase
     15 from .nabu_config import nabu_config, renamed_keys

~/.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'

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

Getting dataset (downloading if necessary) ...
... OK

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-2-90f06e822805> in <module>
     10
     11 # Parse this configuration file
---> 12 conf = ProcessConfig(conf_fname)

NameError: name 'ProcessConfig' is not defined

Note that ProcessConfig will do quite a few things under the hood: - Parse the configuration file and check parameters correctness - Browse the dataset - Get or compute the reduced flats/darks - Estimate the center of rotation

The resulting object contains all necessary information to process the dataset.

[3]:

# We can easily get information on the processing steps.
nabu_config = conf.nabu_config
from pprint import pprint
pprint(nabu_config)
# The same can be done with the dataset structure
dataset_info = conf.dataset_info
# print([getattr(dataset_info, attr) for attr in ["energy", "distance", "n_angles", "radio_dims"]])

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-3-26803947f63b> in <module>
      1 # We can easily get information on the processing steps.
----> 2 nabu_config = conf.nabu_config
      3 from pprint import pprint
      4 pprint(nabu_config)
      5 # The same can be done with the dataset structure

NameError: name 'conf' is not defined

2 - Chunk processing

Nabu processes data by chunks of radios (see the documentation for more explanations).
In a first step, we define how to read chunks of radios.
[4]:

from nabu.io.reader import NXTomoReader

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-4-ae2d3bf3c298> in <module>
----> 1 from nabu.io.reader import NXTomoReader

ImportError: cannot import name 'NXTomoReader' from 'nabu.io.reader' (/var/lib/gitlab-runner/.local/lib/python3.8/site-packages/nabu/io/reader.py)
What is the largest chunk size we can process ?
The answer is given by inspecting the current GPU memory, and the processing steps.
[5]:

from nabu.cuda.utils import get_gpu_memory
from nabu.pipeline.fullfield.computations import estimate_max_chunk_size

[6]:

chunk_size = estimate_max_chunk_size(
    get_gpu_memory(0),
    conf
)
print("Chunk_size = %d" % chunk_size)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-6-85ce6a589ee8> in <module>
      1 chunk_size = estimate_max_chunk_size(
----> 2     get_gpu_memory(0),
      3     conf
      4 )
      5 print("Chunk_size = %d" % chunk_size)

~/.local/lib/python3.8/site-packages/nabu/cuda/utils.py in get_gpu_memory(device_id)
     69     Return the total memory (in GigaBytes) of a device.
     70     """
---> 71     cuda.init()
     72     return cuda.Device(device_id).total_memory() / 1e9
     73

NameError: name 'cuda' is not defined

[7]:

# Load the first 'chunk_size' lines of all the radios
# i.e do projections_data[:, 0:chunk_size, :]
sub_region = (
    slice(None),
    slice(0, chunk_size),
    slice(None)
)
projections_reader = NXTomoReader(
    data_path,
    sub_region=sub_region,
)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-7-041bb8d95bc5> in <module>
      3 sub_region = (
      4     slice(None),
----> 5     slice(0, chunk_size),
      6     slice(None)
      7 )

NameError: name 'chunk_size' is not defined

[8]:

# Load the current chunk
projections = projections_reader.load_data() # takes some time

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-8-911b1cf572aa> in <module>
      1 # Load the current chunk
----> 2 projections = projections_reader.load_data() # takes some time

NameError: name 'projections_reader' is not defined

[9]:

print(projections.shape)
print(projections.dtype)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-9-26c3dbef4737> in <module>
----> 1 print(projections.shape)
      2 print(projections.dtype)

NameError: name 'projections' is not defined

3 - Initialize the GPU

Most of the processing can be done on GPU (or many-core CPU if using OpenCL).
With pycuda.gpuarray (or its OpenCL counterpart pyopencl.array), we manipulate array objects with memory residing on device. This allows to avoid extraneous host <-> device copies.
[10]:

import pycuda.gpuarray as garray
from nabu.cuda.utils import get_cuda_context
import numpy as np

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-10-faa9caeb2ba0> in <module>
----> 1 import pycuda.gpuarray as garray
      2 from nabu.cuda.utils import get_cuda_context
      3 import numpy as np

ModuleNotFoundError: No module named 'pycuda'

[11]:

# Create a Cuda context on device ID 0
# By default, all following GPU processings will be bound on this context
ctx = get_cuda_context(device_id=0)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-11-e674a806a181> in <module>
      1 # Create a Cuda context on device ID 0
      2 # By default, all following GPU processings will be bound on this context
----> 3 ctx = get_cuda_context(device_id=0)

NameError: name 'get_cuda_context' is not defined

[12]:

n_angles, n_z, n_x = projections.shape
# transfer the chunk on GPU
d_radios = garray.to_gpu(projections)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-12-a12b819cbf3c> in <module>
----> 1 n_angles, n_z, n_x = projections.shape
      2 # transfer the chunk on GPU
      3 d_radios = garray.to_gpu(projections)

NameError: name 'projections' is not defined

4 - Pre-processing

Pre-processing utilities are available in the nabu.preproc module.
Utilities available with the cuda backend are implemented in a module with a _cuda suffix.

4.1 - Flat-field

[13]:

from nabu.preproc.flatfield_cuda import CudaFlatField

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-13-b96aa7a236ae> in <module>
----> 1 from nabu.preproc.flatfield_cuda import CudaFlatField

~/.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'

[14]:

radios_indices = sorted(conf.dataset_info.projections.keys())
# Configure the `FlatField` processor
cuda_flatfield = CudaFlatField(
    d_radios.shape,
    dataset_info.get_reduced_flats(sub_region=sub_region),
    dataset_info.get_reduced_darks(sub_region=sub_region),
    radios_indices=radios_indices,
)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-14-fd8007a9887a> in <module>
----> 1 radios_indices = sorted(conf.dataset_info.projections.keys())
      2 # Configure the `FlatField` processor
      3 cuda_flatfield = CudaFlatField(
      4     d_radios.shape,
      5     dataset_info.get_reduced_flats(sub_region=sub_region),

NameError: name 'conf' is not defined

[15]:

# Perform the normalization on GPU
if nabu_config["preproc"]["flatfield"]:
    print("Doing flat-field")
    cuda_flatfield.normalize_radios(d_radios)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-15-6b9a691df582> in <module>
      1 # Perform the normalization on GPU
----> 2 if nabu_config["preproc"]["flatfield"]:
      3     print("Doing flat-field")
      4     cuda_flatfield.normalize_radios(d_radios)

NameError: name 'nabu_config' is not defined

4.2 - Phase retrieval

[16]:

from nabu.preproc.phase_cuda import CudaPaganinPhaseRetrieval

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-16-e2d0f529fecf> in <module>
----> 1 from nabu.preproc.phase_cuda import CudaPaganinPhaseRetrieval

~/.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'

[17]:

energy = dataset_info.energy
# Phase retrieval is done on each radio individually, with the sub-region specified above
if (nabu_config["phase"]["method"] or "").lower() == "paganin":
    print("Doing phase retrieval")
    cudapaganin = CudaPaganinPhaseRetrieval(
        (n_z, n_x),
        distance=dataset_info.distance,
        energy=energy,
        delta_beta=nabu_config["phase"]["delta_beta"],
        pixel_size=dataset_info.pixel_size * 1e6,
    )
    for i in range(n_angles):
        cudapaganin.apply_filter(d_radios[i], output=d_radios[i])

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-17-dee969e8f870> in <module>
----> 1 energy = dataset_info.energy
      2 # Phase retrieval is done on each radio individually, with the sub-region specified above
      3 if (nabu_config["phase"]["method"] or "").lower() == "paganin":
      4     print("Doing phase retrieval")
      5     cudapaganin = CudaPaganinPhaseRetrieval(

NameError: name 'dataset_info' is not defined

4.3 - Logarithm

[18]:

from nabu.preproc.ccd_cuda import CudaLog

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-18-c5f2185d7784> in <module>
----> 1 from nabu.preproc.ccd_cuda import CudaLog

~/.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'

[19]:

if nabu_config["preproc"]["take_logarithm"]:
    print("Taking logarithm")
    cuda_log = CudaLog(d_radios.shape, clip_min=0.01)
    cuda_log.take_logarithm(d_radios)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-19-6efa48f64963> in <module>
----> 1 if nabu_config["preproc"]["take_logarithm"]:
      2     print("Taking logarithm")
      3     cuda_log = CudaLog(d_radios.shape, clip_min=0.01)
      4     cuda_log.take_logarithm(d_radios)

NameError: name 'nabu_config' is not defined

5 - Reconstruction

We use the filtered backprojection with nabu.reconstruction.fbp

[20]:

from nabu.reconstruction.fbp import Backprojector

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-20-02853a40d52f> in <module>
----> 1 from nabu.reconstruction.fbp import Backprojector

~/.local/lib/python3.8/site-packages/nabu/reconstruction/fbp.py in <module>
      1 import numpy as np
----> 2 import pycuda.driver as cuda
      3 from pycuda import gpuarray as garray
      4 from ..utils import updiv, get_cuda_srcfile, nextpow2, convert_index, deprecation_warning
      5 from ..cuda.utils import copy_array

ModuleNotFoundError: No module named 'pycuda'

[21]:

rec_options = conf.processing_options["reconstruction"]
B = Backprojector(
    (n_angles, n_x),
    angles=rec_options["angles"],
    rot_center=rec_options["rotation_axis_position"],
    padding_mode="edges",
    # extra_options={"use_textures": False}
)
d_recs = garray.zeros((n_z, n_x, n_x), "f")

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-21-d7e489038aee> in <module>
----> 1 rec_options = conf.processing_options["reconstruction"]
      2 B = Backprojector(
      3     (n_angles, n_x),
      4     angles=rec_options["angles"],
      5     rot_center=rec_options["rotation_axis_position"],

NameError: name 'conf' is not defined

[22]:

print("Reconstructing...", end="")
for i in range(n_z):
    B.fbp(d_radios[:, i, :], output=d_recs[i])
recs = d_recs.get()
print(" ... OK")

Reconstructing...

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-22-450afed7e212> in <module>
      1 print("Reconstructing...", end="")
----> 2 for i in range(n_z):
      3     B.fbp(d_radios[:, i, :], output=d_recs[i])
      4 recs = d_recs.get()
      5 print(" ... OK")

NameError: name 'n_z' is not defined

6 - Visualize

[23]:

%pylab nbagg

Populating the interactive namespace from numpy and matplotlib

[24]:

figure()
imshow(recs[0], cmap="gray")

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-24-2dab5e30a3e1> in <module>
      1 figure()
----> 2 imshow(recs[0], cmap="gray")

NameError: name 'recs' is not defined

[ ]: