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
print(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-fda645532af7> in <module>
1 # We can easily get information on the processing steps.
----> 2 nabu_config = conf.nabu_config
3 print(nabu_config)
4 # The same can be done with the dataset structure
5 dataset_info = conf.dataset_info
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 ChunkReader
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
sub_region = (None, None, 0, chunk_size) # start_x, end_x, start_z, end_z
chunk_reader = ChunkReader(
dataset_info.projections,
sub_region=sub_region,
convert_float=True
)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-7-a143757dd95f> in <module>
1 # Load the first 'chunk_size' lines of all the radios
----> 2 sub_region = (None, None, 0, chunk_size) # start_x, end_x, start_z, end_z
3 chunk_reader = ChunkReader(
4 dataset_info.projections,
5 sub_region=sub_region,
NameError: name 'chunk_size' is not defined
[8]:
# Load the current chunk
chunk_reader.load_files() # takes some time
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-8-8761496ebca1> in <module>
1 # Load the current chunk
----> 2 chunk_reader.load_files() # takes some time
NameError: name 'chunk_reader' is not defined
[9]:
print(chunk_reader.files_data.shape)
print(chunk_reader.files_data.dtype)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-9-e7b0d32fada6> in <module>
----> 1 print(chunk_reader.files_data.shape)
2 print(chunk_reader.files_data.dtype)
NameError: name 'chunk_reader' 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]:
radios = chunk_reader.files_data
n_angles, n_z, n_x = radios.shape
# transfer the chunk on GPU
d_radios = garray.to_gpu(radios)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-12-31ac0d7dc08c> in <module>
----> 1 radios = chunk_reader.files_data
2 n_angles, n_z, n_x = radios.shape
3 # transfer the chunk on GPU
4 d_radios = garray.to_gpu(radios)
NameError: name 'chunk_reader' 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 CudaFlatFieldDataUrls
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
<ipython-input-13-78f7f91bc045> in <module>
----> 1 from nabu.preproc.flatfield_cuda import CudaFlatFieldDataUrls
~/.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 = CudaFlatFieldDataUrls(
d_radios.shape,
dataset_info.flats,
dataset_info.darks,
radios_indices=radios_indices,
sub_region=sub_region,
convert_float=True,
)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-14-63069d9b3241> in <module>
----> 1 radios_indices = sorted(conf.dataset_info.projections.keys())
2 # Configure the `FlatField` processor
3 cuda_flatfield = CudaFlatFieldDataUrls(
4 d_radios.shape,
5 dataset_info.flats,
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"
)
d_recs = garray.zeros((n_z, n_x, n_x), "f")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-21-b4e2f99eb8be> 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(radios[:, i, :], output=d_recs[i])
recs = d_recs.get()
print(" ... OK")
Reconstructing...
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-22-e196b3e05d89> in <module>
1 print("Reconstructing...", end="")
----> 2 for i in range(n_z):
3 B.fbp(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
[ ]: